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 m_FixedImageSamplesIntensityThreshold=0;
38 m_UseFixedImageSamplesIntensityThreshold=false;
39 m_FixedImageMask=NULL;
43 //=========================================================================================================================
45 template <class args_info_type, class FixedImageType, class MovingImageType>
46 typename GenericMetric<args_info_type, FixedImageType, MovingImageType>::MetricPointer
47 GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer()
49 //============================================================================
51 // The image is required for the region and spacing
52 if( ! this->m_FixedImage ) {
53 itkExceptionMacro( "Fixed Image has not been set" );
56 // If the image come from a filter, then update that filter.
57 if( this->m_FixedImage->GetSource() ) {
58 this->m_FixedImage->Update();
62 if( ! m_FixedImageRegionGiven ) {
63 m_FixedImageRegion=m_FixedImage->GetLargestPossibleRegion();
67 //============================================================================
68 //switch on the metric type chosen adn set metric specific members
69 switch (m_ArgsInfo.metric_arg) {
72 typename itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::Pointer m =
73 itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::New();
75 //Set Parameters for this metric
77 if (m_Verbose) std::cout<<"Using the mean squares metric..."<<std::endl;
83 typename clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
84 =clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::New();
86 //Set Parameters for this metric
87 m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
90 if ( !m_ArgsInfo.subtractMean_flag) std::cout<<"Using the normalized correlation metric without subtracting the mean..."<<std::endl;
91 else std::cout<<"Using the normalized correlation metric with subtraction of mean..."<<std::endl;
98 typename itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
99 = itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
101 //Set parameters for this metric
102 typename itk::CorrelationCoefficientHistogramImageToImageMetric<FixedImageType, MovingImageType>::HistogramSizeType size;
103 for (unsigned int i=0; i < FixedImageDimension; i++)size[i]=m_ArgsInfo.bins_arg;
104 m->SetHistogramSize(size);
106 if (m_Verbose) std::cout<<"Using the histogram correlation coefficient metric..."<<std::endl;
112 typename itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
113 = itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::New();
115 //Set parameters for this metric
117 if (m_Verbose) std::cout<<"Using the gradient difference metric..."<<std::endl;
123 typename itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
124 = itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
126 //Set parameters for this metric
127 m->SetFixedImageStandardDeviation( m_ArgsInfo.stdDev_arg );
128 m->SetMovingImageStandardDeviation( m_ArgsInfo.stdDev_arg);
131 //Randomize samples if demanded
132 if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
133 else m->ReinitializeSeed(0);
134 if (m_Verbose) std::cout<<"Using Viola-Wells MI..."<<std::endl;
140 typename itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
141 = itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType>::New();
143 //Set parameters for this metric
144 typename itk::MutualInformationHistogramImageToImageMetric<FixedImageType, MovingImageType>::HistogramSizeType size;
145 for (unsigned int i=0; i < FixedImageDimension; i++)size[i]=m_ArgsInfo.bins_arg;
146 m->SetHistogramSize(size);
148 if (m_Verbose) std::cout<<"Using the histogram MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
154 typename itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
155 = itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
157 //Set parameters for this metric
159 m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
161 //Randomize samples if demanded
162 if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
163 else m->ReinitializeSeed(0);
165 // Two ways of calculating the derivatives
166 m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
170 std::cout<<"Using Mattes' MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
171 if (m_ArgsInfo.explicitPDFDerivatives_flag) std::cout<<"Calculating PDFs explicitly..."<<std::endl;
178 typename itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
179 =itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
181 //Set parameters for this metric
182 typename itk::NormalizedMutualInformationHistogramImageToImageMetric<FixedImageType, MovingImageType>::HistogramSizeType size;
183 for (unsigned int i=0; i < FixedImageDimension; i++)size[i]=m_ArgsInfo.bins_arg;
184 m->SetHistogramSize(size);
186 if (m_Verbose) std::cout<<"Using the normalized MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
192 typename clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
193 = clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::New();
195 //Set parameters for this metric
197 m->SetNumberOfBins(m_ArgsInfo.bins_arg);
198 if (m_Verbose) std::cout<<"Using the correlation ratio..."<<std::endl;
204 typename itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
205 = itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
207 //Set Parameters for this metric
209 if (m_Verbose) std::cout<<"Using the mean squares metric for 3D BLUT FFD..."<<std::endl;
215 typename clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
216 =clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
218 //Set Parameters for this metric
219 m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
222 if ( !m_ArgsInfo.subtractMean_flag) std::cout<<"Using the normalized correlation metric for 3D BLUT FFD without subtracting the mean..."<<std::endl;
223 else std::cout<<"Using the normalized correlation metric 3D BLUT FFD with subtraction of mean..."<<std::endl;
230 typename itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
231 = itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
233 //Set parameters for this metric
235 m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
237 //Randomize samples if demanded
238 if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
239 else m->ReinitializeSeed(0);
241 // Two ways of calculating the derivatives
242 m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
246 std::cout<<"Using Mattes' MI for 3D BLUT FFD with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
247 if (m_ArgsInfo.explicitPDFDerivatives_flag) std::cout<<"Calculating PDFs explicitly..."<<std::endl;
257 if( m_FixedImageMask.IsNotNull() ) m_Metric->SetFixedImageMask(m_FixedImageMask);
258 m_Metric->SetFixedImageRegion(m_FixedImageRegion);
261 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
263 //============================================================================
264 // Set the lower intensity threshold
265 if (m_ArgsInfo.intThreshold_given) {
266 m_UseFixedImageSamplesIntensityThreshold=true;
267 m_FixedImageSamplesIntensityThreshold=m_ArgsInfo.intThreshold_arg;
268 m_Metric->SetFixedImageSamplesIntensityThreshold(m_FixedImageSamplesIntensityThreshold);
269 if (m_Verbose) std::cout<<"Setting the fixed image intensity threshold to "<<m_FixedImageSamplesIntensityThreshold<<"..."<<std::endl;
273 //============================================================================
274 // Set the number of samples
277 if ( ( m_ArgsInfo.samples_arg==1.0) && (m_FixedImageMask.IsNull()) && (!m_UseFixedImageSamplesIntensityThreshold ) ) {
278 m_Metric->SetUseAllPixels(true);
279 if (m_Verbose) std::cout<<"Using all pixels (a fraction of "<<m_ArgsInfo.samples_arg<<")..."<<std::endl;
281 // JV the optimized metric will resample points to obtain the number of pixels:
282 // Pass the correct number of spatial samples, with indexes
284 std::vector<typename FixedImageType::IndexType> fiic;// fixedImageindexContainer
285 unsigned int numberOfValidPixels=0;
286 FixedImageIndexType index;
287 FixedImagePointType inputPoint;
289 // Calculate the number
290 const unsigned int totalNumberOfPixels = m_FixedImageRegion.GetNumberOfPixels();
291 const unsigned int numberOfDemandedPixels = static_cast< unsigned int >( (double) totalNumberOfPixels *m_ArgsInfo.samples_arg );
293 // --------------------------------------------------
294 // Sample whole image sequentially and pass the indexes
295 if (m_ArgsInfo.samples_arg==1.0) {
296 if (m_Verbose) std::cout<<"Sequentially scanning the image for valid pixels..."<<std::endl;
298 // Set up a region interator within the user specified fixed image region.
299 typedef ImageRegionConstIteratorWithIndex<FixedImageType> RegionIterator;
300 RegionIterator regionIter( m_FixedImage, m_FixedImageRegion );
302 // go over the whole region
303 regionIter.GoToBegin();
304 while(!regionIter.IsAtEnd() ) {
307 index = regionIter.GetIndex();
310 if( m_FixedImageMask.IsNotNull() ) {
311 // Check if the Index is inside the mask, translate index to point
312 m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
314 // If not inside the mask, ignore the point
315 if( !m_FixedImageMask->IsInside( inputPoint ) ) {
316 ++regionIter; // jump to next pixel
322 if( m_UseFixedImageSamplesIntensityThreshold &&
323 ( regionIter.Get() < m_FixedImageSamplesIntensityThreshold) ) {
324 ++regionIter; // jump to next pixel
328 // Add point to the numbers
329 fiic.push_back(index);
330 ++numberOfValidPixels;
335 // --------------------------------------------------
338 if (m_Verbose) std::cout<<"Randomly scanning the image for valid pixels..."<<std::endl;
340 // Set up a random interator within the user specified fixed image region.
341 typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
342 RandomIterator randIter( m_FixedImage, m_FixedImageRegion );
344 // Randomly sample the image
345 randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
346 randIter.GoToBegin();
347 while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels) ) {
349 index = randIter.GetIndex();
352 if( m_FixedImageMask.IsNotNull() ) {
354 // Check if the Index is inside the mask, translate index to point
355 m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
357 // If not inside the mask, ignore the point
358 if( !m_FixedImageMask->IsInside( inputPoint ) ) {
359 ++randIter; // jump to next pixel
366 if( m_UseFixedImageSamplesIntensityThreshold &&
367 randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
372 // Add point to the numbers
373 fiic.push_back(index);
374 ++numberOfValidPixels;
380 // Set the indexes of valid samples
381 m_Metric->SetFixedImageIndexes(fiic);
382 // m_Metric->SetNumberOfSpatialSamples( numberOfValidPixels);
383 if (m_Verbose) std::cout<<"A fraction of "<<m_ArgsInfo.samples_arg<<" spatial samples was requested..."<<std::endl;
384 double fraction=(double)numberOfValidPixels/ (double) totalNumberOfPixels;
385 if (m_Verbose) std::cout<<"Found "<<numberOfValidPixels <<" valid pixels for a total of "<<totalNumberOfPixels<<" (a fraction of "<<fraction<<")..."<<std::endl;
390 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;
394 //============================================================================