X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkCorrelationRatioImageToImageMetric.txx;fp=registration%2FclitkCorrelationRatioImageToImageMetric.txx;h=595a347eff10748a7b277a9b9d3990831cf7e076;hb=c18059db4f507fd31b5898667f57eced7d48c5f7;hp=0000000000000000000000000000000000000000;hpb=6e6c4206bc13210d5dd5dabf1bd24f17de316a7a;p=clitk.git

diff --git a/registration/clitkCorrelationRatioImageToImageMetric.txx b/registration/clitkCorrelationRatioImageToImageMetric.txx
new file mode 100644
index 0000000..595a347
--- /dev/null
+++ b/registration/clitkCorrelationRatioImageToImageMetric.txx
@@ -0,0 +1,445 @@
+/*=========================================================================
+  Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
+
+  Authors belong to:
+  - University of LYON              http://www.universite-lyon.fr/
+  - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
+  - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
+
+  This software is distributed WITHOUT ANY WARRANTY; without even
+  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+  PURPOSE.  See the copyright notices for more information.
+
+  It is distributed under dual licence
+
+  - BSD        See included LICENSE.txt file
+  - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+======================================================================-====*/
+
+#ifndef _clitkCorrelationRatioImageToImageMetric_txx
+#define _clitkCorrelationRatioImageToImageMetric_txx
+
+/**
+ * @file   clitkCorrelationRatioImageToImageMetric.txx
+ * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
+ * @date   July 30  18:14:53 2007
+ *
+ * @brief  Compute the correlation ratio between 2 images
+ *
+ *
+ */
+
+#include "clitkCorrelationRatioImageToImageMetric.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionIterator.h"
+
+namespace clitk
+{
+
+/*
+ * Constructor
+ */
+template <class TFixedImage, class TMovingImage>
+CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
+::CorrelationRatioImageToImageMetric()
+{
+  m_NumberOfBins = 50;
+
+}
+
+template <class TFixedImage, class TMovingImage>
+void
+CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
+::Initialize(void) throw ( ExceptionObject )
+{
+
+  this->Superclass::Initialize();
+
+  // Compute the minimum and maximum for the FixedImage over the FixedImageRegion.
+  // We can't use StatisticsImageFilter to do this because the filter computes the min/max for the largest possible region
+  double fixedImageMin = NumericTraits<double>::max();
+  double fixedImageMax = NumericTraits<double>::NonpositiveMin();
+
+  typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
+  FixedIteratorType fixedImageIterator(
+    this->m_FixedImage, this->GetFixedImageRegion() );
+
+  for ( fixedImageIterator.GoToBegin();
+        !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
+
+    double sample = static_cast<double>( fixedImageIterator.Get() );
+
+    if ( sample < fixedImageMin ) {
+      fixedImageMin = sample;
+    }
+
+    if ( sample > fixedImageMax ) {
+      fixedImageMax = sample;
+    }
+  }
+
+  // Compute binsize for the fixedImage
+  m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / m_NumberOfBins;
+  m_FixedImageMin=fixedImageMin;
+  //Allocate mempry and initialise the fixed image bin
+  m_NumberOfPixelsCountedPerBin.resize( m_NumberOfBins, 0 );
+  m_mMSVPB.resize( m_NumberOfBins, 0.0 );
+  m_mSMVPB.resize( m_NumberOfBins, 0.0 );
+}
+
+
+/*
+ * Get the match Measure
+ */
+template <class TFixedImage, class TMovingImage>
+typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
+CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
+::GetValue( const TransformParametersType & parameters ) const
+{
+
+  itkDebugMacro("GetValue( " << parameters << " ) ");
+
+  FixedImageConstPointer fixedImage = this->m_FixedImage;
+
+  if( !fixedImage ) {
+    itkExceptionMacro( << "Fixed image has not been assigned" );
+  }
+
+  typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
+
+
+  FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
+
+  typename FixedImageType::IndexType index;
+
+  MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
+
+  this->m_NumberOfPixelsCounted = 0;
+  this->SetTransformParameters( parameters );
+
+
+  //temporary measures for the calculation
+  RealType mSMV=0;
+  RealType mMSV=0;
+
+  while(!ti.IsAtEnd()) {
+
+    index = ti.GetIndex();
+
+    typename Superclass::InputPointType inputPoint;
+    fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
+
+    // Verify that the point is in the fixed Image Mask
+    if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
+      ++ti;
+      continue;
+    }
+
+    typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
+
+    //Verify that the point is in the moving Image Mask
+    if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
+      ++ti;
+      continue;
+    }
+
+    // Verify is the interpolated value is in the buffer
+    if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
+      //Accumulate calculations for the correlation ratio
+      //For each pixel the is in both masks and the buffer we adapt the following measures:
+      //movingMeanSquaredValue mMSV; movingSquaredMeanValue mSMV;
+      //movingMeanSquaredValuePerBin[i] mSMVPB; movingSquaredMeanValuePerBin[i] mSMVPB
+      //NumberOfPixelsCounted, NumberOfPixelsCountedPerBin[i]
+
+      //get the value of the moving image
+      const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
+      // for the variance of the overlapping moving image we accumulate the following measures
+      const RealType movingSquaredValue=movingValue*movingValue;
+      mMSV+=movingSquaredValue;
+      mSMV+=movingValue;
+
+      //get the fixed value
+      const RealType fixedValue   = ti.Get();
+
+      //check in which bin the fixed value belongs, get the index
+      const double fixedImageBinTerm =        (fixedValue - m_FixedImageMin) / m_FixedImageBinSize;
+      const unsigned int fixedImageBinIndex = static_cast<unsigned int>( vcl_floor(fixedImageBinTerm ) );
+      //adapt the measures per bin
+      this->m_mMSVPB[fixedImageBinIndex]+=movingSquaredValue;
+      this->m_mSMVPB[fixedImageBinIndex]+=movingValue;
+      //increase the fixed image bin and the total pixel count
+      this->m_NumberOfPixelsCountedPerBin[fixedImageBinIndex]+=1;
+      this->m_NumberOfPixelsCounted++;
+    }
+
+    ++ti;
+  }
+
+  if( !this->m_NumberOfPixelsCounted ) {
+    itkExceptionMacro(<<"All the points mapped to outside of the moving image");
+  } else {
+
+    //apdapt the measures per bin
+    for (unsigned int i=0; i< m_NumberOfBins; i++ ) {
+      if (this->m_NumberOfPixelsCountedPerBin[i]>0) {
+        measure+=(this->m_mMSVPB[i]-((this->m_mSMVPB[i]*this->m_mSMVPB[i])/this->m_NumberOfPixelsCountedPerBin[i]));
+      }
+    }
+
+    //Normalize with the global measures
+    measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
+    return measure;
+
+  }
+}
+
+
+
+
+
+/*
+ * Get the Derivative Measure
+ */
+template < class TFixedImage, class TMovingImage>
+void
+CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
+::GetDerivative( const TransformParametersType & parameters,
+                 DerivativeType & derivative  ) const
+{
+
+  itkDebugMacro("GetDerivative( " << parameters << " ) ");
+
+  if( !this->GetGradientImage() ) {
+    itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
+  }
+
+  FixedImageConstPointer fixedImage = this->m_FixedImage;
+
+  if( !fixedImage ) {
+    itkExceptionMacro( << "Fixed image has not been assigned" );
+  }
+
+  const unsigned int ImageDimension = FixedImageType::ImageDimension;
+
+
+  typedef  itk::ImageRegionConstIteratorWithIndex<
+  FixedImageType> FixedIteratorType;
+
+  typedef  itk::ImageRegionConstIteratorWithIndex<
+  ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
+
+
+  FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
+
+  typename FixedImageType::IndexType index;
+
+  this->m_NumberOfPixelsCounted = 0;
+
+  this->SetTransformParameters( parameters );
+
+  const unsigned int ParametersDimension = this->GetNumberOfParameters();
+  derivative = DerivativeType( ParametersDimension );
+  derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
+
+  ti.GoToBegin();
+
+  while(!ti.IsAtEnd()) {
+
+    index = ti.GetIndex();
+
+    typename Superclass::InputPointType inputPoint;
+    fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
+
+    if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
+      ++ti;
+      continue;
+    }
+
+    typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
+
+    if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
+      ++ti;
+      continue;
+    }
+
+    if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
+      const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
+
+      const TransformJacobianType & jacobian =
+        this->m_Transform->GetJacobian( inputPoint );
+
+
+      const RealType fixedValue     = ti.Value();
+      this->m_NumberOfPixelsCounted++;
+      const RealType diff = movingValue - fixedValue;
+
+      // Get the gradient by NearestNeighboorInterpolation:
+      // which is equivalent to round up the point components.
+      typedef typename Superclass::OutputPointType OutputPointType;
+      typedef typename OutputPointType::CoordRepType CoordRepType;
+      typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
+      MovingImageContinuousIndexType;
+
+      MovingImageContinuousIndexType tempIndex;
+      this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
+
+      typename MovingImageType::IndexType mappedIndex;
+      for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
+        mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
+      }
+
+      const GradientPixelType gradient =
+        this->GetGradientImage()->GetPixel( mappedIndex );
+
+      for(unsigned int par=0; par<ParametersDimension; par++) {
+        RealType sum = NumericTraits< RealType >::Zero;
+        for(unsigned int dim=0; dim<ImageDimension; dim++) {
+          sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
+        }
+        derivative[par] += sum;
+      }
+    }
+
+    ++ti;
+  }
+
+  if( !this->m_NumberOfPixelsCounted ) {
+    itkExceptionMacro(<<"All the points mapped to outside of the moving image");
+  } else {
+    for(unsigned int i=0; i<ParametersDimension; i++) {
+      derivative[i] /= this->m_NumberOfPixelsCounted;
+    }
+  }
+
+}
+
+
+/*
+ * Get both the match Measure and the Derivative Measure
+ */
+template <class TFixedImage, class TMovingImage>
+void
+CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
+::GetValueAndDerivative(const TransformParametersType & parameters,
+                        MeasureType & value, DerivativeType  & derivative) const
+{
+
+  itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
+
+  if( !this->GetGradientImage() ) {
+    itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
+  }
+
+  FixedImageConstPointer fixedImage = this->m_FixedImage;
+
+  if( !fixedImage ) {
+    itkExceptionMacro( << "Fixed image has not been assigned" );
+  }
+
+  const unsigned int ImageDimension = FixedImageType::ImageDimension;
+
+  typedef  itk::ImageRegionConstIteratorWithIndex<
+  FixedImageType> FixedIteratorType;
+
+  typedef  itk::ImageRegionConstIteratorWithIndex<
+  ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
+
+
+  FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
+
+  typename FixedImageType::IndexType index;
+
+  MeasureType measure = NumericTraits< MeasureType >::Zero;
+
+  this->m_NumberOfPixelsCounted = 0;
+
+  this->SetTransformParameters( parameters );
+
+  const unsigned int ParametersDimension = this->GetNumberOfParameters();
+  derivative = DerivativeType( ParametersDimension );
+  derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
+
+  ti.GoToBegin();
+
+  while(!ti.IsAtEnd()) {
+
+    index = ti.GetIndex();
+
+    typename Superclass::InputPointType inputPoint;
+    fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
+
+    if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
+      ++ti;
+      continue;
+    }
+
+    typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
+
+    if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
+      ++ti;
+      continue;
+    }
+
+    if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
+      const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
+
+      const TransformJacobianType & jacobian =
+        this->m_Transform->GetJacobian( inputPoint );
+
+
+      const RealType fixedValue     = ti.Value();
+      this->m_NumberOfPixelsCounted++;
+
+      const RealType diff = movingValue - fixedValue;
+
+      measure += diff * diff;
+
+      // Get the gradient by NearestNeighboorInterpolation:
+      // which is equivalent to round up the point components.
+      typedef typename Superclass::OutputPointType OutputPointType;
+      typedef typename OutputPointType::CoordRepType CoordRepType;
+      typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
+      MovingImageContinuousIndexType;
+
+      MovingImageContinuousIndexType tempIndex;
+      this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
+
+      typename MovingImageType::IndexType mappedIndex;
+      for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
+        mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
+      }
+
+      const GradientPixelType gradient =
+        this->GetGradientImage()->GetPixel( mappedIndex );
+
+      for(unsigned int par=0; par<ParametersDimension; par++) {
+        RealType sum = NumericTraits< RealType >::Zero;
+        for(unsigned int dim=0; dim<ImageDimension; dim++) {
+          sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
+        }
+        derivative[par] += sum;
+      }
+    }
+
+    ++ti;
+  }
+
+  if( !this->m_NumberOfPixelsCounted ) {
+    itkExceptionMacro(<<"All the points mapped to outside of the moving image");
+  } else {
+    for(unsigned int i=0; i<ParametersDimension; i++) {
+      derivative[i] /= this->m_NumberOfPixelsCounted;
+    }
+    measure /= this->m_NumberOfPixelsCounted;
+  }
+
+  value = measure;
+
+}
+
+} // end namespace itk
+
+
+#endif
+