registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18
  19 /*=========================================================================
  20
  21   Program:   Insight Segmentation & Registration Toolkit
  22   Module:    $RCSfile: itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx,v $
  23   Language:  C++
  24   Date:      $Date: 2010/06/14 17:32:07 $
  25   Version:   $Revision: 1.1 $
  26
  27   Copyright (c) Insight Software Consortium. All rights reserved.
  28   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
  29
  30      This software is distributed WITHOUT ANY WARRANTY; without even
  31      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  32      PURPOSE.  See the above copyright notices for more information.
  33
  34 =========================================================================*/
  35 #ifndef __itkMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
  36 #define __itkMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
  37
  38 // First make sure that the configuration is available.
  39 // This line can be removed once the optimized versions
  40 // gets integrated into the main directories.
  41 #include "itkConfigure.h"
  42
  43 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
  44 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx"
  45 #else
  46
  47 #include "itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h"
  48 #include "itkImageRegionConstIteratorWithIndex.h"
  49
  50 namespace itk
  51 {
  52
  53 /**
  54  * Constructor
  55  */
  56 template <class TFixedImage, class TMovingImage>
  57 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
  58 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
  59 {
  60   itkDebugMacro("Constructor");
  61 }
  62
  63 /**
  64  * Get the match Measure
  65  */
  66 template <class TFixedImage, class TMovingImage>
  67 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>::MeasureType
  68 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
  69 ::GetValue( const TransformParametersType & parameters ) const
  70 {
  71
  72   itkDebugMacro("GetValue( " << parameters << " ) ");
  73
  74   FixedImageConstPointer fixedImage = this->m_FixedImage;
  75
  76   if( !fixedImage ) {
  77     itkExceptionMacro( << "Fixed image has not been assigned" );
  78   }
  79
  80   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
  81
  82
  83   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
  84
  85   typename FixedImageType::IndexType index;
  86
  87   MeasureType measure = NumericTraits< MeasureType >::Zero;
  88
  89   this->m_NumberOfPixelsCounted = 0;
  90
  91   this->SetTransformParameters( parameters );
  92
  93   while(!ti.IsAtEnd()) {
  94
  95     index = ti.GetIndex();
  96
  97     InputPointType inputPoint;
  98     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
  99
 100     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 101       ++ti;
 102       continue;
 103     }
 104
 105     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 106
 107     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 108       ++ti;
 109       continue;
 110     }
 111
 112     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 113       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 114       const RealType fixedValue   = ti.Get();
 115       this->m_NumberOfPixelsCounted++;
 116       const RealType diff = movingValue - fixedValue;
 117       measure += diff * diff;
 118     }
 119
 120     ++ti;
 121   }
 122
 123   if( !this->m_NumberOfPixelsCounted ) {
 124     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 125   } else {
 126     measure /= this->m_NumberOfPixelsCounted;
 127   }
 128
 129   return measure;
 130
 131 }
 132
 133 /**
 134  * Get the Derivative Measure
 135  */
 136 template < class TFixedImage, class TMovingImage>
 137 void
 138 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 139 ::GetDerivative( const TransformParametersType & parameters,
 140                  DerivativeType & derivative  ) const
 141 {
 142
 143   itkDebugMacro("GetDerivative( " << parameters << " ) ");
 144
 145   if( !this->GetGradientImage() ) {
 146     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 147   }
 148
 149   FixedImageConstPointer fixedImage = this->m_FixedImage;
 150
 151   if( !fixedImage ) {
 152     itkExceptionMacro( << "Fixed image has not been assigned" );
 153   }
 154
 155   const unsigned int ImageDimension = FixedImageType::ImageDimension;
 156
 157
 158   typedef  itk::ImageRegionConstIteratorWithIndex<
 159   FixedImageType> FixedIteratorType;
 160
 161   typedef  itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
 162
 163
 164   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 165
 166   typename FixedImageType::IndexType index;
 167
 168   this->m_NumberOfPixelsCounted = 0;
 169
 170   this->SetTransformParameters( parameters );
 171
 172   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 173   derivative = DerivativeType( ParametersDimension );
 174   derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 175
 176   ti.GoToBegin();
 177
 178   while(!ti.IsAtEnd()) {
 179
 180     index = ti.GetIndex();
 181
 182     InputPointType inputPoint;
 183     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 184
 185     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 186       ++ti;
 187       continue;
 188     }
 189
 190     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 191
 192     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 193       ++ti;
 194       continue;
 195     }
 196
 197     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 198       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 199
 200       const TransformJacobianType & jacobian =
 201         this->m_Transform->GetJacobian( inputPoint );
 202
 203
 204       const RealType fixedValue     = ti.Value();
 205       this->m_NumberOfPixelsCounted++;
 206       const RealType diff = movingValue - fixedValue;
 207
 208       // Get the gradient by NearestNeighboorInterpolation:
 209       // which is equivalent to round up the point components.
 210       typedef typename OutputPointType::CoordRepType CoordRepType;
 211       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 212       MovingImageContinuousIndexType;
 213
 214       MovingImageContinuousIndexType tempIndex;
 215       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 216
 217       typename MovingImageType::IndexType mappedIndex;
 218       mappedIndex.CopyWithRound( tempIndex );
 219
 220       const GradientPixelType gradient =
 221         this->GetGradientImage()->GetPixel( mappedIndex );
 222
 223       for(unsigned int par=0; par<ParametersDimension; par++) {
 224         RealType sum = NumericTraits< RealType >::Zero;
 225         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 226           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 227         }
 228         derivative[par] += sum;
 229       }
 230     }
 231
 232     ++ti;
 233   }
 234
 235   if( !this->m_NumberOfPixelsCounted ) {
 236     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 237   } else {
 238     for(unsigned int i=0; i<ParametersDimension; i++) {
 239       derivative[i] /= this->m_NumberOfPixelsCounted;
 240     }
 241   }
 242
 243 }
 244
 245
 246 /*
 247  * Get both the match Measure and theDerivative Measure
 248  */
 249 template <class TFixedImage, class TMovingImage>
 250 void
 251 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 252 ::GetValueAndDerivative(const TransformParametersType & parameters,
 253                         MeasureType & value, DerivativeType  & derivative) const
 254 {
 255
 256   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
 257
 258   if( !this->GetGradientImage() ) {
 259     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 260   }
 261
 262   FixedImageConstPointer fixedImage = this->m_FixedImage;
 263
 264   if( !fixedImage ) {
 265     itkExceptionMacro( << "Fixed image has not been assigned" );
 266   }
 267
 268   const unsigned int ImageDimension = FixedImageType::ImageDimension;
 269
 270   typedef  itk::ImageRegionConstIteratorWithIndex<
 271   FixedImageType> FixedIteratorType;
 272
 273   typedef  itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
 274
 275
 276   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 277
 278   typename FixedImageType::IndexType index;
 279
 280   MeasureType measure = NumericTraits< MeasureType >::Zero;
 281
 282   this->m_NumberOfPixelsCounted = 0;
 283
 284   this->SetTransformParameters( parameters );
 285
 286   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 287   derivative = DerivativeType( ParametersDimension );
 288   derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 289
 290   ti.GoToBegin();
 291
 292   while(!ti.IsAtEnd()) {
 293
 294     index = ti.GetIndex();
 295
 296     InputPointType inputPoint;
 297     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 298
 299     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 300       ++ti;
 301       continue;
 302     }
 303
 304     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 305
 306     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 307       ++ti;
 308       continue;
 309     }
 310
 311     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 312       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 313
 314       const TransformJacobianType & jacobian =
 315         this->m_Transform->GetJacobian( inputPoint );
 316
 317
 318       const RealType fixedValue     = ti.Value();
 319       this->m_NumberOfPixelsCounted++;
 320
 321       const RealType diff = movingValue - fixedValue;
 322
 323       measure += diff * diff;
 324
 325       // Get the gradient by NearestNeighboorInterpolation:
 326       // which is equivalent to round up the point components.
 327       typedef typename OutputPointType::CoordRepType CoordRepType;
 328       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 329       MovingImageContinuousIndexType;
 330
 331       MovingImageContinuousIndexType tempIndex;
 332       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 333
 334       typename MovingImageType::IndexType mappedIndex;
 335       mappedIndex.CopyWithRound( tempIndex );
 336
 337       const GradientPixelType gradient =
 338         this->GetGradientImage()->GetPixel( mappedIndex );
 339
 340       for(unsigned int par=0; par<ParametersDimension; par++) {
 341         RealType sum = NumericTraits< RealType >::Zero;
 342         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 343           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 344         }
 345         derivative[par] += sum;
 346       }
 347     }
 348
 349     ++ti;
 350   }
 351
 352   if( !this->m_NumberOfPixelsCounted ) {
 353     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 354   } else {
 355     for(unsigned int i=0; i<ParametersDimension; i++) {
 356       derivative[i] /= this->m_NumberOfPixelsCounted;
 357     }
 358     measure /= this->m_NumberOfPixelsCounted;
 359   }
 360
 361   value = measure;
 362
 363 }
 364
 365 } // end namespace itk
 366
 367
 368 #endif
 369
 370 #endif