registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.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: itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.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 __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
  36 #define __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
  37
  38 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.h"
  39 #include "itkCovariantVector.h"
  40 #include "itkImageRandomConstIteratorWithIndex.h"
  41 #include "itkImageRegionConstIterator.h"
  42 #include "itkImageRegionIterator.h"
  43 #include "itkImageIterator.h"
  44 #include "vnl/vnl_math.h"
  45
  46 namespace itk
  47 {
  48
  49 /**
  50  * Constructor
  51  */
  52 template < class TFixedImage, class TMovingImage >
  53 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
  54 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
  55 {
  56   this->SetComputeGradient(true);
  57
  58   m_ThreaderMSE = NULL;
  59   m_ThreaderMSEDerivatives = NULL;
  60   this->m_WithinThreadPreProcess = false;
  61   this->m_WithinThreadPostProcess = false;
  62
  63   //  For backward compatibility, the default behavior is to use all the pixels
  64   //  in the fixed image.
  65   this->UseAllPixelsOn();
  66 }
  67
  68 template < class TFixedImage, class TMovingImage >
  69 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
  70 ::~MeanSquaresImageToImageMetricFor3DBLUTFFD()
  71 {
  72   if(m_ThreaderMSE != NULL) {
  73     delete [] m_ThreaderMSE;
  74   }
  75   m_ThreaderMSE = NULL;
  76
  77   if(m_ThreaderMSEDerivatives != NULL) {
  78     delete [] m_ThreaderMSEDerivatives;
  79   }
  80   m_ThreaderMSEDerivatives = NULL;
  81 }
  82
  83 /**
  84  * Print out internal information about this class
  85  */
  86 template < class TFixedImage, class TMovingImage  >
  87 void
  88 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
  89 ::PrintSelf(std::ostream& os, Indent indent) const
  90 {
  91
  92   Superclass::PrintSelf(os, indent);
  93
  94 }
  95
  96
  97 /**
  98  * Initialize
  99  */
 100 template <class TFixedImage, class TMovingImage>
 101 void
 102 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 103 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
 104 ::Initialize(void)
 105 #else
 106 ::Initialize(void) throw ( ExceptionObject )
 107 #endif
 108 {
 109
 110   this->Superclass::Initialize();
 111   this->Superclass::MultiThreadingInitialize();
 112
 113   if(m_ThreaderMSE != NULL) {
 114     delete [] m_ThreaderMSE;
 115   }
 116 #if ITK_VERSION_MAJOR <= 4
 117   m_ThreaderMSE = new double[this->m_NumberOfThreads];
 118 #else
 119   m_ThreaderMSE = new double[this->m_NumberOfWorkUnits];
 120 #endif
 121
 122   if(m_ThreaderMSEDerivatives != NULL) {
 123     delete [] m_ThreaderMSEDerivatives;
 124   }
 125 #if ITK_VERSION_MAJOR <= 4
 126   m_ThreaderMSEDerivatives = new DerivativeType[this->m_NumberOfThreads];
 127   for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
 128 #else
 129   m_ThreaderMSEDerivatives = new DerivativeType[this->m_NumberOfWorkUnits];
 130   for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
 131 #endif
 132     m_ThreaderMSEDerivatives[threadID].SetSize( this->m_NumberOfParameters );
 133   }
 134 }
 135
 136 template < class TFixedImage, class TMovingImage  >
 137 inline bool
 138 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 139 ::GetValueThreadProcessSample( unsigned int threadID,
 140                                unsigned long fixedImageSample,
 141                                const MovingImagePointType & itkNotUsed(mappedPoint),
 142                                double movingImageValue) const
 143 {
 144   double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
 145
 146   m_ThreaderMSE[threadID] += diff*diff;
 147
 148   return true;
 149 }
 150
 151 template < class TFixedImage, class TMovingImage  >
 152 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 153 ::MeasureType
 154 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 155 ::GetValue( const ParametersType & parameters ) const
 156 {
 157   itkDebugMacro("GetValue( " << parameters << " ) ");
 158
 159   if( !this->m_FixedImage ) {
 160     itkExceptionMacro( << "Fixed image has not been assigned" );
 161   }
 162
 163 #if ITK_VERSION_MAJOR <= 4
 164   memset( m_ThreaderMSE, 0, this->m_NumberOfThreads * sizeof(MeasureType) );
 165 #else
 166   memset( m_ThreaderMSE, 0, this->m_NumberOfWorkUnits * sizeof(MeasureType) );
 167 #endif
 168
 169   // Set up the parameters in the transform
 170   this->m_Transform->SetParameters( parameters );
 171
 172   // MUST BE CALLED TO INITIATE PROCESSING
 173   this->GetValueMultiThreadedInitiate();
 174
 175   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
 176                  << this->m_NumberOfPixelsCounted << " / "
 177                  << this->m_NumberOfFixedImageSamples
 178                  << std::endl );
 179
 180   if( this->m_NumberOfPixelsCounted <
 181       this->m_NumberOfFixedImageSamples / 4 ) {
 182     itkExceptionMacro( "Too many samples map outside moving image buffer: "
 183                        << this->m_NumberOfPixelsCounted << " / "
 184                        << this->m_NumberOfFixedImageSamples
 185                        << std::endl );
 186   }
 187
 188   double mse = m_ThreaderMSE[0];
 189 #if ITK_VERSION_MAJOR <= 4
 190   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
 191 #else
 192   for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
 193 #endif
 194     mse += m_ThreaderMSE[t];
 195   }
 196   mse /= this->m_NumberOfPixelsCounted;
 197
 198   return mse;
 199 }
 200
 201
 202 template < class TFixedImage, class TMovingImage  >
 203 inline bool
 204 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 205 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
 206     unsigned long fixedImageSample,
 207     const MovingImagePointType & itkNotUsed(mappedPoint),
 208     double movingImageValue,
 209     const ImageDerivativesType &
 210     movingImageGradientValue ) const
 211 {
 212   double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
 213
 214   m_ThreaderMSE[threadID] += diff*diff;
 215
 216   //JV
 217   //FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
 218
 219   // Need to use one of the threader transforms if we're
 220   // not in thread 0.
 221   //
 222   // Use a raw pointer here to avoid the overhead of smart pointers.
 223   // For instance, Register and UnRegister have mutex locks around
 224   // the reference counts.
 225   TransformType* transform;
 226
 227   if (threadID > 0) {
 228     transform = this->m_ThreaderTransform[threadID - 1];
 229   } else {
 230     transform = this->m_Transform;
 231   }
 232
 233   // Jacobian should be evaluated at the unmapped (fixed image) point.
 234   TransformJacobianType jacobian;
 235   transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[fixedImageSample].point, jacobian);
 236   //double sum;
 237   unsigned int par, dim;
 238   for( par=0; par<this->m_NumberOfParameters; par+=3) {
 239     //     double sum = 0.0;
 240     //     for(unsigned int dim=0; dim<MovingImageDimension; dim++)
 241     //       {
 242     //       sum += 2.0 * diff * jacobian( dim, par ) * movingImageGradientValue[dim];
 243     //       }
 244     //     m_ThreaderMSEDerivatives[threadID][par] += sum;
 245
 246     // JV only for 3D BLUT FFD
 247     if (jacobian( 0, par ) )
 248       for( dim=0; dim<3; dim++)
 249         m_ThreaderMSEDerivatives[threadID][par+dim  ]  += 2.0 * diff* jacobian( dim, par+dim   ) * movingImageGradientValue[dim];
 250
 251   }
 252   return true;
 253 }
 254
 255 /**
 256  * Get the both Value and Derivative Measure
 257  */
 258 template < class TFixedImage, class TMovingImage  >
 259 void
 260 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 261 ::GetValueAndDerivative( const ParametersType & parameters,
 262                          MeasureType & value,
 263                          DerivativeType & derivative) const
 264 {
 265
 266   if( !this->m_FixedImage ) {
 267     itkExceptionMacro( << "Fixed image has not been assigned" );
 268   }
 269
 270   // Set up the parameters in the transform
 271   this->m_Transform->SetParameters( parameters );
 272
 273   // Reset the joint pdfs to zero
 274 #if ITK_VERSION_MAJOR <= 4
 275   memset( m_ThreaderMSE, 0, this->m_NumberOfThreads * sizeof(MeasureType) );
 276 #else
 277   memset( m_ThreaderMSE, 0, this->m_NumberOfWorkUnits * sizeof(MeasureType) );
 278 #endif
 279
 280   // Set output values to zero
 281   if(derivative.GetSize() != this->m_NumberOfParameters) {
 282     derivative = DerivativeType( this->m_NumberOfParameters );
 283   }
 284   memset( derivative.data_block(),
 285           0,
 286           this->m_NumberOfParameters * sizeof(double) );
 287
 288 #if ITK_VERSION_MAJOR <= 4
 289   for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
 290 #else
 291   for( unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits; threadID++ ) {
 292 #endif
 293     memset( m_ThreaderMSEDerivatives[threadID].data_block(),
 294             0,
 295             this->m_NumberOfParameters * sizeof(double) );
 296   }
 297
 298   // MUST BE CALLED TO INITIATE PROCESSING
 299   this->GetValueAndDerivativeMultiThreadedInitiate();
 300
 301   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
 302                  << this->m_NumberOfPixelsCounted << " / "
 303                  << this->m_NumberOfFixedImageSamples
 304                  << std::endl );
 305
 306   if( this->m_NumberOfPixelsCounted <
 307       this->m_NumberOfFixedImageSamples / 4 ) {
 308     itkExceptionMacro( "Too many samples map outside moving image buffer: "
 309                        << this->m_NumberOfPixelsCounted << " / "
 310                        << this->m_NumberOfFixedImageSamples
 311                        << std::endl );
 312   }
 313
 314   value = 0;
 315 #if ITK_VERSION_MAJOR <= 4
 316   for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
 317 #else
 318   for(unsigned int t=0; t<this->m_NumberOfWorkUnits; t++) {
 319 #endif
 320     value += m_ThreaderMSE[t];
 321     for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
 322         parameter++) {
 323       derivative[parameter] += m_ThreaderMSEDerivatives[t][parameter];
 324     }
 325   }
 326
 327   value /= this->m_NumberOfPixelsCounted;
 328   for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
 329       parameter++) {
 330     derivative[parameter] /= this->m_NumberOfPixelsCounted;
 331     // JV
 332     //DD(parameter<<"\t"<<derivative[parameter]);
 333   }
 334
 335 }
 336
 337
 338 /**
 339  * Get the match measure derivative
 340  */
 341 template < class TFixedImage, class TMovingImage  >
 342 void
 343 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 344 ::GetDerivative( const ParametersType & parameters,
 345                  DerivativeType & derivative ) const
 346 {
 347   if( !this->m_FixedImage ) {
 348     itkExceptionMacro( << "Fixed image has not been assigned" );
 349   }
 350
 351   MeasureType value;
 352   // call the combined version
 353   this->GetValueAndDerivative( parameters, value, derivative );
 354 }
 355
 356 } // end namespace itk
 357
 358
 359 #endif