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