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 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
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 #include "itkBSplineDeformableTransform.h"
56 #include "itkArray2D.h"
61 /** \class MattesMutualInformationImageToImageMetricFor3DBLUTFFD
62 * \brief Computes the mutual information between two images to be
63 * registered using the method of Mattes et al.
65 * MattesMutualInformationImageToImageMetricFor3DBLUTFFD computes the mutual
66 * information between a fixed and moving image to be registered.
68 * This class is templated over the FixedImage type and the MovingImage
71 * The fixed and moving images are set via methods SetFixedImage() and
72 * SetMovingImage(). This metric makes use of user specified Transform and
73 * Interpolator. The Transform is used to map points from the fixed image to
74 * the moving image domain. The Interpolator is used to evaluate the image
75 * intensity at user specified geometric points in the moving image.
76 * The Transform and Interpolator are set via methods SetTransform() and
79 * If a BSplineInterpolationFunction is used, this class obtain
80 * image derivatives from the BSpline interpolator. Otherwise,
81 * image derivatives are computed using central differencing.
83 * \warning This metric assumes that the moving image has already been
84 * connected to the interpolator outside of this class.
86 * The method GetValue() computes of the mutual information
87 * while method GetValueAndDerivative() computes
88 * both the mutual information and its derivatives with respect to the
89 * transform parameters.
91 * The calculations are based on the method of Mattes et al [1,2]
92 * where the probability density distribution are estimated using
93 * Parzen histograms. Since the fixed image PDF does not contribute
94 * to the derivatives, it does not need to be smooth. Hence,
95 * a zero order (box car) BSpline kernel is used
96 * for the fixed image intensity PDF. On the other hand, to ensure
97 * smoothness a third order BSpline kernel is used for the
98 * moving image intensity PDF.
100 * On Initialize(), the FixedImage is uniformly sampled within
101 * the FixedImageRegion. The number of samples used can be set
102 * via SetNumberOfSpatialSamples(). Typically, the number of
103 * spatial samples used should increase with the image size.
105 * The option UseAllPixelOn() disables the random sampling and uses
106 * all the pixels of the FixedImageRegion in order to estimate the
107 * joint intensity PDF.
109 * During each call of GetValue(), GetDerivatives(),
110 * GetValueAndDerivatives(), marginal and joint intensity PDF's
111 * values are estimated at discrete position or bins.
112 * The number of bins used can be set via SetNumberOfHistogramBins().
113 * To handle data with arbitray magnitude and dynamic range,
114 * the image intensity is scale such that any contribution to the
115 * histogram will fall into a valid bin.
117 * One the PDF's have been contructed, the mutual information
118 * is obtained by doubling summing over the discrete PDF values.
122 * 1. This class returns the negative mutual information value.
123 * 2. This class in not thread safe due the private data structures
124 * used to the store the sampled points and the marginal and joint pdfs.
127 * [1] "Nonrigid multimodality image registration"
128 * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
129 * Medical Imaging 2001: Image Processing, 2001, pp. 1609-1620.
130 * [2] "PET-CT Image Registration in the Chest Using Free-form Deformations"
131 * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
132 * IEEE Transactions in Medical Imaging. Vol.22, No.1,
133 January 2003. pp.120-128.
134 * [3] "Optimization of Mutual Information for MultiResolution Image
136 * P. Thevenaz and M. Unser
137 * IEEE Transactions in Image Processing, 9(12) December 2000.
139 * \ingroup RegistrationMetrics
140 * \ingroup ThreadUnSafe
142 template <class TFixedImage,class TMovingImage >
143 class ITK_EXPORT MattesMutualInformationImageToImageMetricFor3DBLUTFFD :
144 public ImageToImageMetric< TFixedImage, TMovingImage >
148 /** Standard class typedefs. */
149 typedef MattesMutualInformationImageToImageMetricFor3DBLUTFFD Self;
150 typedef ImageToImageMetric< TFixedImage, TMovingImage > Superclass;
151 typedef SmartPointer<Self> Pointer;
152 typedef SmartPointer<const Self> ConstPointer;
154 /** Method for creation through the object factory. */
157 /** Run-time type information (and related methods). */
158 itkTypeMacro(MattesMutualInformationImageToImageMetricFor3DBLUTFFD, ImageToImageMetric);
160 /** Types inherited from Superclass. */
161 typedef typename Superclass::TransformType TransformType;
162 typedef typename Superclass::TransformPointer TransformPointer;
163 typedef typename Superclass::TransformJacobianType TransformJacobianType;
164 typedef typename Superclass::InterpolatorType InterpolatorType;
165 typedef typename Superclass::MeasureType MeasureType;
166 typedef typename Superclass::DerivativeType DerivativeType;
167 typedef typename Superclass::ParametersType ParametersType;
168 typedef typename Superclass::FixedImageType FixedImageType;
169 typedef typename Superclass::MovingImageType MovingImageType;
170 typedef typename Superclass::FixedImageConstPointer FixedImageConstPointer;
171 typedef typename Superclass::MovingImageConstPointer MovingImageCosntPointer;
172 typedef typename Superclass::InputPointType InputPointType;
173 typedef typename Superclass::OutputPointType OutputPointType;
175 typedef typename Superclass::CoordinateRepresentationType
176 CoordinateRepresentationType;
178 /** Index and Point typedef support. */
179 typedef typename FixedImageType::IndexType FixedImageIndexType;
180 typedef typename FixedImageIndexType::IndexValueType FixedImageIndexValueType;
181 typedef typename MovingImageType::IndexType MovingImageIndexType;
182 typedef typename TransformType::InputPointType FixedImagePointType;
183 typedef typename TransformType::OutputPointType MovingImagePointType;
185 /** The moving image dimension. */
186 itkStaticConstMacro( MovingImageDimension, unsigned int,
187 MovingImageType::ImageDimension );
190 * Initialize the Metric by
191 * (1) making sure that all the components are present and plugged
192 * together correctly,
193 * (2) uniformly select NumberOfSpatialSamples within
194 * the FixedImageRegion, and
195 * (3) allocate memory for pdf data structures. */
196 virtual void Initialize(void) throw ( ExceptionObject );
198 /** Get the derivatives of the match measure. */
199 void GetDerivative( const ParametersType& parameters,
200 DerivativeType & Derivative ) const;
202 /** Get the value. */
203 MeasureType GetValue( const ParametersType& parameters ) const;
205 /** Get the value and derivatives for single valued optimizers. */
206 void GetValueAndDerivative( const ParametersType& parameters,
208 DerivativeType& Derivative ) const;
210 /** Number of spatial samples to used to compute metric */
211 itkSetClampMacro( NumberOfSpatialSamples, unsigned long,
212 1, NumericTraits<unsigned long>::max() );
213 itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long);
215 /** Number of bins to used in the histogram. Typical value is 50. */
216 itkSetClampMacro( NumberOfHistogramBins, unsigned long,
217 1, NumericTraits<unsigned long>::max() );
218 itkGetConstReferenceMacro( NumberOfHistogramBins, unsigned long);
220 /** Reinitialize the seed of the random number generator that selects the
221 * sample of pixels used for estimating the image histograms and the joint
222 * histogram. By nature, this metric is not deterministic, since at each run
223 * it may select a different set of pixels. By initializing the random number
224 * generator seed to the same value you can restore determinism. On the other
225 * hand, calling the method ReinitializeSeed() without arguments will use the
226 * clock from your machine in order to have a very random initialization of
227 * the seed. This will indeed increase the non-deterministic behavior of the
229 void ReinitializeSeed();
230 void ReinitializeSeed(int);
232 /** Select whether the metric will be computed using all the pixels on the
233 * fixed image region, or only using a set of randomly selected pixels. */
234 itkSetMacro(UseAllPixels,bool);
235 itkGetConstReferenceMacro(UseAllPixels,bool);
236 itkBooleanMacro(UseAllPixels);
238 /** This variable selects the method to be used for computing the Metric
239 * derivatives with respect to the Transform parameters. Two modes of
240 * computation are available. The choice between one and the other is a
241 * trade-off between computation speed and memory allocations. The two modes
242 * are described in detail below:
244 * UseExplicitPDFDerivatives = True
245 * will compute the Metric derivative by first calculating the derivatives of
246 * each one of the Joint PDF bins with respect to each one of the Transform
247 * parameters and then accumulating these contributions in the final metric
248 * derivative array by using a bin-specific weight. The memory required for
249 * storing the intermediate derivatives is a 3D array of doubles with size
250 * equals to the product of (number of histogram bins)^2 times number of
251 * transform parameters. This method is well suited for Transform with a small
252 * number of parameters.
254 * UseExplicitPDFDerivatives = False will compute the Metric derivative by
255 * first computing the weights for each one of the Joint PDF bins and caching
256 * them into an array. Then it will revisit each one of the PDF bins for
257 * computing its weighted contribution to the full derivative array. In this
258 * method an extra 2D array is used for storing the weights of each one of
259 * the PDF bins. This is an array of doubles with size equals to (number of
260 * histogram bins)^2. This method is well suited for Transforms with a large
261 * number of parameters, such as, BSplineDeformableTransforms. */
262 itkSetMacro(UseExplicitPDFDerivatives,bool);
263 itkGetConstReferenceMacro(UseExplicitPDFDerivatives,bool);
264 itkBooleanMacro(UseExplicitPDFDerivatives);
266 /** This boolean flag is only relevant when this metric is used along
267 * with a BSplineDeformableTransform. The flag enables/disables the
268 * caching of values computed when a physical point is mapped through
269 * the BSplineDeformableTransform. In particular it will cache the
270 * values of the BSpline weights for that points, and the indexes
271 * indicating what BSpline-grid nodes are relevant for that specific
272 * point. This caching is made optional due to the fact that the
273 * memory arrays used for the caching can reach large sizes even for
274 * moderate image size problems. For example, for a 3D image of
275 * 256^3, using 20% of pixels, these arrays will take about 1
276 * Gigabyte of RAM for storage. The ratio of computing time between
277 * using the cache or not using the cache can reach 1:5, meaning that
278 * using the caching can provide a five times speed up. It is
279 * therefore, interesting to enable the caching, if enough memory is
280 * available for it. The caching is enabled by default, in order to
281 * preserve backward compatibility with previous versions of ITK. */
282 itkSetMacro(UseCachingOfBSplineWeights,bool);
283 itkGetConstReferenceMacro(UseCachingOfBSplineWeights,bool);
284 itkBooleanMacro(UseCachingOfBSplineWeights);
288 MattesMutualInformationImageToImageMetricFor3DBLUTFFD();
289 virtual ~MattesMutualInformationImageToImageMetricFor3DBLUTFFD() {};
290 void PrintSelf(std::ostream& os, Indent indent) const;
292 /** \class FixedImageSpatialSample
293 * A fixed image spatial sample consists of the fixed domain point
294 * and the fixed image value at that point. */
296 class FixedImageSpatialSample
299 FixedImageSpatialSample():FixedImageValue(0.0) {
300 FixedImagePointValue.Fill(0.0);
302 ~FixedImageSpatialSample() {};
304 FixedImagePointType FixedImagePointValue;
305 double FixedImageValue;
306 unsigned int FixedImageParzenWindowIndex;
310 /** FixedImageSpatialSample typedef support. */
311 typedef std::vector<FixedImageSpatialSample>
312 FixedImageSpatialSampleContainer;
314 /** Container to store a set of points and fixed image values. */
315 FixedImageSpatialSampleContainer m_FixedImageSamples;
317 /** Uniformly select a sample set from the fixed image domain. */
318 virtual void SampleFixedImageDomain(
319 FixedImageSpatialSampleContainer& samples);
321 /** Gather all the pixels from the fixed image domain. */
322 virtual void SampleFullFixedImageDomain(
323 FixedImageSpatialSampleContainer& samples);
325 /** Transform a point from FixedImage domain to MovingImage domain.
326 * This function also checks if mapped point is within support region. */
327 virtual void TransformPoint( unsigned int sampleNumber,
328 const ParametersType& parameters,
329 MovingImagePointType& mappedPoint,
330 bool& sampleWithinSupportRegion,
331 double& movingImageValue ) const;
335 //purposely not implemented
336 MattesMutualInformationImageToImageMetricFor3DBLUTFFD(const Self&);
337 //purposely not implemented
338 void operator=(const Self&);
341 /** The marginal PDFs are stored as std::vector. */
342 typedef float PDFValueType;
343 typedef std::vector<PDFValueType> MarginalPDFType;
345 /** The fixed image marginal PDF. */
346 mutable MarginalPDFType m_FixedImageMarginalPDF;
348 /** The moving image marginal PDF. */
349 mutable MarginalPDFType m_MovingImageMarginalPDF;
351 /** Helper array for storing the values of the JointPDF ratios. */
352 typedef double PRatioType;
353 typedef Array2D< PRatioType > PRatioArrayType;
354 mutable PRatioArrayType m_PRatioArray;
356 /** Helper variable for accumulating the derivative of the metric. */
357 mutable DerivativeType m_MetricDerivative;
359 /** Typedef for the joint PDF and PDF derivatives are stored as ITK Images. */
360 typedef Image<PDFValueType,2> JointPDFType;
361 typedef JointPDFType::IndexType JointPDFIndexType;
362 typedef JointPDFType::PixelType JointPDFValueType;
363 typedef JointPDFType::RegionType JointPDFRegionType;
364 typedef JointPDFType::SizeType JointPDFSizeType;
365 typedef Image<PDFValueType,3> JointPDFDerivativesType;
366 typedef JointPDFDerivativesType::IndexType JointPDFDerivativesIndexType;
367 typedef JointPDFDerivativesType::PixelType JointPDFDerivativesValueType;
368 typedef JointPDFDerivativesType::RegionType JointPDFDerivativesRegionType;
369 typedef JointPDFDerivativesType::SizeType JointPDFDerivativesSizeType;
371 /** The joint PDF and PDF derivatives. */
372 typename JointPDFType::Pointer m_JointPDF;
373 typename JointPDFDerivativesType::Pointer m_JointPDFDerivatives;
375 unsigned long m_NumberOfSpatialSamples;
376 unsigned long m_NumberOfParameters;
378 /** Variables to define the marginal and joint histograms. */
379 unsigned long m_NumberOfHistogramBins;
380 double m_MovingImageNormalizedMin;
381 double m_FixedImageNormalizedMin;
382 double m_MovingImageTrueMin;
383 double m_MovingImageTrueMax;
384 double m_FixedImageBinSize;
385 double m_MovingImageBinSize;
387 /** Typedefs for BSpline kernel and derivative functions. */
388 typedef BSplineKernelFunction<3> CubicBSplineFunctionType;
389 typedef BSplineDerivativeKernelFunction<3> CubicBSplineDerivativeFunctionType;
391 /** Cubic BSpline kernel for computing Parzen histograms. */
392 typename CubicBSplineFunctionType::Pointer m_CubicBSplineKernel;
393 typename CubicBSplineDerivativeFunctionType::Pointer
394 m_CubicBSplineDerivativeKernel;
396 /** Precompute fixed image parzen window indices. */
397 virtual void ComputeFixedImageParzenWindowIndices(
398 FixedImageSpatialSampleContainer& samples );
401 * Types and variables related to image derivative calculations.
402 * If a BSplineInterpolationFunction is used, this class obtain
403 * image derivatives from the BSpline interpolator. Otherwise,
404 * image derivatives are computed using central differencing.
406 typedef CovariantVector< double,
407 itkGetStaticConstMacro(MovingImageDimension) >
408 ImageDerivativesType;
410 /** Compute image derivatives at a point. */
411 virtual void ComputeImageDerivatives( const MovingImagePointType& mappedPoint,
412 ImageDerivativesType& gradient ) const;
414 /** Boolean to indicate if the interpolator BSpline. */
415 bool m_InterpolatorIsBSpline;
417 /** Typedefs for using BSpline interpolator. */
419 BSplineInterpolateImageFunction<MovingImageType,
420 CoordinateRepresentationType>
421 BSplineInterpolatorType;
423 /** Pointer to BSplineInterpolator. */
424 typename BSplineInterpolatorType::Pointer m_BSplineInterpolator;
426 /** Typedefs for using central difference calculator. */
427 typedef CentralDifferenceImageFunction<MovingImageType,
428 CoordinateRepresentationType>
429 DerivativeFunctionType;
431 /** Pointer to central difference calculator. */
432 typename DerivativeFunctionType::Pointer m_DerivativeCalculator;
435 /** Compute PDF derivative contribution for each parameter. */
436 virtual void ComputePDFDerivatives( unsigned int sampleNumber,
437 int movingImageParzenWindowIndex,
438 const ImageDerivativesType&
439 movingImageGradientValue,
440 double cubicBSplineDerivativeValue
444 * Types and variables related to BSpline deformable transforms.
445 * If the transform is of type third order BSplineDeformableTransform,
446 * then we can speed up the metric derivative calculation by
447 * only inspecting the parameters within the support region
451 /** Boolean to indicate if the transform is BSpline deformable. */
452 bool m_TransformIsBSpline;
454 /** The number of BSpline parameters per image dimension. */
455 long m_NumParametersPerDim;
458 * The number of BSpline transform weights is the number of
459 * of parameter in the support region (per dimension ). */
460 unsigned long m_NumBSplineWeights;
462 /** The fixed image dimension. */
463 itkStaticConstMacro( FixedImageDimension, unsigned int,
464 FixedImageType::ImageDimension );
467 * Enum of the deformabtion field spline order.
469 enum { DeformationSplineOrder = 3 };
472 * Typedefs for the BSplineDeformableTransform.
474 typedef BSplineDeformableTransform<
475 CoordinateRepresentationType,
476 ::itk::GetImageDimension<FixedImageType>::ImageDimension,
477 DeformationSplineOrder> BSplineTransformType;
478 typedef typename BSplineTransformType::WeightsType
479 BSplineTransformWeightsType;
480 typedef typename BSplineTransformType::ParameterIndexArrayType
481 BSplineTransformIndexArrayType;
484 * Variables used when transform is of type BSpline deformable.
486 typename BSplineTransformType::Pointer m_BSplineTransform;
489 * Cache pre-transformed points, weights, indices and
490 * within support region flag.
492 typedef typename BSplineTransformWeightsType::ValueType WeightsValueType;
493 typedef Array2D<WeightsValueType> BSplineTransformWeightsArrayType;
494 typedef typename BSplineTransformIndexArrayType::ValueType IndexValueType;
495 typedef Array2D<IndexValueType> BSplineTransformIndicesArrayType;
496 typedef std::vector<MovingImagePointType> MovingImagePointArrayType;
497 typedef std::vector<bool> BooleanArrayType;
499 BSplineTransformWeightsArrayType m_BSplineTransformWeightsArray;
500 BSplineTransformIndicesArrayType m_BSplineTransformIndicesArray;
501 MovingImagePointArrayType m_PreTransformPointsArray;
502 BooleanArrayType m_WithinSupportRegionArray;
504 typedef FixedArray<unsigned long,
505 ::itk::GetImageDimension<FixedImageType>::ImageDimension>
506 ParametersOffsetType;
507 ParametersOffsetType m_ParametersOffset;
511 virtual void PreComputeTransformValues();
513 bool m_ReseedIterator;
516 // Selection of explicit or implicit computation of PDF derivatives
517 // with respect to Transform parameters.
518 bool m_UseExplicitPDFDerivatives;
520 // Variables needed for optionally caching values when using a BSpline transform.
521 bool m_UseCachingOfBSplineWeights;
522 mutable BSplineTransformWeightsType m_Weights;
523 mutable BSplineTransformIndexArrayType m_Indices;
527 } // end namespace itk
529 #ifndef ITK_MANUAL_INSTANTIATION
530 #include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx"