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