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 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx"
  44
  45 #include "itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h"
  46 #include "itkImageRegionConstIteratorWithIndex.h"
  47
  48 namespace itk
  49 {
  50
  51 /**
  52  * Constructor
  53  */
  54 template <class TFixedImage, class TMovingImage>
  55 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
  56 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
  57 {
  58   itkDebugMacro("Constructor");
  59 }
  60
  61 /**
  62  * Get the match Measure
  63  */
  64 template <class TFixedImage, class TMovingImage>
  65 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>::MeasureType
  66 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
  67 ::GetValue( const TransformParametersType & parameters ) const
  68 {
  69
  70   itkDebugMacro("GetValue( " << parameters << " ) ");
  71
  72   FixedImageConstPointer fixedImage = this->m_FixedImage;
  73
  74   if( !fixedImage ) {
  75     itkExceptionMacro( << "Fixed image has not been assigned" );
  76   }
  77
  78   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
  79
  80
  81   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
  82
  83   typename FixedImageType::IndexType index;
  84
  85   MeasureType measure = NumericTraits< MeasureType >::Zero;
  86
  87   this->m_NumberOfPixelsCounted = 0;
  88
  89   this->SetTransformParameters( parameters );
  90
  91   while(!ti.IsAtEnd()) {
  92
  93     index = ti.GetIndex();
  94
  95     InputPointType inputPoint;
  96     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
  97
  98     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
  99       ++ti;
 100       continue;
 101     }
 102
 103     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 104
 105     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 106       ++ti;
 107       continue;
 108     }
 109
 110     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 111       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 112       const RealType fixedValue   = ti.Get();
 113       this->m_NumberOfPixelsCounted++;
 114       const RealType diff = movingValue - fixedValue;
 115       measure += diff * diff;
 116     }
 117
 118     ++ti;
 119   }
 120
 121   if( !this->m_NumberOfPixelsCounted ) {
 122     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 123   } else {
 124     measure /= this->m_NumberOfPixelsCounted;
 125   }
 126
 127   return measure;
 128
 129 }
 130
 131 /**
 132  * Get the Derivative Measure
 133  */
 134 template < class TFixedImage, class TMovingImage>
 135 void
 136 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 137 ::GetDerivative( const TransformParametersType & parameters,
 138                  DerivativeType & derivative  ) const
 139 {
 140
 141   itkDebugMacro("GetDerivative( " << parameters << " ) ");
 142
 143   if( !this->GetGradientImage() ) {
 144     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 145   }
 146
 147   FixedImageConstPointer fixedImage = this->m_FixedImage;
 148
 149   if( !fixedImage ) {
 150     itkExceptionMacro( << "Fixed image has not been assigned" );
 151   }
 152
 153   const unsigned int ImageDimension = FixedImageType::ImageDimension;
 154
 155
 156   typedef  itk::ImageRegionConstIteratorWithIndex<
 157   FixedImageType> FixedIteratorType;
 158
 159   typedef  itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
 160
 161
 162   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 163
 164   typename FixedImageType::IndexType index;
 165
 166   this->m_NumberOfPixelsCounted = 0;
 167
 168   this->SetTransformParameters( parameters );
 169
 170   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 171   derivative = DerivativeType( ParametersDimension );
 172   derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
 173
 174   ti.GoToBegin();
 175
 176   while(!ti.IsAtEnd()) {
 177
 178     index = ti.GetIndex();
 179
 180     InputPointType inputPoint;
 181     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 182
 183     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 184       ++ti;
 185       continue;
 186     }
 187
 188     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 189
 190     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 191       ++ti;
 192       continue;
 193     }
 194
 195     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 196       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 197
 198       TransformJacobianType jacobian;
 199       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
 200
 201       const RealType fixedValue     = ti.Value();
 202       this->m_NumberOfPixelsCounted++;
 203       const RealType diff = movingValue - fixedValue;
 204
 205       // Get the gradient by NearestNeighboorInterpolation:
 206       // which is equivalent to round up the point components.
 207       typedef typename OutputPointType::CoordRepType CoordRepType;
 208       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 209       MovingImageContinuousIndexType;
 210
 211       MovingImageContinuousIndexType tempIndex;
 212       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 213
 214       typename MovingImageType::IndexType mappedIndex;
 215       mappedIndex.CopyWithRound( tempIndex );
 216
 217       const GradientPixelType gradient =
 218         this->GetGradientImage()->GetPixel( mappedIndex );
 219
 220       for(unsigned int par=0; par<ParametersDimension; par++) {
 221         RealType sum = NumericTraits< RealType >::Zero;
 222         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 223           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 224         }
 225         derivative[par] += sum;
 226       }
 227     }
 228
 229     ++ti;
 230   }
 231
 232   if( !this->m_NumberOfPixelsCounted ) {
 233     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 234   } else {
 235     for(unsigned int i=0; i<ParametersDimension; i++) {
 236       derivative[i] /= this->m_NumberOfPixelsCounted;
 237     }
 238   }
 239
 240 }
 241
 242
 243 /*
 244  * Get both the match Measure and theDerivative Measure
 245  */
 246 template <class TFixedImage, class TMovingImage>
 247 void
 248 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 249 ::GetValueAndDerivative(const TransformParametersType & parameters,
 250                         MeasureType & value, DerivativeType  & derivative) const
 251 {
 252
 253   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
 254
 255   if( !this->GetGradientImage() ) {
 256     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 257   }
 258
 259   FixedImageConstPointer fixedImage = this->m_FixedImage;
 260
 261   if( !fixedImage ) {
 262     itkExceptionMacro( << "Fixed image has not been assigned" );
 263   }
 264
 265   const unsigned int ImageDimension = FixedImageType::ImageDimension;
 266
 267   typedef  itk::ImageRegionConstIteratorWithIndex<
 268   FixedImageType> FixedIteratorType;
 269
 270   typedef  itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
 271
 272
 273   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 274
 275   typename FixedImageType::IndexType index;
 276
 277   MeasureType measure = NumericTraits< MeasureType >::Zero;
 278
 279   this->m_NumberOfPixelsCounted = 0;
 280
 281   this->SetTransformParameters( parameters );
 282
 283   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 284   derivative = DerivativeType( ParametersDimension );
 285   derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
 286
 287   ti.GoToBegin();
 288
 289   while(!ti.IsAtEnd()) {
 290
 291     index = ti.GetIndex();
 292
 293     InputPointType inputPoint;
 294     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 295
 296     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 297       ++ti;
 298       continue;
 299     }
 300
 301     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 302
 303     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 304       ++ti;
 305       continue;
 306     }
 307
 308     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 309       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 310
 311       TransformJacobianType jacobian;
 312       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
 313
 314       const RealType fixedValue     = ti.Value();
 315       this->m_NumberOfPixelsCounted++;
 316
 317       const RealType diff = movingValue - fixedValue;
 318
 319       measure += diff * diff;
 320
 321       // Get the gradient by NearestNeighboorInterpolation:
 322       // which is equivalent to round up the point components.
 323       typedef typename OutputPointType::CoordRepType CoordRepType;
 324       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 325       MovingImageContinuousIndexType;
 326
 327       MovingImageContinuousIndexType tempIndex;
 328       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 329
 330       typename MovingImageType::IndexType mappedIndex;
 331       mappedIndex.CopyWithRound( tempIndex );
 332
 333       const GradientPixelType gradient =
 334         this->GetGradientImage()->GetPixel( mappedIndex );
 335
 336       for(unsigned int par=0; par<ParametersDimension; par++) {
 337         RealType sum = NumericTraits< RealType >::Zero;
 338         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 339           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 340         }
 341         derivative[par] += sum;
 342       }
 343     }
 344
 345     ++ti;
 346   }
 347
 348   if( !this->m_NumberOfPixelsCounted ) {
 349     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 350   } else {
 351     for(unsigned int i=0; i<ParametersDimension; i++) {
 352       derivative[i] /= this->m_NumberOfPixelsCounted;
 353     }
 354     measure /= this->m_NumberOfPixelsCounted;
 355   }
 356
 357   value = measure;
 358
 359 }
 360
 361 } // end namespace itk
 362
 363
 364 #endif
 365
 366 #endif