]> Creatis software - clitk.git/blob - registration/clitkGenericMetric.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / 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 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
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 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
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         /*
341         if( m_UseFixedImageSamplesIntensityThreshold &&
342             ( regionIter.Get() < m_FixedImageSamplesIntensityThreshold) ) {
343           ++regionIter; // jump to next pixel
344           continue;
345         }
346         */
347
348         // Add point to the numbers
349         fiic.push_back(index);
350         ++numberOfValidPixels;
351         ++regionIter;
352       }
353     }
354
355     // --------------------------------------------------
356     // Sample randomly
357     else {
358       if (m_Verbose) std::cout<<"Randomly scanning the image for valid pixels..."<<std::endl;
359
360       // Set up a random interator within the user specified fixed image region.
361       typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
362       
363       
364       RandomIterator randIter( m_FixedImage, m_FixedImageRegion );
365       //RandomIterator randIter( m_FixedImage, mask_region );
366       
367       if (m_Verbose) std::cout << "Search region " << m_FixedImageRegion << std::endl;
368       if (m_Verbose) std::cout << "Mask search region " << mask_region << std::endl;
369
370       // Randomly sample the image
371       short att = 1;
372       short natts = 5;
373       while (att <= natts) {
374         if (m_Verbose) std::cout << "Attempt " << att << std::endl;
375           
376         int count_out = 0;
377         int count_not_thres = 0;
378         clock_t c = clock();
379         randIter.ReinitializeSeed(c);
380         randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
381         randIter.GoToBegin();
382         //int cnt = 0;
383         while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels)  ) {
384           // Get sampled index
385           index = randIter.GetIndex();
386           //if (m_Verbose) std::cout << "testing pixel " << index << std::endl;
387
388           // Mask?
389           if( m_FixedImageMask.IsNotNull() ) {
390
391             // Check if the Index is inside the mask, translate index to point
392             m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
393
394             // If not inside the mask, ignore the point
395             if( !m_FixedImageMask->IsInside( inputPoint ) ) {
396               ++randIter; // jump to next pixel
397               //if (m_Verbose) std::cout << "not inside " << inputPoint << std::endl;
398               count_out++;
399               continue;
400             }
401
402           }
403
404           // Intensity?
405 //           if( m_UseFixedImageSamplesIntensityThreshold &&
406 //               randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
407 //             ++randIter;
408 //               //if (m_Verbose) std::cout << "not in threshold" << std::endl;
409 //               count_not_thres++;
410 //             continue;
411 //           }
412
413           // Add point to the numbers
414           fiic.push_back(index);
415           ++numberOfValidPixels;
416           ++randIter;
417         }
418         
419         if (m_Verbose) std::cout << "points outside = " << count_out << std::endl;
420         if (m_Verbose) std::cout << "points not in threshold = " << count_not_thres << std::endl;
421           
422         if (fiic.size())
423            break;
424         
425         att++;
426       }
427     }
428
429
430     // Set the indexes of valid samples
431     m_Metric->SetFixedImageIndexes(fiic);
432     // m_Metric->SetNumberOfSpatialSamples( numberOfValidPixels);
433     if (m_Verbose) std::cout<<"A fraction of "<<m_ArgsInfo.samples_arg<<" spatial samples was requested..."<<std::endl;
434     double fraction=(double)numberOfValidPixels/ (double) totalNumberOfPixels;
435     if (m_Verbose) std::cout<<"Found "<<numberOfValidPixels <<" valid pixels for a total of "<<totalNumberOfPixels<<" (a fraction of "<<fraction<<")..."<<std::endl;
436     if (m_Verbose) std::cout<<"number of mask pixels "<<totalNumberOfMaskPixels<<std::endl;
437
438   }
439
440 #else
441   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;
442
443
444 #endif
445   //============================================================================
446   //return the pointer
447   return m_Metric;
448 }
449
450 }
451
452 #endif