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