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