]> Creatis software - clitk.git/blob - registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
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
8
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.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18
19 /*=========================================================================
20
21   Program:   Insight Segmentation & Registration Toolkit
22   Module:    $RCSfile: itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h,v $
23   Language:  C++
24   Date:      $Date: 2010/06/14 17:32:07 $
25   Version:   $Revision: 1.1 $
26
27   Copyright (c) Insight Software Consortium. All rights reserved.
28   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
29
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.
33
34 =========================================================================*/
35 #ifndef __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_h
36 #define __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_h
37
38 #if ITK_VERSION_MAJOR >= 4
39   #include "itkImageToImageMetric.h"
40 #else
41   #include "itkOptImageToImageMetric.h"
42 #endif
43 #include "itkCovariantVector.h"
44 #include "itkPoint.h"
45 #include "itkIndex.h"
46 #include "itkBSplineKernelFunction.h"
47 #include "itkBSplineDerivativeKernelFunction.h"
48 #include "itkCentralDifferenceImageFunction.h"
49 #include "itkBSplineInterpolateImageFunction.h"
50 #include "itkBSplineDeformableTransform.h"
51 #include "itkArray2D.h"
52
53 #include "itkMultiThreader.h"
54
55 namespace itk
56 {
57
58 /** \class MattesMutualInformationImageToImageMetricFor3DBLUTFFD
59  * \brief Computes the mutual information between two images to be
60  * registered using the method of Mattes et al.
61  *
62  * MattesMutualInformationImageToImageMetricFor3DBLUTFFD computes the mutual
63  * information between a fixed and moving image to be registered.
64  *
65  * This class is templated over the FixedImage type and the MovingImage
66  * type.
67  *
68  * The fixed and moving images are set via methods SetFixedImage() and
69  * SetMovingImage(). This metric makes use of user specified Transform and
70  * Interpolator. The Transform is used to map points from the fixed image to
71  * the moving image domain. The Interpolator is used to evaluate the image
72  * intensity at user specified geometric points in the moving image.
73  * The Transform and Interpolator are set via methods SetTransform() and
74  * SetInterpolator().
75  *
76  * If a BSplineInterpolationFunction is used, this class obtain
77  * image derivatives from the BSpline interpolator. Otherwise,
78  * image derivatives are computed using central differencing.
79  *
80  * \warning This metric assumes that the moving image has already been
81  * connected to the interpolator outside of this class.
82  *
83  * The method GetValue() computes of the mutual information
84  * while method GetValueAndDerivative() computes
85  * both the mutual information and its derivatives with respect to the
86  * transform parameters.
87  *
88  * The calculations are based on the method of Mattes et al [1,2]
89  * where the probability density distribution are estimated using
90  * Parzen histograms. Since the fixed image PDF does not contribute
91  * to the derivatives, it does not need to be smooth. Hence,
92  * a zero order (box car) BSpline kernel is used
93  * for the fixed image intensity PDF. On the other hand, to ensure
94  * smoothness a third order BSpline kernel is used for the
95  * moving image intensity PDF.
96  *
97  * On Initialize(), the FixedImage is uniformly sampled within
98  * the FixedImageRegion. The number of samples used can be set
99  * via SetNumberOfSpatialSamples(). Typically, the number of
100  * spatial samples used should increase with the image size.
101  *
102  * The option UseAllPixelOn() disables the random sampling and uses
103  * all the pixels of the FixedImageRegion in order to estimate the
104  * joint intensity PDF.
105  *
106  * During each call of GetValue(), GetDerivatives(),
107  * GetValueAndDerivatives(), marginal and joint intensity PDF's
108  * values are estimated at discrete position or bins.
109  * The number of bins used can be set via SetNumberOfHistogramBins().
110  * To handle data with arbitray magnitude and dynamic range,
111  * the image intensity is scale such that any contribution to the
112  * histogram will fall into a valid bin.
113  *
114  * One the PDF's have been contructed, the mutual information
115  * is obtained by doubling summing over the discrete PDF values.
116  *
117  *
118  * Notes:
119  * 1. This class returns the negative mutual information value.
120  * 2. This class in not thread safe due the private data structures
121  *     used to the store the sampled points and the marginal and joint pdfs.
122  *
123  * References:
124  * [1] "Nonrigid multimodality image registration"
125  *      D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
126  *      Medical Imaging 2001: Image Processing, 2001, pp. 1609-1620.
127  * [2] "PET-CT Image Registration in the Chest Using Free-form Deformations"
128  *      D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
129  *      IEEE Transactions in Medical Imaging. Vol.22, No.1,
130         January 2003. pp.120-128.
131  * [3] "Optimization of Mutual Information for MultiResolution Image
132  *      Registration"
133  *      P. Thevenaz and M. Unser
134  *      IEEE Transactions in Image Processing, 9(12) December 2000.
135  *
136  * \ingroup RegistrationMetrics
137  */
138 template <class TFixedImage,class TMovingImage >
139 class ITK_EXPORT MattesMutualInformationImageToImageMetricFor3DBLUTFFD :
140   public ImageToImageMetric< TFixedImage, TMovingImage >
141 {
142 public:
143
144   /** Standard class typedefs. */
145   typedef MattesMutualInformationImageToImageMetricFor3DBLUTFFD           Self;
146   typedef ImageToImageMetric< TFixedImage, TMovingImage >     Superclass;
147   typedef SmartPointer<Self>                                  Pointer;
148   typedef SmartPointer<const Self>                            ConstPointer;
149
150   /** Method for creation through the object factory. */
151   itkNewMacro(Self);
152
153   /** Run-time type information (and related methods). */
154   itkTypeMacro(MattesMutualInformationImageToImageMetricFor3DBLUTFFD,
155                ImageToImageMetric);
156
157   /** Types inherited from Superclass. */
158   typedef typename Superclass::TransformType            TransformType;
159   typedef typename Superclass::TransformPointer         TransformPointer;
160   typedef typename Superclass::TransformJacobianType    TransformJacobianType;
161   typedef typename Superclass::InterpolatorType         InterpolatorType;
162   typedef typename Superclass::MeasureType              MeasureType;
163   typedef typename Superclass::DerivativeType           DerivativeType;
164   typedef typename Superclass::ParametersType           ParametersType;
165   typedef typename Superclass::FixedImageType           FixedImageType;
166   typedef typename Superclass::MovingImageType          MovingImageType;
167   typedef typename Superclass::MovingImagePointType     MovingImagePointType;
168   typedef typename Superclass::FixedImageConstPointer   FixedImageConstPointer;
169   typedef typename Superclass::MovingImageConstPointer  MovingImageConstPointer;
170   typedef typename Superclass::BSplineTransformWeightsType
171   BSplineTransformWeightsType;
172   typedef typename Superclass::BSplineTransformIndexArrayType
173   BSplineTransformIndexArrayType;
174
175   typedef typename Superclass::CoordinateRepresentationType
176   CoordinateRepresentationType;
177   typedef typename Superclass::FixedImageSampleContainer
178   FixedImageSampleContainer;
179   typedef typename Superclass::ImageDerivativesType     ImageDerivativesType;
180   typedef typename Superclass::WeightsValueType         WeightsValueType;
181   typedef typename Superclass::IndexValueType           IndexValueType;
182
183   /** The moving image dimension. */
184   itkStaticConstMacro( MovingImageDimension, unsigned int,
185                        MovingImageType::ImageDimension );
186
187   /**
188    *  Initialize the Metric by
189    *  (1) making sure that all the components are present and plugged
190    *      together correctly,
191    *  (2) uniformly select NumberOfSpatialSamples within
192    *      the FixedImageRegion, and
193    *  (3) allocate memory for pdf data structures. */
194   virtual void Initialize(void) throw ( ExceptionObject );
195
196   /**  Get the value. */
197   MeasureType GetValue( const ParametersType & parameters ) const;
198
199   /** Get the derivatives of the match measure. */
200   void GetDerivative( const ParametersType & parameters,
201                       DerivativeType & Derivative ) const;
202
203   /**  Get the value and derivatives for single valued optimizers. */
204   void GetValueAndDerivative( const ParametersType & parameters,
205                               MeasureType & Value,
206                               DerivativeType & Derivative ) const;
207
208   /** Number of bins to used in the histogram. Typical value is 50. */
209   itkSetClampMacro( NumberOfHistogramBins, unsigned long,
210                     1, NumericTraits<unsigned long>::max() );
211   itkGetConstReferenceMacro( NumberOfHistogramBins, unsigned long);
212
213   /** This variable selects the method to be used for computing the Metric
214    * derivatives with respect to the Transform parameters. Two modes of
215    * computation are available. The choice between one and the other is a
216    * trade-off between computation speed and memory allocations. The two modes
217    * are described in detail below:
218    *
219    * UseExplicitPDFDerivatives = True
220    * will compute the Metric derivative by first calculating the derivatives of
221    * each one of the Joint PDF bins with respect to each one of the Transform
222    * parameters and then accumulating these contributions in the final metric
223    * derivative array by using a bin-specific weight.  The memory required for
224    * storing the intermediate derivatives is a 3D array of doubles with size
225    * equals to the product of (number of histogram bins)^2 times number of
226    * transform parameters. This method is well suited for Transform with a small
227    * number of parameters.
228    *
229    * UseExplicitPDFDerivatives = False will compute the Metric derivative by
230    * first computing the weights for each one of the Joint PDF bins and caching
231    * them into an array. Then it will revisit each one of the PDF bins for
232    * computing its weighted contribution to the full derivative array. In this
233    * method an extra 2D array is used for storing the weights of each one of
234    * the PDF bins. This is an array of doubles with size equals to (number of
235    * histogram bins)^2. This method is well suited for Transforms with a large
236    * number of parameters, such as, BSplineDeformableTransforms. */
237   itkSetMacro(UseExplicitPDFDerivatives,bool);
238   itkGetConstReferenceMacro(UseExplicitPDFDerivatives,bool);
239   itkBooleanMacro(UseExplicitPDFDerivatives);
240
241 protected:
242
243   MattesMutualInformationImageToImageMetricFor3DBLUTFFD();
244   virtual ~MattesMutualInformationImageToImageMetricFor3DBLUTFFD();
245   void PrintSelf(std::ostream& os, Indent indent) const;
246
247 private:
248
249   //purposely not implemented
250   MattesMutualInformationImageToImageMetricFor3DBLUTFFD(const Self &);
251   //purposely not implemented
252   void operator=(const Self &);
253
254
255   /** The marginal PDFs are stored as std::vector. */
256   typedef float                       PDFValueType;
257   typedef float *                     MarginalPDFType;
258
259   mutable MarginalPDFType             m_FixedImageMarginalPDF;
260
261   /** The moving image marginal PDF. */
262   mutable MarginalPDFType             m_MovingImageMarginalPDF;
263
264   /** Helper array for storing the values of the JointPDF ratios. */
265   typedef double                      PRatioType;
266   typedef Array2D< PRatioType >       PRatioArrayType;
267   mutable PRatioArrayType             m_PRatioArray;
268
269   /** Helper variable for accumulating the derivative of the metric. */
270   mutable DerivativeType              m_MetricDerivative;
271   mutable DerivativeType            * m_ThreaderMetricDerivative;
272
273   /** Typedef for the joint PDF and PDF derivatives are stored as ITK Images. */
274   typedef Image<PDFValueType,2>                 JointPDFType;
275   typedef Image<PDFValueType,3>                 JointPDFDerivativesType;
276   typedef JointPDFType::IndexType               JointPDFIndexType;
277   typedef JointPDFType::PixelType               JointPDFValueType;
278   typedef JointPDFType::RegionType              JointPDFRegionType;
279   typedef JointPDFType::SizeType                JointPDFSizeType;
280   typedef JointPDFDerivativesType::IndexType    JointPDFDerivativesIndexType;
281   typedef JointPDFDerivativesType::PixelType    JointPDFDerivativesValueType;
282   typedef JointPDFDerivativesType::RegionType   JointPDFDerivativesRegionType;
283   typedef JointPDFDerivativesType::SizeType     JointPDFDerivativesSizeType;
284
285   /** The joint PDF and PDF derivatives. */
286   typename JointPDFType::Pointer                m_JointPDF;
287   unsigned long                                 m_JointPDFBufferSize;
288   typename JointPDFDerivativesType::Pointer     m_JointPDFDerivatives;
289   unsigned long                                 m_JointPDFDerivativesBufferSize;
290
291   /** Variables to define the marginal and joint histograms. */
292   unsigned long  m_NumberOfHistogramBins;
293   double         m_MovingImageNormalizedMin;
294   double         m_FixedImageNormalizedMin;
295   double         m_FixedImageTrueMin;
296   double         m_FixedImageTrueMax;
297   double         m_MovingImageTrueMin;
298   double         m_MovingImageTrueMax;
299   double         m_FixedImageBinSize;
300   double         m_MovingImageBinSize;
301
302   /** Typedefs for BSpline kernel and derivative functions. */
303   typedef BSplineKernelFunction<3>           CubicBSplineFunctionType;
304   typedef BSplineDerivativeKernelFunction<3> CubicBSplineDerivativeFunctionType;
305
306   /** Cubic BSpline kernel for computing Parzen histograms. */
307   typename CubicBSplineFunctionType::Pointer   m_CubicBSplineKernel;
308   typename CubicBSplineDerivativeFunctionType::Pointer
309   m_CubicBSplineDerivativeKernel;
310
311   /** Precompute fixed image parzen window indices. */
312   virtual void ComputeFixedImageParzenWindowIndices(
313     FixedImageSampleContainer & samples );
314
315   /** Compute PDF derivative contribution for each parameter. */
316   virtual void ComputePDFDerivatives( unsigned int threadID,
317                                       unsigned int sampleNumber,
318                                       int movingImageParzenWindowIndex,
319                                       const ImageDerivativesType
320                                       &  movingImageGradientValue,
321                                       double cubicBSplineDerivativeValue
322                                     ) const;
323
324   PDFValueType                                * m_ThreaderFixedImageMarginalPDF;
325   typename JointPDFType::Pointer              * m_ThreaderJointPDF;
326   typename JointPDFDerivativesType::Pointer   * m_ThreaderJointPDFDerivatives;
327   int                                         * m_ThreaderJointPDFStartBin;
328   int                                         * m_ThreaderJointPDFEndBin;
329   mutable double                              * m_ThreaderJointPDFSum;
330   mutable double                                m_JointPDFSum;
331
332   bool                                          m_UseExplicitPDFDerivatives;
333   mutable bool                                  m_ImplicitDerivativesSecondPass;
334
335
336   virtual inline void GetValueThreadPreProcess( unsigned int threadID,
337       bool withinSampleThread ) const;
338   virtual inline bool GetValueThreadProcessSample( unsigned int threadID,
339       unsigned long fixedImageSample,
340       const MovingImagePointType & mappedPoint,
341       double movingImageValue ) const;
342   virtual inline void GetValueThreadPostProcess( unsigned int threadID,
343       bool withinSampleThread ) const;
344
345   virtual inline void GetValueAndDerivativeThreadPreProcess(
346     unsigned int threadID,
347     bool withinSampleThread ) const;
348   virtual inline bool GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
349       unsigned long fixedImageSample,
350       const MovingImagePointType & mappedPoint,
351       double movingImageValue,
352       const ImageDerivativesType &
353       movingImageGradientValue ) const;
354   virtual inline void GetValueAndDerivativeThreadPostProcess(
355     unsigned int threadID,
356     bool withinSampleThread ) const;
357
358 };
359
360 } // end namespace itk
361
362 #ifndef ITK_MANUAL_INSTANTIATION
363 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx"
364 #endif
365
366 #endif