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://www.centreleonberard.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"
23 #include <itkImageMaskSpatialObject.h>
29 //=========================================================================================================================
31 template <class args_info_type, class FixedImageType, class MovingImageType>
32 GenericMetric<args_info_type, FixedImageType, MovingImageType>::GenericMetric()
37 m_FixedImageRegionGiven=false;
38 m_FixedImageSamplesIntensityThreshold=0;
39 m_UseFixedImageSamplesIntensityThreshold=false;
40 m_FixedImageMask=ITK_NULLPTR;
44 //=========================================================================================================================
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()
50 //============================================================================
52 // The image is required for the region and spacing
53 if( ! this->m_FixedImage ) {
54 itkExceptionMacro( "Fixed Image has not been set" );
57 // If the image come from a filter, then update that filter.
58 if( this->m_FixedImage->GetSource() ) {
59 this->m_FixedImage->Update();
63 if( ! m_FixedImageRegionGiven ) {
64 m_FixedImageRegion=m_FixedImage->GetLargestPossibleRegion();
68 //============================================================================
69 //switch on the metric type chosen adn set metric specific members
70 switch (m_ArgsInfo.metric_arg) {
73 typename itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::Pointer m =
74 itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType >::New();
76 //Set Parameters for this metric
78 if (m_Verbose) std::cout<<"Using the mean squares metric..."<<std::endl;
84 typename clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
85 =clitk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType >::New();
87 //Set Parameters for this metric
88 m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
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;
99 typename itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
100 = itk::CorrelationCoefficientHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
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);
107 if (m_Verbose) std::cout<<"Using the histogram correlation coefficient metric..."<<std::endl;
113 typename itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
114 = itk::GradientDifferenceImageToImageMetric< FixedImageType, MovingImageType >::New();
116 //Set parameters for this metric
118 if (m_Verbose) std::cout<<"Using the gradient difference metric..."<<std::endl;
124 typename itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
125 = itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
127 //Set parameters for this metric
128 m->SetFixedImageStandardDeviation( m_ArgsInfo.stdDev_arg );
129 m->SetMovingImageStandardDeviation( m_ArgsInfo.stdDev_arg);
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;
141 typename itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
142 = itk::MutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType>::New();
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);
149 if (m_Verbose) std::cout<<"Using the histogram MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
155 typename itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
156 = itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType >::New();
158 //Set parameters for this metric
160 m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
162 //Randomize samples if demanded
163 if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
164 else m->ReinitializeSeed(0);
166 // Two ways of calculating the derivatives
167 m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
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;
179 typename itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
180 =itk::NormalizedMutualInformationHistogramImageToImageMetric< FixedImageType, MovingImageType >::New();
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);
187 if (m_Verbose) std::cout<<"Using the normalized MI with "<<m_ArgsInfo.bins_arg<<" bins..."<<std::endl;
193 typename clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::Pointer m
194 = clitk::CorrelationRatioImageToImageMetric< FixedImageType, MovingImageType >::New();
196 //Set parameters for this metric
198 m->SetNumberOfBins(m_ArgsInfo.bins_arg);
199 if (m_Verbose) std::cout<<"Using the correlation ratio..."<<std::endl;
205 typename itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
206 = itk::MeanSquaresImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
208 //Set Parameters for this metric
210 if (m_Verbose) std::cout<<"Using the mean squares metric for 3D BLUT FFD..."<<std::endl;
216 typename clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
217 =clitk::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
219 //Set Parameters for this metric
220 m->SetSubtractMean(m_ArgsInfo.subtractMean_flag);
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;
231 typename itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::Pointer m
232 = itk::MattesMutualInformationImageToImageMetricFor3DBLUTFFD< FixedImageType, MovingImageType >::New();
234 //Set parameters for this metric
236 m->SetNumberOfHistogramBins(m_ArgsInfo.bins_arg);
238 //Randomize samples if demanded
239 if (m_ArgsInfo.random_flag ) m->ReinitializeSeed();
240 else m->ReinitializeSeed(0);
242 // Two ways of calculating the derivatives
243 m->SetUseExplicitPDFDerivatives(m_ArgsInfo.explicitPDFDerivatives_flag);
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;
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());
262 typedef typename ImageMaskSpatialObjectType::RegionType ImageMaskRegionType;
263 ImageMaskRegionType mask_region;
264 if (mask.IsNotNull())
265 mask_region = mask->GetAxisAlignedBoundingBoxRegion();
268 if( m_FixedImageMask.IsNotNull() )
269 m_Metric->SetFixedImageMask(m_FixedImageMask);
271 m_Metric->SetFixedImageRegion(m_FixedImageRegion);
272 //m_Metric->SetFixedImageRegion(mask_region);
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;
284 //============================================================================
285 // Set the number of samples
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;
292 // JV the optimized metric will resample points to obtain the number of pixels:
293 // Pass the correct number of spatial samples, with indexes
295 std::vector<typename FixedImageType::IndexType> fiic;// fixedImageindexContainer
296 unsigned int numberOfValidPixels=0;
297 FixedImageIndexType index;
298 FixedImagePointType inputPoint;
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 );
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;
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 );
315 // go over the whole region
316 regionIter.GoToBegin();
317 while(!regionIter.IsAtEnd() ) {
320 index = regionIter.GetIndex();
323 if( m_FixedImageMask.IsNotNull() ) {
324 // Check if the Index is inside the mask, translate index to point
325 m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
327 // If not inside the mask, ignore the point
328 if( !m_FixedImageMask->IsInside( inputPoint ) ) {
329 ++regionIter; // jump to next pixel
336 if( m_UseFixedImageSamplesIntensityThreshold &&
337 ( regionIter.Get() < m_FixedImageSamplesIntensityThreshold) ) {
338 ++regionIter; // jump to next pixel
343 // Add point to the numbers
344 fiic.push_back(index);
345 ++numberOfValidPixels;
350 // --------------------------------------------------
353 if (m_Verbose) std::cout<<"Randomly scanning the image for valid pixels..."<<std::endl;
355 // Set up a random interator within the user specified fixed image region.
356 typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
359 RandomIterator randIter( m_FixedImage, m_FixedImageRegion );
360 //RandomIterator randIter( m_FixedImage, mask_region );
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;
365 // Randomly sample the image
368 while (att <= natts) {
369 if (m_Verbose) std::cout << "Attempt " << att << std::endl;
372 int count_not_thres = 0;
374 randIter.ReinitializeSeed(c);
375 randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
376 randIter.GoToBegin();
378 while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels) ) {
380 index = randIter.GetIndex();
381 //if (m_Verbose) std::cout << "testing pixel " << index << std::endl;
384 if( m_FixedImageMask.IsNotNull() ) {
386 // Check if the Index is inside the mask, translate index to point
387 m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
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;
400 // if( m_UseFixedImageSamplesIntensityThreshold &&
401 // randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
403 // //if (m_Verbose) std::cout << "not in threshold" << std::endl;
404 // count_not_thres++;
408 // Add point to the numbers
409 fiic.push_back(index);
410 ++numberOfValidPixels;
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;
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;
434 //============================================================================