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