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 #if ITK_VERSION_MAJOR >= 4
 201       TransformJacobianType jacobian;
 202       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
 203 #else
 204       const TransformJacobianType & jacobian =
 205         this->m_Transform->GetJacobian( inputPoint );
 206 #endif
 207
 208       const RealType fixedValue     = ti.Value();
 209       this->m_NumberOfPixelsCounted++;
 210       const RealType diff = movingValue - fixedValue;
 211
 212       // Get the gradient by NearestNeighboorInterpolation:
 213       // which is equivalent to round up the point components.
 214       typedef typename OutputPointType::CoordRepType CoordRepType;
 215       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 216       MovingImageContinuousIndexType;
 217
 218       MovingImageContinuousIndexType tempIndex;
 219       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 220
 221       typename MovingImageType::IndexType mappedIndex;
 222       mappedIndex.CopyWithRound( tempIndex );
 223
 224       const GradientPixelType gradient =
 225         this->GetGradientImage()->GetPixel( mappedIndex );
 226
 227       for(unsigned int par=0; par<ParametersDimension; par++) {
 228         RealType sum = NumericTraits< RealType >::Zero;
 229         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 230           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 231         }
 232         derivative[par] += sum;
 233       }
 234     }
 235
 236     ++ti;
 237   }
 238
 239   if( !this->m_NumberOfPixelsCounted ) {
 240     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 241   } else {
 242     for(unsigned int i=0; i<ParametersDimension; i++) {
 243       derivative[i] /= this->m_NumberOfPixelsCounted;
 244     }
 245   }
 246
 247 }
 248
 249
 250 /*
 251  * Get both the match Measure and theDerivative Measure
 252  */
 253 template <class TFixedImage, class TMovingImage>
 254 void
 255 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 256 ::GetValueAndDerivative(const TransformParametersType & parameters,
 257                         MeasureType & value, DerivativeType  & derivative) const
 258 {
 259
 260   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
 261
 262   if( !this->GetGradientImage() ) {
 263     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 264   }
 265
 266   FixedImageConstPointer fixedImage = this->m_FixedImage;
 267
 268   if( !fixedImage ) {
 269     itkExceptionMacro( << "Fixed image has not been assigned" );
 270   }
 271
 272   const unsigned int ImageDimension = FixedImageType::ImageDimension;
 273
 274   typedef  itk::ImageRegionConstIteratorWithIndex<
 275   FixedImageType> FixedIteratorType;
 276
 277   typedef  itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
 278
 279
 280   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 281
 282   typename FixedImageType::IndexType index;
 283
 284   MeasureType measure = NumericTraits< MeasureType >::Zero;
 285
 286   this->m_NumberOfPixelsCounted = 0;
 287
 288   this->SetTransformParameters( parameters );
 289
 290   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 291   derivative = DerivativeType( ParametersDimension );
 292   derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 293
 294   ti.GoToBegin();
 295
 296   while(!ti.IsAtEnd()) {
 297
 298     index = ti.GetIndex();
 299
 300     InputPointType inputPoint;
 301     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 302
 303     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 304       ++ti;
 305       continue;
 306     }
 307
 308     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 309
 310     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 311       ++ti;
 312       continue;
 313     }
 314
 315     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 316       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 317
 318 #if ITK_VERSION_MAJOR >= 4
 319       TransformJacobianType jacobian;
 320       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
 321 #else
 322       const TransformJacobianType & jacobian =
 323         this->m_Transform->GetJacobian( inputPoint );
 324 #endif
 325
 326       const RealType fixedValue     = ti.Value();
 327       this->m_NumberOfPixelsCounted++;
 328
 329       const RealType diff = movingValue - fixedValue;
 330
 331       measure += diff * diff;
 332
 333       // Get the gradient by NearestNeighboorInterpolation:
 334       // which is equivalent to round up the point components.
 335       typedef typename OutputPointType::CoordRepType CoordRepType;
 336       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 337       MovingImageContinuousIndexType;
 338
 339       MovingImageContinuousIndexType tempIndex;
 340       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 341
 342       typename MovingImageType::IndexType mappedIndex;
 343       mappedIndex.CopyWithRound( tempIndex );
 344
 345       const GradientPixelType gradient =
 346         this->GetGradientImage()->GetPixel( mappedIndex );
 347
 348       for(unsigned int par=0; par<ParametersDimension; par++) {
 349         RealType sum = NumericTraits< RealType >::Zero;
 350         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 351           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 352         }
 353         derivative[par] += sum;
 354       }
 355     }
 356
 357     ++ti;
 358   }
 359
 360   if( !this->m_NumberOfPixelsCounted ) {
 361     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 362   } else {
 363     for(unsigned int i=0; i<ParametersDimension; i++) {
 364       derivative[i] /= this->m_NumberOfPixelsCounted;
 365     }
 366     measure /= this->m_NumberOfPixelsCounted;
 367   }
 368
 369   value = measure;
 370
 371 }
 372
 373 } // end namespace itk
 374
 375
 376 #endif
 377
 378 #endif