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 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
 191   virtual void Initialize(void) ITK_OVERRIDE;
 192 #else
 193   virtual void Initialize(void) throw ( ExceptionObject ) ITK_OVERRIDE;
 194 #endif
 195
 196   /**  Get the value. */
 197   MeasureType GetValue( const ParametersType & parameters ) const ITK_OVERRIDE;
 198
 199   /** Get the derivatives of the match measure. */
 200   void GetDerivative( const ParametersType & parameters,
 201                       DerivativeType & Derivative ) const ITK_OVERRIDE;
 202
 203   /**  Get the value and derivatives for single valued optimizers. */
 204   void GetValueAndDerivative( const ParametersType & parameters,
 205                               MeasureType & Value,
 206                               DerivativeType & Derivative ) const ITK_OVERRIDE;
 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 ITK_OVERRIDE;
 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 ITK_OVERRIDE;
 338   virtual inline bool GetValueThreadProcessSample( unsigned int threadID,
 339       unsigned long fixedImageSample,
 340       const MovingImagePointType & mappedPoint,
 341       double movingImageValue ) const ITK_OVERRIDE;
 342   virtual inline void GetValueThreadPostProcess( unsigned int threadID,
 343       bool withinSampleThread ) const ITK_OVERRIDE;
 344
 345   virtual inline void GetValueAndDerivativeThreadPreProcess(
 346     unsigned int threadID,
 347     bool withinSampleThread ) const ITK_OVERRIDE;
 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 ITK_OVERRIDE;
 354   virtual inline void GetValueAndDerivativeThreadPostProcess(
 355     unsigned int threadID,
 356     bool withinSampleThread ) const ITK_OVERRIDE;
 357
 358 };
 359
 360 } // end namespace itk
 361
 362 #ifndef ITK_MANUAL_INSTANTIATION
 363 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx"
 364 #endif
 365
 366 #endif