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 #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 = dynamic_cast<const ImageMaskSpatialObjectType*>(m_FixedImageMask.GetPointer());
 261
 262   typedef typename ImageMaskSpatialObjectType::RegionType ImageMaskRegionType;
 263   ImageMaskRegionType mask_region = mask->GetAxisAlignedBoundingBoxRegion();
 264
 265   // Common properties
 266   if( m_FixedImageMask.IsNotNull() )
 267     m_Metric->SetFixedImageMask(m_FixedImageMask);
 268
 269   m_Metric->SetFixedImageRegion(m_FixedImageRegion);
 270   //m_Metric->SetFixedImageRegion(mask_region);
 271
 272
 273 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
 274
 275   //============================================================================
 276   // Set the lower intensity threshold
 277   if (m_ArgsInfo.intThreshold_given) {
 278     m_UseFixedImageSamplesIntensityThreshold=true;
 279     m_FixedImageSamplesIntensityThreshold=m_ArgsInfo.intThreshold_arg;
 280     m_Metric->SetFixedImageSamplesIntensityThreshold(m_FixedImageSamplesIntensityThreshold);
 281     if (m_Verbose) std::cout<<"Setting the fixed image intensity threshold to "<<m_FixedImageSamplesIntensityThreshold<<"..."<<std::endl;
 282   }
 283
 284
 285   //============================================================================
 286   // Set the number of samples
 287
 288   // Sample all pixel
 289   if ( ( m_ArgsInfo.samples_arg==1.0) && (m_FixedImageMask.IsNull()) && (!m_UseFixedImageSamplesIntensityThreshold ) ) {
 290     m_Metric->SetUseAllPixels(true);
 291     if (m_Verbose) std::cout<<"Using all pixels (a fraction of "<<m_ArgsInfo.samples_arg<<")..."<<std::endl;
 292   }
 293   // JV the optimized metric will resample points to obtain the number of pixels:
 294   // Pass the correct number of spatial samples, with indexes
 295   else {
 296     std::vector<typename FixedImageType::IndexType> fiic;// fixedImageindexContainer
 297     unsigned int numberOfValidPixels=0;
 298     FixedImageIndexType index;
 299     FixedImagePointType inputPoint;
 300
 301     // Calculate the number
 302     const unsigned int totalNumberOfPixels = m_FixedImageRegion.GetNumberOfPixels();
 303     //const unsigned int totalNumberOfPixels = mask_region.GetNumberOfPixels();
 304     const unsigned int numberOfDemandedPixels =  static_cast< unsigned int >( (double) totalNumberOfPixels *m_ArgsInfo.samples_arg );
 305
 306     // --------------------------------------------------
 307     // Sample whole image sequentially and pass the indexes
 308     if (m_ArgsInfo.samples_arg==1.0) {
 309       if (m_Verbose) std::cout<<"Sequentially scanning the image for valid pixels..."<<std::endl;
 310
 311       // Set up a region interator within the user specified fixed image region.
 312       typedef ImageRegionConstIteratorWithIndex<FixedImageType> RegionIterator;
 313       RegionIterator regionIter( m_FixedImage, m_FixedImageRegion );
 314       //RegionIterator regionIter( m_FixedImage, mask_region );
 315
 316       // go over the whole region
 317       regionIter.GoToBegin();
 318       while(!regionIter.IsAtEnd() ) {
 319
 320         // Get sampled index
 321         index = regionIter.GetIndex();
 322
 323         // Mask?
 324         if( m_FixedImageMask.IsNotNull() ) {
 325           // Check if the Index is inside the mask, translate index to point
 326           m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 327
 328           // If not inside the mask, ignore the point
 329           if( !m_FixedImageMask->IsInside( inputPoint ) ) {
 330             ++regionIter; // jump to next pixel
 331             continue;
 332           }
 333         }
 334
 335         // Intensity?
 336         if( m_UseFixedImageSamplesIntensityThreshold &&
 337             ( regionIter.Get() < m_FixedImageSamplesIntensityThreshold) ) {
 338           ++regionIter; // jump to next pixel
 339           continue;
 340         }
 341
 342         // Add point to the numbers
 343         fiic.push_back(index);
 344         ++numberOfValidPixels;
 345         ++regionIter;
 346       }
 347     }
 348
 349     // --------------------------------------------------
 350     // Sample randomly
 351     else {
 352       if (m_Verbose) std::cout<<"Randomly scanning the image for valid pixels..."<<std::endl;
 353
 354       // Set up a random interator within the user specified fixed image region.
 355       typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
 356
 357
 358       RandomIterator randIter( m_FixedImage, m_FixedImageRegion );
 359       //RandomIterator randIter( m_FixedImage, mask_region );
 360
 361       if (m_Verbose) std::cout << "Search region " << m_FixedImageRegion << std::endl;
 362       if (m_Verbose) std::cout << "Mask search region " << mask_region << std::endl;
 363
 364       // Randomly sample the image
 365       short att = 1;
 366       short natts = 5;
 367       while (att <= natts) {
 368         if (m_Verbose) std::cout << "Attempt " << att << std::endl;
 369
 370         int count_out = 0;
 371         int count_not_thres = 0;
 372         clock_t c = clock();
 373         randIter.ReinitializeSeed(c);
 374         randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
 375         randIter.GoToBegin();
 376         //int cnt = 0;
 377         while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels)  ) {
 378           // Get sampled index
 379           index = randIter.GetIndex();
 380           //if (m_Verbose) std::cout << "testing pixel " << index << std::endl;
 381
 382           // Mask?
 383           if( m_FixedImageMask.IsNotNull() ) {
 384
 385             // Check if the Index is inside the mask, translate index to point
 386             m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 387
 388             // If not inside the mask, ignore the point
 389             if( !m_FixedImageMask->IsInside( inputPoint ) ) {
 390               ++randIter; // jump to next pixel
 391               //if (m_Verbose) std::cout << "not inside " << inputPoint << std::endl;
 392               count_out++;
 393               continue;
 394             }
 395
 396           }
 397
 398           // Intensity?
 399           if( m_UseFixedImageSamplesIntensityThreshold &&
 400               randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
 401             ++randIter;
 402               //if (m_Verbose) std::cout << "not in threshold" << std::endl;
 403               count_not_thres++;
 404             continue;
 405           }
 406
 407           // Add point to the numbers
 408           fiic.push_back(index);
 409           ++numberOfValidPixels;
 410           ++randIter;
 411         }
 412
 413         if (m_Verbose) std::cout << "points outside = " << count_out << std::endl;
 414         if (m_Verbose) std::cout << "points not in threshold = " << count_not_thres << std::endl;
 415
 416         if (fiic.size())
 417            break;
 418
 419         att++;
 420       }
 421     }
 422
 423
 424     // Set the indexes of valid samples
 425     m_Metric->SetFixedImageIndexes(fiic);
 426     // m_Metric->SetNumberOfSpatialSamples( numberOfValidPixels);
 427     if (m_Verbose) std::cout<<"A fraction of "<<m_ArgsInfo.samples_arg<<" spatial samples was requested..."<<std::endl;
 428     double fraction=(double)numberOfValidPixels/ (double) totalNumberOfPixels;
 429     if (m_Verbose) std::cout<<"Found "<<numberOfValidPixels <<" valid pixels for a total of "<<totalNumberOfPixels<<" (a fraction of "<<fraction<<")..."<<std::endl;
 430
 431   }
 432
 433 #else
 434   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;
 435
 436
 437 #endif
 438   //============================================================================
 439   //return the pointer
 440   return m_Metric;
 441 }
 442
 443 }
 444
 445 #endif