X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkCorrelationRatioImageToImageMetric.txx;h=2b45c2ffc8a84b72e25d3e2f812ae3810b0092e2;hb=65315d4598295e8581a22004d2cce85ae7a598fb;hp=0596ef8bd783a026eb9db7b6b396f771a40325ed;hpb=0083c3fb2c66812489631c7551709d121de51625;p=clitk.git diff --git a/tools/clitkCorrelationRatioImageToImageMetric.txx b/tools/clitkCorrelationRatioImageToImageMetric.txx index 0596ef8..2b45c2f 100644 --- a/tools/clitkCorrelationRatioImageToImageMetric.txx +++ b/tools/clitkCorrelationRatioImageToImageMetric.txx @@ -1,3 +1,20 @@ +/*========================================================================= + 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 @@ -5,10 +22,10 @@ * @file clitkCorrelationRatioImageToImageMetric.txx * @author Jef Vandemeulebroucke * @date July 30 18:14:53 2007 - * + * * @brief Compute the correlation ratio between 2 images - * - * + * + * */ #include "clitkCorrelationRatioImageToImageMetric.h" @@ -22,7 +39,7 @@ namespace clitk /* * Constructor */ -template +template CorrelationRatioImageToImageMetric ::CorrelationRatioImageToImageMetric() { @@ -30,39 +47,36 @@ CorrelationRatioImageToImageMetric } -template +template void CorrelationRatioImageToImageMetric ::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::max(); double fixedImageMax = NumericTraits::NonpositiveMin(); typedef ImageRegionConstIterator FixedIteratorType; - FixedIteratorType fixedImageIterator( + FixedIteratorType fixedImageIterator( this->m_FixedImage, this->GetFixedImageRegion() ); - for ( fixedImageIterator.GoToBegin(); - !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) - { + for ( fixedImageIterator.GoToBegin(); + !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) { double sample = static_cast( fixedImageIterator.Get() ); - if ( sample < fixedImageMin ) - { + if ( sample < fixedImageMin ) { fixedImageMin = sample; - } + } - if ( sample > fixedImageMax ) - { + if ( sample > fixedImageMax ) { fixedImageMax = sample; - } } + } // Compute binsize for the fixedImage m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / m_NumberOfBins; @@ -77,7 +91,7 @@ CorrelationRatioImageToImageMetric /* * Get the match Measure */ -template +template typename CorrelationRatioImageToImageMetric::MeasureType CorrelationRatioImageToImageMetric ::GetValue( const TransformParametersType & parameters ) const @@ -87,10 +101,9 @@ CorrelationRatioImageToImageMetric FixedImageConstPointer fixedImage = this->m_FixedImage; - if( !fixedImage ) - { + if( !fixedImage ) { itkExceptionMacro( << "Fixed image has not been assigned" ); - } + } typedef itk::ImageRegionConstIteratorWithIndex FixedIteratorType; @@ -105,86 +118,79 @@ CorrelationRatioImageToImageMetric this->SetTransformParameters( parameters ); - //temporary measures for the calculation + //temporary measures for the calculation RealType mSMV=0; RealType mMSV=0; - while(!ti.IsAtEnd()) - { + 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 ) ) - { + 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 ) ) - { + 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( 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_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( 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++; } - 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])); - } - } + ++ti; + } - //Normalize with the global measures - measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted)); - return measure; + 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; + + } } @@ -194,7 +200,7 @@ CorrelationRatioImageToImageMetric /* * Get the Derivative Measure */ -template < class TFixedImage, class TMovingImage> +template < class TFixedImage, class TMovingImage> void CorrelationRatioImageToImageMetric ::GetDerivative( const TransformParametersType & parameters, @@ -202,27 +208,25 @@ CorrelationRatioImageToImageMetric { itkDebugMacro("GetDerivative( " << parameters << " ) "); - - if( !this->GetGradientImage() ) - { + + if( !this->GetGradientImage() ) { itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()"); - } + } FixedImageConstPointer fixedImage = this->m_FixedImage; - if( !fixedImage ) - { + if( !fixedImage ) { itkExceptionMacro( << "Fixed image has not been assigned" ); - } + } const unsigned int ImageDimension = FixedImageType::ImageDimension; typedef itk::ImageRegionConstIteratorWithIndex< - FixedImageType> FixedIteratorType; + FixedImageType> FixedIteratorType; typedef itk::ImageRegionConstIteratorWithIndex< - ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType; + ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType; FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); @@ -239,119 +243,106 @@ CorrelationRatioImageToImageMetric ti.GoToBegin(); - while(!ti.IsAtEnd()) - { + while(!ti.IsAtEnd()) { index = ti.GetIndex(); - + typename Superclass::InputPointType inputPoint; fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( 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 ) ) - { + if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) { ++ti; continue; - } + } - if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) - { + if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) { const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint ); const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); + this->m_Transform->GetJacobian( inputPoint ); + - const RealType fixedValue = ti.Value(); this->m_NumberOfPixelsCounted++; - const RealType diff = movingValue - fixedValue; + const RealType diff = movingValue - fixedValue; - // Get the gradient by NearestNeighboorInterpolation: + // 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 - MovingImageContinuousIndexType; + MovingImageContinuousIndexType; MovingImageContinuousIndexType tempIndex; this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex ); - typename MovingImageType::IndexType mappedIndex; - for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) - { + typename MovingImageType::IndexType mappedIndex; + for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) { mappedIndex[j] = static_cast( vnl_math_rnd( tempIndex[j] ) ); - } + } - const GradientPixelType gradient = + const GradientPixelType gradient = this->GetGradientImage()->GetPixel( mappedIndex ); - for(unsigned int par=0; par::Zero; - for(unsigned int dim=0; dimm_NumberOfPixelsCounted ) - { + if( !this->m_NumberOfPixelsCounted ) { itkExceptionMacro(<<"All the points mapped to outside of the moving image"); - } - else - { - for(unsigned int i=0; im_NumberOfPixelsCounted; - } } + } } /* - * Get both the match Measure and the Derivative Measure + * Get both the match Measure and the Derivative Measure */ -template +template void CorrelationRatioImageToImageMetric -::GetValueAndDerivative(const TransformParametersType & parameters, +::GetValueAndDerivative(const TransformParametersType & parameters, MeasureType & value, DerivativeType & derivative) const { itkDebugMacro("GetValueAndDerivative( " << parameters << " ) "); - if( !this->GetGradientImage() ) - { + if( !this->GetGradientImage() ) { itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()"); - } + } FixedImageConstPointer fixedImage = this->m_FixedImage; - if( !fixedImage ) - { + if( !fixedImage ) { itkExceptionMacro( << "Fixed image has not been assigned" ); - } + } const unsigned int ImageDimension = FixedImageType::ImageDimension; typedef itk::ImageRegionConstIteratorWithIndex< - FixedImageType> FixedIteratorType; + FixedImageType> FixedIteratorType; typedef itk::ImageRegionConstIteratorWithIndex< - ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType; + ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType; FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); @@ -370,88 +361,77 @@ CorrelationRatioImageToImageMetric ti.GoToBegin(); - while(!ti.IsAtEnd()) - { + while(!ti.IsAtEnd()) { index = ti.GetIndex(); - + typename Superclass::InputPointType inputPoint; fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( 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 ) ) - { + if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) { ++ti; continue; - } + } - if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) - { + if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) { const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint ); const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); + this->m_Transform->GetJacobian( inputPoint ); + - const RealType fixedValue = ti.Value(); this->m_NumberOfPixelsCounted++; - const RealType diff = movingValue - fixedValue; - + const RealType diff = movingValue - fixedValue; + measure += diff * diff; - // Get the gradient by NearestNeighboorInterpolation: + // 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 - MovingImageContinuousIndexType; + MovingImageContinuousIndexType; MovingImageContinuousIndexType tempIndex; this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex ); - typename MovingImageType::IndexType mappedIndex; - for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) - { + typename MovingImageType::IndexType mappedIndex; + for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) { mappedIndex[j] = static_cast( vnl_math_rnd( tempIndex[j] ) ); - } + } - const GradientPixelType gradient = + const GradientPixelType gradient = this->GetGradientImage()->GetPixel( mappedIndex ); - for(unsigned int par=0; par::Zero; - for(unsigned int dim=0; dimm_NumberOfPixelsCounted ) - { + if( !this->m_NumberOfPixelsCounted ) { itkExceptionMacro(<<"All the points mapped to outside of the moving image"); - } - else - { - for(unsigned int i=0; im_NumberOfPixelsCounted; - } - measure /= this->m_NumberOfPixelsCounted; } + measure /= this->m_NumberOfPixelsCounted; + } value = measure;