1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
19 #ifndef __clitkGenericMetric_txx
20 #define __clitkGenericMetric_txx
22 #include "clitkGenericMetric.h"
28 //=========================================================================================================================
30 template <class args_info_type, class FixedImageType, class MovingImageType>
31 GenericMetric<args_info_type, FixedImageType, MovingImageType>::GenericMetric()
36 m_FixedImageRegionGiven=false;
37 #ifdef ITK_USE_OPTIMISED_REGISTRATION_METHODS
38 m_FixedImageSamplesIntensityThreshold=0;
39 m_UseFixedImageSamplesIntensityThreshold=false;
41 m_FixedImageMask=NULL;
45 //=========================================================================================================================
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()
51 //============================================================================
53 // The image is required for the region and spacing
54 if( ! this->m_FixedImage ) {
55 itkExceptionMacro( "Fixed Image has not been set" );
58 // If the image come from a filter, then update that filter.
59 if( this->m_FixedImage->GetSource() ) {
60 this->m_FixedImage->Update();
64 if( ! m_FixedImageRegionGiven ) {
65 m_FixedImageRegion=m_FixedImage->GetLargestPossibleRegion();
69 //============================================================================
70 //switch on the metric type chosen adn set metric specific members
71 switch (m_ArgsInfo.metric_arg) {
74 typename itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::Pointer m =
75 itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::New();
77 //Set Parameters for this metric
79 if (m_Verbose) std::cout<<"Using the mean squares metric..."<<std::endl;
85 typename clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
86 =clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::New();
88 //Set Parameters for this metric
89 m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
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;
100 typename itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
101 = itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
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);
108 if (m_Verbose) std::cout<<"Using the histogram correlation coefficient metric..."<<std::endl;
114 typename itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
115 = itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::New();
117 //Set parameters for this metric
119 if (m_Verbose) std::cout<<"Using the gradient difference metric..."<<std::endl;
125 typename itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
126 = itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
128 //Set parameters for this metric
129 m->SetFixedImageStandardDeviation( m_ArgsInfo.stdDev_arg );
130 m->SetMovingImageStandardDeviation( m_ArgsInfo.stdDev_arg);
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;
142 typename itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
143 = itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType>::New();
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);
150 if (m_Verbose) std::cout<<"Using the histogram MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
156 typename itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
157 = itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
159 //Set parameters for this metric
161 m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
163 //Randomize samples if demanded
164 if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
165 else m->ReinitializeSeed(0);
167 // Two ways of calculating the derivatives
168 m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
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;
180 typename itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
181 =itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
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);
188 if (m_Verbose) std::cout<<"Using the normalized MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
194 typename clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
195 = clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::New();
197 //Set parameters for this metric
199 m->SetNumberOfBins(m_ArgsInfo.bins_arg);
200 if (m_Verbose) std::cout<<"Using the correlation ratio..."<<std::endl;
206 typename itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
207 = itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
209 //Set Parameters for this metric
211 if (m_Verbose) std::cout<<"Using the mean squares metric for 3D BLUT FFD..."<<std::endl;
217 typename clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
218 =clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
220 //Set Parameters for this metric
221 m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
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;
232 typename itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
233 = itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
235 //Set parameters for this metric
237 m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
239 //Randomize samples if demanded
240 if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
241 else m->ReinitializeSeed(0);
243 // Two ways of calculating the derivatives
244 m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
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;
259 if( m_FixedImageMask.IsNotNull() ) m_Metric->SetFixedImageMask(m_FixedImageMask);
260 m_Metric->SetFixedImageRegion(m_FixedImageRegion);
263 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
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;
275 //============================================================================
276 // Set the number of samples
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;
283 // JV the optimized metric will resample points to obtain the number of pixels:
284 // Pass the correct number of spatial samples, with indexes
286 std::vector<typename FixedImageType::IndexType> fiic;// fixedImageindexContainer
287 unsigned int numberOfValidPixels=0;
288 FixedImageIndexType index;
289 FixedImagePointType inputPoint;
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 );
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;
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 );
304 // go over the whole region
305 regionIter.GoToBegin();
306 while(!regionIter.IsAtEnd() ) {
309 index = regionIter.GetIndex();
312 if( m_FixedImageMask.IsNotNull() ) {
313 // Check if the Index is inside the mask, translate index to point
314 m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
316 // If not inside the mask, ignore the point
317 if( !m_FixedImageMask->IsInside( inputPoint ) ) {
318 ++regionIter; // jump to next pixel
324 if( m_UseFixedImageSamplesIntensityThreshold &&
325 ( regionIter.Get() < m_FixedImageSamplesIntensityThreshold) ) {
326 ++regionIter; // jump to next pixel
330 // Add point to the numbers
331 fiic.push_back(index);
332 ++numberOfValidPixels;
337 // --------------------------------------------------
340 if (m_Verbose) std::cout<<"Randomly scanning the image for valid pixels..."<<std::endl;
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 );
346 // Randomly sample the image
347 randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
348 randIter.GoToBegin();
349 while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels) ) {
351 index = randIter.GetIndex();
354 if( m_FixedImageMask.IsNotNull() ) {
356 // Check if the Index is inside the mask, translate index to point
357 m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
359 // If not inside the mask, ignore the point
360 if( !m_FixedImageMask->IsInside( inputPoint ) ) {
361 ++randIter; // jump to next pixel
368 if( m_UseFixedImageSamplesIntensityThreshold &&
369 randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
374 // Add point to the numbers
375 fiic.push_back(index);
376 ++numberOfValidPixels;
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;
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;
396 //============================================================================