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 /*=========================================================================
21 Program: Insight Segmentation & Registration Toolkit
22 Module: $RCSfile: itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h,v $
24 Date: $Date: 2010/06/14 17:32:07 $
25 Version: $Revision: 1.1 $
27 Copyright (c) Insight Software Consortium. All rights reserved.
28 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
30 This software is distributed WITHOUT ANY WARRANTY; without even
31 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
32 PURPOSE. See the above copyright notices for more information.
34 =========================================================================*/
35 #ifndef __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_h
36 #define __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_h
38 // First make sure that the configuration is available.
39 // This line can be removed once the optimized versions
40 // gets integrated into the main directories.
41 #include "itkConfigure.h"
43 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
44 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h"
47 #include "itkImageToImageMetric.h"
48 #include "itkCovariantVector.h"
51 #include "itkBSplineKernelFunction.h"
52 #include "itkBSplineDerivativeKernelFunction.h"
53 #include "itkCentralDifferenceImageFunction.h"
54 #include "itkBSplineInterpolateImageFunction.h"
55 #if ITK_VERSION_MAJOR >= 4
56 #include "itkBSplineTransform.h"
58 #include "itkBSplineDeformableTransform.h"
60 #include "itkArray2D.h"
65 /** \class MattesMutualInformationImageToImageMetricFor3DBLUTFFD
66 * \brief Computes the mutual information between two images to be
67 * registered using the method of Mattes et al.
69 * MattesMutualInformationImageToImageMetricFor3DBLUTFFD computes the mutual
70 * information between a fixed and moving image to be registered.
72 * This class is templated over the FixedImage type and the MovingImage
75 * The fixed and moving images are set via methods SetFixedImage() and
76 * SetMovingImage(). This metric makes use of user specified Transform and
77 * Interpolator. The Transform is used to map points from the fixed image to
78 * the moving image domain. The Interpolator is used to evaluate the image
79 * intensity at user specified geometric points in the moving image.
80 * The Transform and Interpolator are set via methods SetTransform() and
83 * If a BSplineInterpolationFunction is used, this class obtain
84 * image derivatives from the BSpline interpolator. Otherwise,
85 * image derivatives are computed using central differencing.
87 * \warning This metric assumes that the moving image has already been
88 * connected to the interpolator outside of this class.
90 * The method GetValue() computes of the mutual information
91 * while method GetValueAndDerivative() computes
92 * both the mutual information and its derivatives with respect to the
93 * transform parameters.
95 * The calculations are based on the method of Mattes et al [1,2]
96 * where the probability density distribution are estimated using
97 * Parzen histograms. Since the fixed image PDF does not contribute
98 * to the derivatives, it does not need to be smooth. Hence,
99 * a zero order (box car) BSpline kernel is used
100 * for the fixed image intensity PDF. On the other hand, to ensure
101 * smoothness a third order BSpline kernel is used for the
102 * moving image intensity PDF.
104 * On Initialize(), the FixedImage is uniformly sampled within
105 * the FixedImageRegion. The number of samples used can be set
106 * via SetNumberOfSpatialSamples(). Typically, the number of
107 * spatial samples used should increase with the image size.
109 * The option UseAllPixelOn() disables the random sampling and uses
110 * all the pixels of the FixedImageRegion in order to estimate the
111 * joint intensity PDF.
113 * During each call of GetValue(), GetDerivatives(),
114 * GetValueAndDerivatives(), marginal and joint intensity PDF's
115 * values are estimated at discrete position or bins.
116 * The number of bins used can be set via SetNumberOfHistogramBins().
117 * To handle data with arbitray magnitude and dynamic range,
118 * the image intensity is scale such that any contribution to the
119 * histogram will fall into a valid bin.
121 * One the PDF's have been contructed, the mutual information
122 * is obtained by doubling summing over the discrete PDF values.
126 * 1. This class returns the negative mutual information value.
127 * 2. This class in not thread safe due the private data structures
128 * used to the store the sampled points and the marginal and joint pdfs.
131 * [1] "Nonrigid multimodality image registration"
132 * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
133 * Medical Imaging 2001: Image Processing, 2001, pp. 1609-1620.
134 * [2] "PET-CT Image Registration in the Chest Using Free-form Deformations"
135 * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
136 * IEEE Transactions in Medical Imaging. Vol.22, No.1,
137 January 2003. pp.120-128.
138 * [3] "Optimization of Mutual Information for MultiResolution Image
140 * P. Thevenaz and M. Unser
141 * IEEE Transactions in Image Processing, 9(12) December 2000.
143 * \ingroup RegistrationMetrics
144 * \ingroup ThreadUnSafe
146 template <class TFixedImage,class TMovingImage >
147 class ITK_EXPORT MattesMutualInformationImageToImageMetricFor3DBLUTFFD :
148 public ImageToImageMetric< TFixedImage, TMovingImage >
152 /** Standard class typedefs. */
153 typedef MattesMutualInformationImageToImageMetricFor3DBLUTFFD Self;
154 typedef ImageToImageMetric< TFixedImage, TMovingImage > Superclass;
155 typedef SmartPointer<Self> Pointer;
156 typedef SmartPointer<const Self> ConstPointer;
158 /** Method for creation through the object factory. */
161 /** Run-time type information (and related methods). */
162 itkTypeMacro(MattesMutualInformationImageToImageMetricFor3DBLUTFFD, ImageToImageMetric);
164 /** Types inherited from Superclass. */
165 typedef typename Superclass::TransformType TransformType;
166 typedef typename Superclass::TransformPointer TransformPointer;
167 typedef typename Superclass::TransformJacobianType TransformJacobianType;
168 typedef typename Superclass::InterpolatorType InterpolatorType;
169 typedef typename Superclass::MeasureType MeasureType;
170 typedef typename Superclass::DerivativeType DerivativeType;
171 typedef typename Superclass::ParametersType ParametersType;
172 typedef typename Superclass::FixedImageType FixedImageType;
173 typedef typename Superclass::MovingImageType MovingImageType;
174 typedef typename Superclass::FixedImageConstPointer FixedImageConstPointer;
175 typedef typename Superclass::MovingImageConstPointer MovingImageCosntPointer;
176 typedef typename Superclass::InputPointType InputPointType;
177 typedef typename Superclass::OutputPointType OutputPointType;
179 typedef typename Superclass::CoordinateRepresentationType
180 CoordinateRepresentationType;
182 /** Index and Point typedef support. */
183 typedef typename FixedImageType::IndexType FixedImageIndexType;
184 typedef typename FixedImageIndexType::IndexValueType FixedImageIndexValueType;
185 typedef typename MovingImageType::IndexType MovingImageIndexType;
186 typedef typename TransformType::InputPointType FixedImagePointType;
187 typedef typename TransformType::OutputPointType MovingImagePointType;
189 /** The moving image dimension. */
190 itkStaticConstMacro( MovingImageDimension, unsigned int,
191 MovingImageType::ImageDimension );
194 * Initialize the Metric by
195 * (1) making sure that all the components are present and plugged
196 * together correctly,
197 * (2) uniformly select NumberOfSpatialSamples within
198 * the FixedImageRegion, and
199 * (3) allocate memory for pdf data structures. */
200 virtual void Initialize(void) throw ( ExceptionObject );
202 /** Get the derivatives of the match measure. */
203 void GetDerivative( const ParametersType& parameters,
204 DerivativeType & Derivative ) const;
206 /** Get the value. */
207 MeasureType GetValue( const ParametersType& parameters ) const;
209 /** Get the value and derivatives for single valued optimizers. */
210 void GetValueAndDerivative( const ParametersType& parameters,
212 DerivativeType& Derivative ) const;
214 /** Number of spatial samples to used to compute metric */
215 itkSetClampMacro( NumberOfSpatialSamples, unsigned long,
216 1, NumericTraits<unsigned long>::max() );
217 itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long);
219 /** Number of bins to used in the histogram. Typical value is 50. */
220 itkSetClampMacro( NumberOfHistogramBins, unsigned long,
221 1, NumericTraits<unsigned long>::max() );
222 itkGetConstReferenceMacro( NumberOfHistogramBins, unsigned long);
224 /** Reinitialize the seed of the random number generator that selects the
225 * sample of pixels used for estimating the image histograms and the joint
226 * histogram. By nature, this metric is not deterministic, since at each run
227 * it may select a different set of pixels. By initializing the random number
228 * generator seed to the same value you can restore determinism. On the other
229 * hand, calling the method ReinitializeSeed() without arguments will use the
230 * clock from your machine in order to have a very random initialization of
231 * the seed. This will indeed increase the non-deterministic behavior of the
233 void ReinitializeSeed();
234 void ReinitializeSeed(int);
236 /** Select whether the metric will be computed using all the pixels on the
237 * fixed image region, or only using a set of randomly selected pixels. */
238 itkSetMacro(UseAllPixels,bool);
239 itkGetConstReferenceMacro(UseAllPixels,bool);
240 itkBooleanMacro(UseAllPixels);
242 /** This variable selects the method to be used for computing the Metric
243 * derivatives with respect to the Transform parameters. Two modes of
244 * computation are available. The choice between one and the other is a
245 * trade-off between computation speed and memory allocations. The two modes
246 * are described in detail below:
248 * UseExplicitPDFDerivatives = True
249 * will compute the Metric derivative by first calculating the derivatives of
250 * each one of the Joint PDF bins with respect to each one of the Transform
251 * parameters and then accumulating these contributions in the final metric
252 * derivative array by using a bin-specific weight. The memory required for
253 * storing the intermediate derivatives is a 3D array of doubles with size
254 * equals to the product of (number of histogram bins)^2 times number of
255 * transform parameters. This method is well suited for Transform with a small
256 * number of parameters.
258 * UseExplicitPDFDerivatives = False will compute the Metric derivative by
259 * first computing the weights for each one of the Joint PDF bins and caching
260 * them into an array. Then it will revisit each one of the PDF bins for
261 * computing its weighted contribution to the full derivative array. In this
262 * method an extra 2D array is used for storing the weights of each one of
263 * the PDF bins. This is an array of doubles with size equals to (number of
264 * histogram bins)^2. This method is well suited for Transforms with a large
265 * number of parameters, such as, BSplineDeformableTransforms. */
266 itkSetMacro(UseExplicitPDFDerivatives,bool);
267 itkGetConstReferenceMacro(UseExplicitPDFDerivatives,bool);
268 itkBooleanMacro(UseExplicitPDFDerivatives);
270 /** This boolean flag is only relevant when this metric is used along
271 * with a BSplineDeformableTransform. The flag enables/disables the
272 * caching of values computed when a physical point is mapped through
273 * the BSplineDeformableTransform. In particular it will cache the
274 * values of the BSpline weights for that points, and the indexes
275 * indicating what BSpline-grid nodes are relevant for that specific
276 * point. This caching is made optional due to the fact that the
277 * memory arrays used for the caching can reach large sizes even for
278 * moderate image size problems. For example, for a 3D image of
279 * 256^3, using 20% of pixels, these arrays will take about 1
280 * Gigabyte of RAM for storage. The ratio of computing time between
281 * using the cache or not using the cache can reach 1:5, meaning that
282 * using the caching can provide a five times speed up. It is
283 * therefore, interesting to enable the caching, if enough memory is
284 * available for it. The caching is enabled by default, in order to
285 * preserve backward compatibility with previous versions of ITK. */
286 itkSetMacro(UseCachingOfBSplineWeights,bool);
287 itkGetConstReferenceMacro(UseCachingOfBSplineWeights,bool);
288 itkBooleanMacro(UseCachingOfBSplineWeights);
292 MattesMutualInformationImageToImageMetricFor3DBLUTFFD();
293 virtual ~MattesMutualInformationImageToImageMetricFor3DBLUTFFD() {};
294 void PrintSelf(std::ostream& os, Indent indent) const;
296 /** \class FixedImageSpatialSample
297 * A fixed image spatial sample consists of the fixed domain point
298 * and the fixed image value at that point. */
300 class FixedImageSpatialSample
303 FixedImageSpatialSample():FixedImageValue(0.0) {
304 FixedImagePointValue.Fill(0.0);
306 ~FixedImageSpatialSample() {};
308 FixedImagePointType FixedImagePointValue;
309 double FixedImageValue;
310 unsigned int FixedImageParzenWindowIndex;
314 /** FixedImageSpatialSample typedef support. */
315 typedef std::vector<FixedImageSpatialSample>
316 FixedImageSpatialSampleContainer;
318 /** Container to store a set of points and fixed image values. */
319 FixedImageSpatialSampleContainer m_FixedImageSamples;
321 /** Uniformly select a sample set from the fixed image domain. */
322 virtual void SampleFixedImageDomain(
323 FixedImageSpatialSampleContainer& samples);
325 /** Gather all the pixels from the fixed image domain. */
326 virtual void SampleFullFixedImageDomain(
327 FixedImageSpatialSampleContainer& samples);
329 /** Transform a point from FixedImage domain to MovingImage domain.
330 * This function also checks if mapped point is within support region. */
331 virtual void TransformPoint( unsigned int sampleNumber,
332 const ParametersType& parameters,
333 MovingImagePointType& mappedPoint,
334 bool& sampleWithinSupportRegion,
335 double& movingImageValue ) const;
339 //purposely not implemented
340 MattesMutualInformationImageToImageMetricFor3DBLUTFFD(const Self&);
341 //purposely not implemented
342 void operator=(const Self&);
345 /** The marginal PDFs are stored as std::vector. */
346 typedef float PDFValueType;
347 typedef std::vector<PDFValueType> MarginalPDFType;
349 /** The fixed image marginal PDF. */
350 mutable MarginalPDFType m_FixedImageMarginalPDF;
352 /** The moving image marginal PDF. */
353 mutable MarginalPDFType m_MovingImageMarginalPDF;
355 /** Helper array for storing the values of the JointPDF ratios. */
356 typedef double PRatioType;
357 typedef Array2D< PRatioType > PRatioArrayType;
358 mutable PRatioArrayType m_PRatioArray;
360 /** Helper variable for accumulating the derivative of the metric. */
361 mutable DerivativeType m_MetricDerivative;
363 /** Typedef for the joint PDF and PDF derivatives are stored as ITK Images. */
364 typedef Image<PDFValueType,2> JointPDFType;
365 typedef JointPDFType::IndexType JointPDFIndexType;
366 typedef JointPDFType::PixelType JointPDFValueType;
367 typedef JointPDFType::RegionType JointPDFRegionType;
368 typedef JointPDFType::SizeType JointPDFSizeType;
369 typedef Image<PDFValueType,3> JointPDFDerivativesType;
370 typedef JointPDFDerivativesType::IndexType JointPDFDerivativesIndexType;
371 typedef JointPDFDerivativesType::PixelType JointPDFDerivativesValueType;
372 typedef JointPDFDerivativesType::RegionType JointPDFDerivativesRegionType;
373 typedef JointPDFDerivativesType::SizeType JointPDFDerivativesSizeType;
375 /** The joint PDF and PDF derivatives. */
376 typename JointPDFType::Pointer m_JointPDF;
377 typename JointPDFDerivativesType::Pointer m_JointPDFDerivatives;
379 unsigned long m_NumberOfSpatialSamples;
380 unsigned long m_NumberOfParameters;
382 /** Variables to define the marginal and joint histograms. */
383 unsigned long m_NumberOfHistogramBins;
384 double m_MovingImageNormalizedMin;
385 double m_FixedImageNormalizedMin;
386 double m_MovingImageTrueMin;
387 double m_MovingImageTrueMax;
388 double m_FixedImageBinSize;
389 double m_MovingImageBinSize;
391 /** Typedefs for BSpline kernel and derivative functions. */
392 typedef BSplineKernelFunction<3> CubicBSplineFunctionType;
393 typedef BSplineDerivativeKernelFunction<3> CubicBSplineDerivativeFunctionType;
395 /** Cubic BSpline kernel for computing Parzen histograms. */
396 typename CubicBSplineFunctionType::Pointer m_CubicBSplineKernel;
397 typename CubicBSplineDerivativeFunctionType::Pointer
398 m_CubicBSplineDerivativeKernel;
400 /** Precompute fixed image parzen window indices. */
401 virtual void ComputeFixedImageParzenWindowIndices(
402 FixedImageSpatialSampleContainer& samples );
405 * Types and variables related to image derivative calculations.
406 * If a BSplineInterpolationFunction is used, this class obtain
407 * image derivatives from the BSpline interpolator. Otherwise,
408 * image derivatives are computed using central differencing.
410 typedef CovariantVector< double,
411 itkGetStaticConstMacro(MovingImageDimension) >
412 ImageDerivativesType;
414 /** Compute image derivatives at a point. */
415 virtual void ComputeImageDerivatives( const MovingImagePointType& mappedPoint,
416 ImageDerivativesType& gradient ) const;
418 /** Boolean to indicate if the interpolator BSpline. */
419 bool m_InterpolatorIsBSpline;
421 /** Typedefs for using BSpline interpolator. */
423 BSplineInterpolateImageFunction<MovingImageType,
424 CoordinateRepresentationType>
425 BSplineInterpolatorType;
427 /** Pointer to BSplineInterpolator. */
428 typename BSplineInterpolatorType::Pointer m_BSplineInterpolator;
430 /** Typedefs for using central difference calculator. */
431 typedef CentralDifferenceImageFunction<MovingImageType,
432 CoordinateRepresentationType>
433 DerivativeFunctionType;
435 /** Pointer to central difference calculator. */
436 typename DerivativeFunctionType::Pointer m_DerivativeCalculator;
439 /** Compute PDF derivative contribution for each parameter. */
440 virtual void ComputePDFDerivatives( unsigned int sampleNumber,
441 int movingImageParzenWindowIndex,
442 const ImageDerivativesType&
443 movingImageGradientValue,
444 double cubicBSplineDerivativeValue
448 * Types and variables related to BSpline deformable transforms.
449 * If the transform is of type third order BSplineDeformableTransform,
450 * then we can speed up the metric derivative calculation by
451 * only inspecting the parameters within the support region
455 /** Boolean to indicate if the transform is BSpline deformable. */
456 bool m_TransformIsBSpline;
458 /** The number of BSpline parameters per image dimension. */
459 long m_NumParametersPerDim;
462 * The number of BSpline transform weights is the number of
463 * of parameter in the support region (per dimension ). */
464 unsigned long m_NumBSplineWeights;
466 /** The fixed image dimension. */
467 itkStaticConstMacro( FixedImageDimension, unsigned int,
468 FixedImageType::ImageDimension );
471 * Enum of the deformabtion field spline order.
473 enum { DeformationSplineOrder = 3 };
476 * Typedefs for the BSplineDeformableTransform.
478 #if ITK_VERSION_MAJOR >= 4
479 typedef BSplineTransform<
480 CoordinateRepresentationType,
481 ::itk::GetImageDimension<FixedImageType>::ImageDimension,
482 DeformationSplineOrder> BSplineTransformType;
484 typedef BSplineDeformableTransform<
485 CoordinateRepresentationType,
486 ::itk::GetImageDimension<FixedImageType>::ImageDimension,
487 DeformationSplineOrder> BSplineTransformType;
489 typedef typename BSplineTransformType::WeightsType
490 BSplineTransformWeightsType;
491 typedef typename BSplineTransformType::ParameterIndexArrayType
492 BSplineTransformIndexArrayType;
495 * Variables used when transform is of type BSpline deformable.
497 typename BSplineTransformType::Pointer m_BSplineTransform;
500 * Cache pre-transformed points, weights, indices and
501 * within support region flag.
503 typedef typename BSplineTransformWeightsType::ValueType WeightsValueType;
504 typedef Array2D<WeightsValueType> BSplineTransformWeightsArrayType;
505 typedef typename BSplineTransformIndexArrayType::ValueType IndexValueType;
506 typedef Array2D<IndexValueType> BSplineTransformIndicesArrayType;
507 typedef std::vector<MovingImagePointType> MovingImagePointArrayType;
508 typedef std::vector<bool> BooleanArrayType;
510 BSplineTransformWeightsArrayType m_BSplineTransformWeightsArray;
511 BSplineTransformIndicesArrayType m_BSplineTransformIndicesArray;
512 MovingImagePointArrayType m_PreTransformPointsArray;
513 BooleanArrayType m_WithinSupportRegionArray;
515 typedef FixedArray<unsigned long,
516 ::itk::GetImageDimension<FixedImageType>::ImageDimension>
517 ParametersOffsetType;
518 ParametersOffsetType m_ParametersOffset;
522 virtual void PreComputeTransformValues();
524 bool m_ReseedIterator;
527 // Selection of explicit or implicit computation of PDF derivatives
528 // with respect to Transform parameters.
529 bool m_UseExplicitPDFDerivatives;
531 // Variables needed for optionally caching values when using a BSpline transform.
532 bool m_UseCachingOfBSplineWeights;
533 mutable BSplineTransformWeightsType m_Weights;
534 mutable BSplineTransformIndexArrayType m_Indices;
538 } // end namespace itk
540 #ifndef ITK_MANUAL_INSTANTIATION
541 #include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx"