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 #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) ITK_OVERRIDE;
 191
 192   /**  Get the value. */
 193   MeasureType GetValue( const ParametersType & parameters ) const ITK_OVERRIDE;
 194
 195   /** Get the derivatives of the match measure. */
 196   void GetDerivative( const ParametersType & parameters,
 197                       DerivativeType & Derivative ) const ITK_OVERRIDE;
 198
 199   /**  Get the value and derivatives for single valued optimizers. */
 200   void GetValueAndDerivative( const ParametersType & parameters,
 201                               MeasureType & Value,
 202                               DerivativeType & Derivative ) const ITK_OVERRIDE;
 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 ITK_OVERRIDE;
 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 ITK_OVERRIDE;
 334   virtual inline bool GetValueThreadProcessSample( unsigned int threadID,
 335       unsigned long fixedImageSample,
 336       const MovingImagePointType & mappedPoint,
 337       double movingImageValue ) const ITK_OVERRIDE;
 338   virtual inline void GetValueThreadPostProcess( unsigned int threadID,
 339       bool withinSampleThread ) const ITK_OVERRIDE;
 340
 341   virtual inline void GetValueAndDerivativeThreadPreProcess(
 342     unsigned int threadID,
 343     bool withinSampleThread ) const ITK_OVERRIDE;
 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 ITK_OVERRIDE;
 350   virtual inline void GetValueAndDerivativeThreadPostProcess(
 351     unsigned int threadID,
 352     bool withinSampleThread ) const ITK_OVERRIDE;
 353
 354 };
 355
 356 } // end namespace itk
 357
 358 #ifndef ITK_MANUAL_INSTANTIATION
 359 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx"
 360 #endif
 361
 362 #endif