registration/clitkGenericMetric.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://oncora1.lyon.fnclcc.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 #ifndef __clitkGenericMetric_txx
  20 #define __clitkGenericMetric_txx
  21
  22 #include "clitkGenericMetric.h"
  23
  24
  25 namespace clitk
  26 {
  27
  28 //=========================================================================================================================
  29 //constructor
  30 template <class args_info_type, class FixedImageType, class MovingImageType>
  31 GenericMetric<args_info_type, FixedImageType, MovingImageType>::GenericMetric()
  32 {
  33   m_Metric=NULL;
  34   m_Maximize=false;
  35   m_Verbose=false;
  36   m_FixedImageRegionGiven=false;
  37   m_FixedImageSamplesIntensityThreshold=0;
  38   m_UseFixedImageSamplesIntensityThreshold=false;
  39   m_FixedImageMask=NULL;
  40 }
  41
  42
  43 //=========================================================================================================================
  44 //Get the pointer
  45 template <class args_info_type, class FixedImageType, class MovingImageType>
  46 typename GenericMetric<args_info_type, FixedImageType, MovingImageType>::MetricPointer
  47 GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer()
  48 {
  49   //============================================================================
  50   // Sanity check:
  51   // The image is required for the region and spacing
  52   if( ! this->m_FixedImage ) {
  53     itkExceptionMacro( "Fixed Image has not been set" );
  54   }
  55
  56   // If the image come from a filter, then update that filter.
  57   if( this->m_FixedImage->GetSource() ) {
  58     this->m_FixedImage->Update();
  59   }
  60
  61   // The metric region
  62   if( !  m_FixedImageRegionGiven ) {
  63     m_FixedImageRegion=m_FixedImage->GetLargestPossibleRegion();
  64   }
  65
  66
  67   //============================================================================
  68   //switch on  the  metric type chosen adn set metric specific members
  69   switch (m_ArgsInfo.metric_arg) {
  70
  71   case 0: {
  72     typename itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::Pointer m =
  73       itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::New();
  74
  75     //Set Parameters for this metric
  76     m_Maximize=false;
  77     if (m_Verbose) std::cout<<"Using the mean squares metric..."<<std::endl;
  78     m_Metric=m;
  79     break;
  80   }
  81
  82   case 1: {
  83     typename  clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
  84     =clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::New();
  85
  86     //Set Parameters for this metric
  87     m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
  88     m_Maximize=false;
  89     if (m_Verbose) {
  90       if ( !m_ArgsInfo.subtractMean_flag) std::cout<<"Using the normalized correlation metric without subtracting the mean..."<<std::endl;
  91       else  std::cout<<"Using the normalized correlation metric with subtraction of mean..."<<std::endl;
  92     }
  93     m_Metric=m;
  94     break;
  95   }
  96
  97   case 2: {
  98     typename itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
  99     = itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
 100
 101     //Set parameters for this metric
 102     typename itk::CorrelationCoefficientHistogramImageToImageMetric<FixedImageType, MovingImageType>::HistogramSizeType size;
 103     for (unsigned int i=0; i < FixedImageDimension; i++)size[i]=m_ArgsInfo.bins_arg;
 104     m->SetHistogramSize(size);
 105     m_Maximize=false;
 106     if (m_Verbose) std::cout<<"Using the histogram correlation coefficient metric..."<<std::endl;
 107     m_Metric = m;
 108     break;
 109   }
 110
 111   case 3: {
 112     typename itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
 113     = itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::New();
 114
 115     //Set parameters for this metric
 116     m_Maximize=false;
 117     if (m_Verbose) std::cout<<"Using the gradient difference metric..."<<std::endl;
 118     m_Metric=m;
 119     break;
 120   }
 121
 122   case 4: {
 123     typename itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
 124     = itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
 125
 126     //Set parameters for this metric
 127     m->SetFixedImageStandardDeviation( m_ArgsInfo.stdDev_arg );
 128     m->SetMovingImageStandardDeviation( m_ArgsInfo.stdDev_arg);
 129     m_Maximize=true;
 130
 131     //Randomize samples if demanded
 132     if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
 133     else m->ReinitializeSeed(0);
 134     if (m_Verbose) std::cout<<"Using Viola-Wells MI..."<<std::endl;
 135     m_Metric=m;
 136     break;
 137   }
 138
 139   case 5: {
 140     typename itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
 141     = itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType>::New();
 142
 143     //Set parameters for this metric
 144     typename itk::MutualInformationHistogramImageToImageMetric<FixedImageType, MovingImageType>::HistogramSizeType size;
 145     for (unsigned int i=0; i < FixedImageDimension; i++)size[i]=m_ArgsInfo.bins_arg;
 146     m->SetHistogramSize(size);
 147     m_Maximize=true;
 148     if (m_Verbose) std::cout<<"Using the histogram MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
 149     m_Metric=m;
 150     break;
 151   }
 152
 153   case 6: {
 154     typename itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
 155     = itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
 156
 157     //Set parameters for this metric
 158     m_Maximize=false;
 159     m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
 160
 161     //Randomize samples if demanded
 162     if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
 163     else m->ReinitializeSeed(0);
 164
 165     // Two ways of calculating the derivatives
 166     m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
 167
 168
 169     if (m_Verbose) {
 170       std::cout<<"Using Mattes' MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
 171       if (m_ArgsInfo.explicitPDFDerivatives_flag) std::cout<<"Calculating PDFs explicitly..."<<std::endl;
 172     }
 173     m_Metric=m;
 174     break;
 175   }
 176
 177   case 7: {
 178     typename itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
 179     =itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
 180
 181     //Set parameters for this metric
 182     typename itk::NormalizedMutualInformationHistogramImageToImageMetric<FixedImageType, MovingImageType>::HistogramSizeType size;
 183     for (unsigned int i=0; i < FixedImageDimension; i++)size[i]=m_ArgsInfo.bins_arg;
 184     m->SetHistogramSize(size);
 185     m_Maximize=false;
 186     if (m_Verbose) std::cout<<"Using the normalized MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
 187     m_Metric=m;
 188     break;
 189   }
 190
 191   case 8: {
 192     typename clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
 193     = clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::New();
 194
 195     //Set parameters for this metric
 196     m_Maximize=false;
 197     m->SetNumberOfBins(m_ArgsInfo.bins_arg);
 198     if (m_Verbose) std::cout<<"Using the correlation ratio..."<<std::endl;
 199     m_Metric =m;
 200     break;
 201   }
 202
 203   case 9: {
 204     typename itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
 205     = itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
 206
 207     //Set Parameters for this metric
 208     m_Maximize=false;
 209     if (m_Verbose) std::cout<<"Using the mean squares metric for 3D BLUT FFD..."<<std::endl;
 210     m_Metric=m;
 211     break;
 212   }
 213
 214   case 10: {
 215     typename  clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
 216     =clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
 217
 218     //Set Parameters for this metric
 219     m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
 220     m_Maximize=false;
 221     if (m_Verbose) {
 222       if ( !m_ArgsInfo.subtractMean_flag) std::cout<<"Using the normalized correlation metric for 3D BLUT FFD without subtracting the mean..."<<std::endl;
 223       else  std::cout<<"Using the normalized correlation metric 3D BLUT FFD with subtraction of mean..."<<std::endl;
 224     }
 225     m_Metric=m;
 226     break;
 227   }
 228
 229   case 11: {
 230     typename itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
 231     = itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
 232
 233     //Set parameters for this metric
 234     m_Maximize=false;
 235     m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
 236
 237     //Randomize samples if demanded
 238     if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
 239     else m->ReinitializeSeed(0);
 240
 241     // Two ways of calculating the derivatives
 242     m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
 243
 244
 245     if (m_Verbose) {
 246       std::cout<<"Using Mattes' MI for 3D BLUT FFD with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
 247       if (m_ArgsInfo.explicitPDFDerivatives_flag) std::cout<<"Calculating PDFs explicitly..."<<std::endl;
 248     }
 249     m_Metric=m;
 250     break;
 251   }
 252
 253   }
 254
 255
 256   // Common properties
 257   if( m_FixedImageMask.IsNotNull() )  m_Metric->SetFixedImageMask(m_FixedImageMask);
 258   m_Metric->SetFixedImageRegion(m_FixedImageRegion);
 259
 260
 261 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
 262
 263   //============================================================================
 264   // Set the lower intensity threshold
 265   if (m_ArgsInfo.intThreshold_given) {
 266     m_UseFixedImageSamplesIntensityThreshold=true;
 267     m_FixedImageSamplesIntensityThreshold=m_ArgsInfo.intThreshold_arg;
 268     m_Metric->SetFixedImageSamplesIntensityThreshold(m_FixedImageSamplesIntensityThreshold);
 269     if (m_Verbose) std::cout<<"Setting the fixed image intensity threshold to "<<m_FixedImageSamplesIntensityThreshold<<"..."<<std::endl;
 270   }
 271
 272
 273   //============================================================================
 274   // Set the number of samples
 275
 276   // Sample all pixel
 277   if ( ( m_ArgsInfo.samples_arg==1.0) && (m_FixedImageMask.IsNull()) && (!m_UseFixedImageSamplesIntensityThreshold ) ) {
 278     m_Metric->SetUseAllPixels(true);
 279     if (m_Verbose) std::cout<<"Using all pixels (a fraction of "<<m_ArgsInfo.samples_arg<<")..."<<std::endl;
 280   }
 281   // JV the optimized metric will resample points to obtain the number of pixels:
 282   // Pass the correct number of spatial samples, with indexes
 283   else {
 284     std::vector<typename FixedImageType::IndexType> fiic;// fixedImageindexContainer
 285     unsigned int numberOfValidPixels=0;
 286     FixedImageIndexType index;
 287     FixedImagePointType inputPoint;
 288
 289     // Calculate the number
 290     const unsigned int totalNumberOfPixels = m_FixedImageRegion.GetNumberOfPixels();
 291     const unsigned int numberOfDemandedPixels =  static_cast< unsigned int >( (double) totalNumberOfPixels *m_ArgsInfo.samples_arg );
 292
 293     // --------------------------------------------------
 294     // Sample whole image sequentially and pass the indexes
 295     if (m_ArgsInfo.samples_arg==1.0) {
 296       if (m_Verbose) std::cout<<"Sequentially scanning the image for valid pixels..."<<std::endl;
 297
 298       // Set up a region interator within the user specified fixed image region.
 299       typedef ImageRegionConstIteratorWithIndex<FixedImageType> RegionIterator;
 300       RegionIterator regionIter( m_FixedImage, m_FixedImageRegion );
 301
 302       // go over the whole region
 303       regionIter.GoToBegin();
 304       while(!regionIter.IsAtEnd() ) {
 305
 306         // Get sampled index
 307         index = regionIter.GetIndex();
 308
 309         // Mask?
 310         if( m_FixedImageMask.IsNotNull() ) {
 311           // Check if the Index is inside the mask, translate index to point
 312           m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 313
 314           // If not inside the mask, ignore the point
 315           if( !m_FixedImageMask->IsInside( inputPoint ) ) {
 316             ++regionIter; // jump to next pixel
 317             continue;
 318           }
 319         }
 320
 321         // Intensity?
 322         if( m_UseFixedImageSamplesIntensityThreshold &&
 323             ( regionIter.Get() < m_FixedImageSamplesIntensityThreshold) ) {
 324           ++regionIter; // jump to next pixel
 325           continue;
 326         }
 327
 328         // Add point to the numbers
 329         fiic.push_back(index);
 330         ++numberOfValidPixels;
 331         ++regionIter;
 332       }
 333     }
 334
 335     // --------------------------------------------------
 336     // Sample randomly
 337     else {
 338       if (m_Verbose) std::cout<<"Randomly scanning the image for valid pixels..."<<std::endl;
 339
 340       // Set up a random interator within the user specified fixed image region.
 341       typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
 342       RandomIterator randIter( m_FixedImage, m_FixedImageRegion );
 343
 344       // Randomly sample the image
 345       randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
 346       randIter.GoToBegin();
 347       while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels)  ) {
 348         // Get sampled index
 349         index = randIter.GetIndex();
 350
 351         // Mask?
 352         if( m_FixedImageMask.IsNotNull() ) {
 353
 354           // Check if the Index is inside the mask, translate index to point
 355           m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 356
 357           // If not inside the mask, ignore the point
 358           if( !m_FixedImageMask->IsInside( inputPoint ) ) {
 359             ++randIter; // jump to next pixel
 360             continue;
 361           }
 362
 363         }
 364
 365         // Intensity?
 366         if( m_UseFixedImageSamplesIntensityThreshold &&
 367             randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
 368           ++randIter;
 369           continue;
 370         }
 371
 372         // Add point to the numbers
 373         fiic.push_back(index);
 374         ++numberOfValidPixels;
 375         ++randIter;
 376       }
 377     }
 378
 379
 380     // Set the indexes of valid samples
 381     m_Metric->SetFixedImageIndexes(fiic);
 382     // m_Metric->SetNumberOfSpatialSamples( numberOfValidPixels);
 383     if (m_Verbose) std::cout<<"A fraction of "<<m_ArgsInfo.samples_arg<<" spatial samples was requested..."<<std::endl;
 384     double fraction=(double)numberOfValidPixels/ (double) totalNumberOfPixels;
 385     if (m_Verbose) std::cout<<"Found "<<numberOfValidPixels <<" valid pixels for a total of "<<totalNumberOfPixels<<" (a fraction of "<<fraction<<")..."<<std::endl;
 386
 387   }
 388
 389 #else
 390   if (m_Verbose) std::cout<<"Not setting the fixed image intensity threshold or the fraction of samples to use (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
 391
 392
 393 #endif
 394   //============================================================================
 395   //return the pointer
 396   return m_Metric;
 397 }
 398
 399 }
 400
 401 #endif