00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef __freDifferenceHistogramImageToImageMetric_txx
00023 #define __freDifferenceHistogramImageToImageMetric_txx
00024
00025 #include "itkArray.h"
00026 #include "itkNumericTraits.h"
00027 #include "itkImageRegionConstIterator.h"
00028 #include "itkImageRegionConstIteratorWithIndex.h"
00029 #include "freDifferenceHistogramImageToImageMetric.h"
00030
00031 namespace itk
00032 {
00033 template <class TFixedImage, class TMovingImage, class TDifImage>
00034 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00035 ::DifferenceHistogramImageToImageMetric()
00036 {
00037 itkDebugMacro("Constructor");
00038
00039 m_HistogramSize.Fill(512);
00040 m_UseMask = false;
00041 m_DerivativeStepLength = 0.1;
00042 m_UpperBoundIncreaseFactor = 0.001;
00043 m_MaskedMaximum = 20;
00044 m_MaskedMinimum = NumericTraits<FixedImagePixelType>::Zero;
00045 m_Histogram = HistogramType::New();
00046 }
00047
00048 template <class TFixedImage, class TMovingImage, class TDifImage>
00049 void
00050 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00051 ::Initialize() throw (ExceptionObject)
00052 {
00053 Superclass::Initialize();
00054
00055
00056 FixedImageConstPointer pFixedImage = this->m_FixedImage;
00057 ImageRegionConstIterator<FixedImageType> fiIt(pFixedImage,
00058 pFixedImage->GetBufferedRegion());
00059 fiIt.GoToBegin();
00060 FixedImagePixelType minFixed = fiIt.Value();
00061 FixedImagePixelType maxFixed = fiIt.Value();
00062 ++fiIt;
00063 while ( !fiIt.IsAtEnd() )
00064 {
00065 FixedImagePixelType value = fiIt.Value();
00066
00067 if (value < minFixed)
00068 {
00069 minFixed = value;
00070 }
00071 else if (value > maxFixed)
00072 {
00073 maxFixed = value;
00074 }
00075
00076 ++fiIt;
00077 }
00078
00079
00080 MovingImageConstPointer pMovingImage = this->m_MovingImage;
00081 ImageRegionConstIterator<MovingImageType> miIt(pMovingImage,
00082 pMovingImage->GetBufferedRegion());
00083 miIt.GoToBegin();
00084 MovingImagePixelType minMoving = miIt.Value();
00085 MovingImagePixelType maxMoving = miIt.Value();
00086 ++miIt;
00087 while ( !miIt.IsAtEnd() )
00088 {
00089 MovingImagePixelType value = miIt.Value();
00090
00091 if (value < minMoving)
00092 {
00093 minMoving = value;
00094 }
00095 else if (value > maxMoving)
00096 {
00097 maxMoving = value;
00098 }
00099 ++miIt;
00100 }
00101
00102
00103
00104
00105 m_LowerBound[0] = minMoving - maxFixed;
00106 m_UpperBound[0] = maxMoving - minFixed;
00107 m_UpperBound[0] = m_UpperBound[0] + (m_UpperBound[0] - m_LowerBound[0] ) * m_UpperBoundIncreaseFactor;
00108 }
00109
00110 template <class TFixedImage, class TMovingImage, class TDifImage>
00111 typename DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>::MeasureType
00112 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00113 ::GetValue(const TransformParametersType& parameters) const
00114 {
00115 itkDebugMacro("GetValue( " << parameters << " ) ");
00116
00117 this->ComputeHistogram(parameters, *m_Histogram);
00118
00119 return this->EvaluateMeasure(*m_Histogram);
00120 }
00121
00122 template <class TFixedImage, class TMovingImage, class TDifImage>
00123 void
00124 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00125 ::GetDerivative(const TransformParametersType& parameters,
00126 DerivativeType& derivative) const
00127 {
00128 itkDebugMacro("GetDerivative( " << parameters << " ) ");
00129
00130
00131 const unsigned int ParametersDimension = this->GetNumberOfParameters();
00132 derivative = DerivativeType(ParametersDimension);
00133 derivative.Fill(NumericTraits<ITK_TYPENAME
00134 DerivativeType::ValueType>::Zero);
00135
00136 typename HistogramType::Pointer pHistogram = HistogramType::New();
00137 this->ComputeHistogram(parameters, *pHistogram);
00138
00139 for (unsigned int i = 0; i < ParametersDimension; i++)
00140 {
00141 typename HistogramType::Pointer pHistogram2 = HistogramType::New();
00142 this->CopyHistogram(*pHistogram2, *pHistogram);
00143
00144 TransformParametersType newParameters = parameters;
00145 newParameters[i] -= m_DerivativeStepLength;
00146 this->ComputeHistogram(newParameters, *pHistogram2);
00147
00148 MeasureType e0 = EvaluateMeasure(*pHistogram2);
00149
00150 pHistogram2 = HistogramType::New();
00151 this->CopyHistogram(*pHistogram2, *pHistogram);
00152
00153 newParameters = parameters;
00154 newParameters[i] += m_DerivativeStepLength;
00155 this->ComputeHistogram(newParameters, *pHistogram2);
00156
00157 MeasureType e1 = EvaluateMeasure(*pHistogram2);
00158
00159 derivative[i] = (e1 - e0)/(2*m_DerivativeStepLength);
00160 }
00161 }
00162
00163 template <class TFixedImage, class TMovingImage, class TDifImage>
00164 void
00165 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00166 ::GetValueAndDerivative(const TransformParametersType& parameters,
00167 MeasureType& value,
00168 DerivativeType& derivative) const
00169 {
00170 value = GetValue(parameters);
00171 this->GetDerivative(parameters, derivative);
00172 }
00173
00174 template <class TFixedImage, class TMovingImage, class TDifImage>
00175 void
00176 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00177 ::ComputeHistogram(TransformParametersType const& parameters,
00178 HistogramType& histogram) const
00179 {
00180 FixedImageConstPointer fixedImage = this->m_FixedImage;
00181
00182 if(!fixedImage)
00183 {
00184 itkExceptionMacro(<< "Fixed image has not been assigned");
00185 }
00186
00187 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType>
00188 FixedIteratorType;
00189
00190 typename FixedImageType::IndexType index;
00191 typename FixedImageType::RegionType fixedRegion;
00192
00193 fixedRegion = this->GetFixedImageRegion();
00194 FixedIteratorType ti( fixedImage, fixedRegion );
00195
00196 this->m_NumberOfPixelsCounted = 0;
00197 this->SetTransformParameters( parameters );
00198
00199 histogram.Initialize( m_HistogramSize, m_LowerBound, m_UpperBound );
00200
00201 ti.GoToBegin();
00202 while ( !ti.IsAtEnd() )
00203 {
00204 index = ti.GetIndex();
00205 const RealType fixedValue = ti.Get();
00206
00207 if (m_UseMask && IsMaskedValue(fixedValue))
00208 {
00209 ++ti;
00210 continue;
00211 }
00212
00213 typename Superclass::InputPointType inputPoint;
00214 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
00215
00216 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
00217
00218 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
00219 {
00220 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
00221
00222 if (m_UseMask && IsMaskedValue(movingValue))
00223 {
00224 ++ti;
00225 continue;
00226 }
00227
00228 this->m_NumberOfPixelsCounted++;
00229
00230 typename HistogramType::MeasurementVectorType sample;
00231 sample[0] = movingValue - fixedValue;
00232 histogram.IncreaseFrequency(sample, 1);
00233 }
00234 ++ti;
00235 }
00236
00237 if (this->m_NumberOfPixelsCounted == 0)
00238 {
00239 itkExceptionMacro(<< "All the points mapped to outside of the moving \
00240 image");
00241 }
00242 }
00243
00244 template <class TFixedImage, class TMovingImage, class TDifImage>
00245 void
00246 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00247 ::CopyHistogram(HistogramType& target, HistogramType& source) const
00248 {
00249
00250 typename HistogramType::MeasurementVectorType min, max;
00251 typename HistogramType::SizeType size = source.GetSize();
00252
00253 for (unsigned int i = 0; i < min.Size(); i++)
00254 {
00255 min[i] = source.GetBinMin(i, 0);
00256 }
00257
00258 for (unsigned int i = 0; i < max.Size(); i++)
00259 {
00260 max[i] = source.GetBinMax(i, size[i] - 1);
00261 }
00262
00263 target.Initialize(size, min, max);
00264
00265
00266 typename HistogramType::Iterator sourceIt = source.Begin();
00267 typename HistogramType::Iterator sourceEnd = source.End();
00268 typename HistogramType::Iterator targetIt = target.Begin();
00269 typename HistogramType::Iterator targetEnd = target.End();
00270
00271 while (sourceIt != sourceEnd && targetIt != targetEnd)
00272 {
00273 typename HistogramType::FrequencyType freq = sourceIt.GetFrequency();
00274
00275 if (freq > 0)
00276 {
00277 targetIt.SetFrequency(freq);
00278 }
00279
00280 ++sourceIt;
00281 ++targetIt;
00282 }
00283 }
00284
00285 template <class TFixedImage, class TMovingImage, class TDifImage>
00286 void
00287 DifferenceHistogramImageToImageMetric<TFixedImage,TMovingImage,TDifImage>
00288 ::PrintSelf(std::ostream& os, Indent indent) const
00289 {
00290 Superclass::PrintSelf(os,indent);
00291 os << indent << "Use mask value?: " << m_UseMask << std::endl;
00292 os << indent << "Derivative step length: " << m_DerivativeStepLength
00293 << std::endl;
00294 os << indent << "Histogram size: ";
00295 os << m_HistogramSize << std::endl;
00296 os << indent << "Histogram upper bound increase factor: ";
00297 os << m_UpperBoundIncreaseFactor << std::endl;
00298 os << indent << "Histogram computed by GetValue(): ";
00299 os << m_Histogram.GetPointer() << std::endl;
00300 }
00301 }
00302
00303 #endif // __freDifferenceHistogramImageToImageMetric_txx
