From: tbaudier Date: Thu, 18 Feb 2016 14:24:54 +0000 (+0100) Subject: Remove ITK3 and ITK4.2 tests and dependencies X-Git-Tag: v1.4.0~58 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=aa66e8e3252e2d42c7a1fd8b21336a5a0a65835d;p=clitk.git Remove ITK3 and ITK4.2 tests and dependencies Now, minimum required version is ITK4.3 --- diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake index 62a6a15..a871fec 100644 --- a/cmake/dependencies.cmake +++ b/cmake/dependencies.cmake @@ -59,15 +59,10 @@ endif() #========================================================= ### Check if ITK was compiled with SYSTEM_GDCM = ON set(CLITK_USE_SYSTEM_GDCM FALSE) -if(ITK_VERSION_MAJOR LESS "4") - if(ITK_USE_SYSTEM_GDCM) - set(CLITK_USE_SYSTEM_GDCM TRUE) - endif(ITK_USE_SYSTEM_GDCM) -else() - # ITK4 creates a target for each gdcm library when it compiles GDCM - get_target_property(GDCMDICTTARG gdcmDICT TYPE ) - if(NOT GDCMDICTTARG) - set(CLITK_USE_SYSTEM_GDCM TRUE) - endif() +# ITK4 creates a target for each gdcm library when it compiles GDCM +get_target_property(GDCMDICTTARG gdcmDICT TYPE ) +if(NOT GDCMDICTTARG) + set(CLITK_USE_SYSTEM_GDCM TRUE) endif() + diff --git a/common/clitkIO.cxx b/common/clitkIO.cxx index aa42db8..63e3be1 100644 --- a/common/clitkIO.cxx +++ b/common/clitkIO.cxx @@ -43,16 +43,13 @@ #include "clitkUSVoxImageIOFactory.h" #include "clitkSvlImageIOFactory.h" #endif -#if ITK_VERSION_MAJOR >= 4 - #include "itkGDCMImageIOFactory.h" - #include "itkPNGImageIOFactory.h" -#endif +#include "itkGDCMImageIOFactory.h" +#include "itkPNGImageIOFactory.h" //-------------------------------------------------------------------- // Register factories void clitk::RegisterClitkFactories() { -#if ITK_VERSION_MAJOR >= 4 std::list< itk::ObjectFactoryBase * > fl = itk::GDCMImageIOFactory::GetRegisteredFactories(); for (std::list< itk::ObjectFactoryBase * >::iterator it = fl.begin(); it != fl.end(); ++it) if (dynamic_cast(*it)) @@ -68,7 +65,6 @@ void clitk::RegisterClitkFactories() itk::PNGImageIOFactory::UnRegisterFactory(*it); break; } -#endif #if CLITK_PRIVATE_FEATURES clitk::UsfImageIOFactory::RegisterOneFactory(); clitk::USVoxImageIOFactory::RegisterOneFactory(); @@ -76,9 +72,6 @@ void clitk::RegisterClitkFactories() #endif clitk::GateAsciiImageIOFactory::RegisterOneFactory(); clitk::DicomRTDoseIOFactory::RegisterOneFactory(); -#if ITK_VERSION_MAJOR <= 3 - itk::ImageIOFactory::RegisterBuiltInFactories(); -#endif clitk::VoxImageIOFactory::RegisterOneFactory(); clitk::VfImageIOFactory::RegisterOneFactory(); clitk::XdrImageIOFactory::RegisterOneFactory(); @@ -88,9 +81,7 @@ void clitk::RegisterClitkFactories() rtk::ImagXImageIOFactory::RegisterOneFactory(); rtk::XRadImageIOFactory::RegisterOneFactory(); clitk::EsrfHstImageIOFactory::RegisterOneFactory(); -#if ITK_VERSION_MAJOR >= 4 itk::GDCMImageIOFactory::RegisterOneFactory(); itk::PNGImageIOFactory::RegisterOneFactory(); -#endif } //// diff --git a/common/vvFromITK.h b/common/vvFromITK.h index 69e41a6..a51cfe2 100644 --- a/common/vvFromITK.h +++ b/common/vvFromITK.h @@ -51,9 +51,7 @@ static inline void ReadTimeSequence (vvImage::Pointer& vv_image, typename itk::I extractedRegion.SetIndex(start); typename FilterType::Pointer filter = FilterType::New(); -#if ITK_VERSION_MAJOR == 4 filter->SetDirectionCollapseToSubmatrix(); -#endif filter->SetExtractionRegion(extractedRegion); filter->SetInput(input); filter->ReleaseDataFlagOn(); diff --git a/common/vvImageReader.txx b/common/vvImageReader.txx index a797bf4..551e426 100644 --- a/common/vvImageReader.txx +++ b/common/vvImageReader.txx @@ -121,9 +121,7 @@ void vvImageReader::UpdateWithDimAndInputPixelType() filter->SetExtractionRegion(extractedRegion); filter->SetInput(reader->GetOutput()); filter->ReleaseDataFlagOn(); -#if ITK_VERSION_MAJOR == 4 filter->SetDirectionCollapseToSubmatrix(); -#endif try { mImage->AddItkImage(filter->GetOutput()); mImage->ComputeScalarRangeBase(filter->GetOutput()); diff --git a/itk/RelativePositionPropImageFilter.h b/itk/RelativePositionPropImageFilter.h index 11797e7..183f260 100644 --- a/itk/RelativePositionPropImageFilter.h +++ b/itk/RelativePositionPropImageFilter.h @@ -81,17 +81,10 @@ namespace itk * This filter is implemented using the propagation algorithm */ -#if ITK_VERSION_MAJOR == 4 template > -#else - template > -#endif class ITK_EXPORT RelativePositionPropImageFilter : public ImageToImageFilter< TInputImage, TOutputImage > { diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx index a3759f4..c97ba0e 100644 --- a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx +++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx @@ -33,11 +33,7 @@ #include #include #include -#if ITK_VERSION_MAJOR >= 4 - #include -#else - #include -#endif +#include // itk [Bloch et al] #include "RelativePositionPropImageFilter.h" @@ -414,15 +410,9 @@ GenerateData() // Divide by the number of relpos if (GetNumberOfAngles() != 1) { -#if ITK_VERSION_MAJOR >= 4 typedef itk::DivideImageFilter DivideFilter; typename DivideFilter::Pointer divideFilter = DivideFilter::New(); divideFilter->SetConstant2(GetNumberOfAngles()); -#else - typedef itk::DivideByConstantImageFilter DivideFilter; - typename DivideFilter::Pointer divideFilter = DivideFilter::New(); - divideFilter->SetConstant(GetNumberOfAngles()); -#endif divideFilter->SetInput(m_FuzzyMap); divideFilter->Update(); m_FuzzyMap = divideFilter->GetOutput(); diff --git a/itk/clitkBackProjectImageFilter.h b/itk/clitkBackProjectImageFilter.h index 689e678..8ff8def 100644 --- a/itk/clitkBackProjectImageFilter.h +++ b/itk/clitkBackProjectImageFilter.h @@ -242,12 +242,7 @@ namespace clitk void BeforeThreadedGenerateData(void ); // Threaded Generate Data -#if ITK_VERSION_MAJOR >= 4 void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId ); -#else - void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, int threadId ); -#endif - //------------------------------------------------ //Member data diff --git a/itk/clitkBackProjectImageFilter.txx b/itk/clitkBackProjectImageFilter.txx index b1f84e2..f981291 100644 --- a/itk/clitkBackProjectImageFilter.txx +++ b/itk/clitkBackProjectImageFilter.txx @@ -302,12 +302,7 @@ namespace clitk //----------------------------------------------------------------------- template void - BackProjectImageFilter -#if ITK_VERSION_MAJOR >= 4 - ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId ) -#else - ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, int threadId ) -#endif + BackProjectImageFilter::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId ) { //Projection pointer InputImageConstPointer inputPtr=this->GetInput(); diff --git a/itk/clitkComposeVFFilter.h b/itk/clitkComposeVFFilter.h index cf1ada2..d043c1a 100644 --- a/itk/clitkComposeVFFilter.h +++ b/itk/clitkComposeVFFilter.h @@ -75,11 +75,7 @@ namespace clitk //======================================================================================== //Threaded execution should implement generate threaded data -#if ITK_VERSION_MAJOR >= 4 - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); -#else - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); -#endif + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); bool m_Verbose; PixelType m_EdgePaddingValue; diff --git a/itk/clitkComposeVFFilter.txx b/itk/clitkComposeVFFilter.txx index a50ca26..5cc5db1 100644 --- a/itk/clitkComposeVFFilter.txx +++ b/itk/clitkComposeVFFilter.txx @@ -35,13 +35,8 @@ namespace clitk //========================================================================================================================= //update the output for the outputRegionForThread -#if ITK_VERSION_MAJOR >= 4 template void ComposeVFFilter::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ) -#else - template - void ComposeVFFilter::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ) -#endif { //Get pointer to the output diff --git a/itk/clitkExtractImageFilter.h b/itk/clitkExtractImageFilter.h index ce31fd9..98a40a6 100644 --- a/itk/clitkExtractImageFilter.h +++ b/itk/clitkExtractImageFilter.h @@ -94,11 +94,7 @@ protected: const OutputImageRegionType &srcRegion); -#if ITK_VERSION_MAJOR >= 4 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); -#else - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); -#endif InputImageRegionType m_ExtractionRegion; OutputImageRegionType m_OutputImageRegion; diff --git a/itk/clitkExtractImageFilter.txx b/itk/clitkExtractImageFilter.txx index 90fee3b..5872469 100644 --- a/itk/clitkExtractImageFilter.txx +++ b/itk/clitkExtractImageFilter.txx @@ -235,12 +235,7 @@ ExtractImageFilter */ template void -ExtractImageFilter -#if ITK_VERSION_MAJOR >= 4 -::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId) -#else -::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId) -#endif +ExtractImageFilter::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId) { itkDebugMacro(<<"Actually executing"); diff --git a/itk/clitkExtractSliceFilter.txx b/itk/clitkExtractSliceFilter.txx index e983be7..af008ec 100644 --- a/itk/clitkExtractSliceFilter.txx +++ b/itk/clitkExtractSliceFilter.txx @@ -105,11 +105,7 @@ GenerateData() { m_size[GetDirection()] = 0; m_region.SetSize(m_size); int start = m_index[GetDirection()]; -#if ITK_VERSION_MAJOR >= 4 this->SetNumberOfIndexedInputs(m_NumberOfSlices); -#else - this->SetNumberOfOutputs(m_NumberOfSlices); -#endif //-------------------------------------------------------------------- // loop ExtractImageFilter with region updated, push_back @@ -122,9 +118,7 @@ GenerateData() { m_index[GetDirection()] = start + i; m_region.SetIndex(m_index); extract->SetExtractionRegion(m_region); -#if ITK_VERSION_MAJOR == 4 extract->SetDirectionCollapseToSubmatrix(); -#endif extract->Update(); this->SetNthOutput(i, extract->GetOutput()); } diff --git a/itk/clitkInvertVFFilter.txx b/itk/clitkInvertVFFilter.txx index 34500a3..525b486 100644 --- a/itk/clitkInvertVFFilter.txx +++ b/itk/clitkInvertVFFilter.txx @@ -75,11 +75,7 @@ protected: //the actual processing void BeforeThreadedGenerateData(); -#if ITK_VERSION_MAJOR >= 4 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); -#else - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); -#endif //member data typename WeightsImageType::Pointer m_Weights; @@ -117,11 +113,7 @@ void HelperClass1::BeforeThreadedGenerateData() //========================================================================================================================= //update the output for the outputRegionForThread template -#if ITK_VERSION_MAJOR >= 4 void HelperClass1::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ) -#else -void HelperClass1::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ) -#endif { // std::cout << "HelperClass1::ThreadedGenerateData - IN " << threadId << std::endl; //Get pointer to the input @@ -297,11 +289,7 @@ protected: //the actual processing -#if ITK_VERSION_MAJOR >= 4 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); -#else - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); -#endif //member data typename WeightsImageType::Pointer m_Weights; @@ -326,11 +314,7 @@ template HelperClass2= 4 template void HelperClass2::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ) -#else -template void HelperClass2::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ) -#endif { // std::cout << "HelperClass2::ThreadedGenerateData - IN " << threadId << std::endl; diff --git a/itk/clitkReconstructThroughDilationImageFilter.h b/itk/clitkReconstructThroughDilationImageFilter.h index 32d31fe..4c7f5eb 100644 --- a/itk/clitkReconstructThroughDilationImageFilter.h +++ b/itk/clitkReconstructThroughDilationImageFilter.h @@ -39,11 +39,7 @@ #include "itkConnectedComponentImageFilter.h" #include "itkStatisticsImageFilter.h" #include "itkCastImageFilter.h" -#if ITK_VERSION_MAJOR >= 4 - #include "itkTestingComparisonImageFilter.h" -#else - #include "itkDifferenceImageFilter.h" -#endif +#include "itkTestingComparisonImageFilter.h" #include "itkThresholdImageFilter.h" namespace clitk diff --git a/itk/clitkReconstructThroughDilationImageFilter.txx b/itk/clitkReconstructThroughDilationImageFilter.txx index 8e30eb3..de290f9 100644 --- a/itk/clitkReconstructThroughDilationImageFilter.txx +++ b/itk/clitkReconstructThroughDilationImageFilter.txx @@ -68,11 +68,7 @@ namespace clitk typedef itk::StatisticsImageFilter StatisticsImageFilterType; typedef itk::BinaryBallStructuringElement KernelType; typedef clitk::ConditionalBinaryDilateImageFilter ConditionalBinaryDilateImageFilterType; -#if ITK_VERSION_MAJOR >= 4 typedef itk::Testing::ComparisonImageFilter DifferenceImageFilterType; -#else - typedef itk::DifferenceImageFilter DifferenceImageFilterType; -#endif typedef itk::CastImageFilter OutputCastImageFilterType; typedef clitk::SetBackgroundImageFilter SetBackgroundImageFilterType; diff --git a/itk/itkFlexibleBinaryFunctorImageFilter.h b/itk/itkFlexibleBinaryFunctorImageFilter.h index 18062a2..6b85a09 100644 --- a/itk/itkFlexibleBinaryFunctorImageFilter.h +++ b/itk/itkFlexibleBinaryFunctorImageFilter.h @@ -145,13 +145,8 @@ protected: * * \sa ImageToImageFilter::ThreadedGenerateData(), * ImageToImageFilter::GenerateData() */ -#if ITK_VERSION_MAJOR >= 4 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); -#else - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, - int threadId ); -#endif private: FlexibleBinaryFunctorImageFilter(const Self&); //purposely not implemented diff --git a/itk/itkFlexibleBinaryFunctorImageFilter.txx b/itk/itkFlexibleBinaryFunctorImageFilter.txx index d8258cc..d7c80e2 100644 --- a/itk/itkFlexibleBinaryFunctorImageFilter.txx +++ b/itk/itkFlexibleBinaryFunctorImageFilter.txx @@ -97,12 +97,7 @@ FlexibleBinaryFunctorImageFilter void FlexibleBinaryFunctorImageFilter -::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, -#if ITK_VERSION_MAJOR >= 4 - itk::ThreadIdType threadId ) -#else - int threadId) -#endif +::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId ) { const unsigned int dim = Input1ImageType::ImageDimension; diff --git a/registration/clitkAffineRegistrationGenericFilter.cxx b/registration/clitkAffineRegistrationGenericFilter.cxx index 6b0ecff..bcf5fd4 100644 --- a/registration/clitkAffineRegistrationGenericFilter.cxx +++ b/registration/clitkAffineRegistrationGenericFilter.cxx @@ -175,11 +175,8 @@ void AffineRegistrationGenericFilter::UpdateWithInputImageType() typedef typename InputImageType::PixelType PixelType; //typedef typename InputImageType::ImageDimension Dimension; - -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 bool threadsGiven=m_ArgsInfo.threads_given; int threads=m_ArgsInfo.threads_arg; -#endif //Coordinate Representation typedef double TCoordRep; @@ -396,11 +393,7 @@ void AffineRegistrationGenericFilter::UpdateWithInputImageType() typename MetricType::Pointer metric=genericMetric->GetMetricPointer(); if (movingMask) metric->SetMovingImageMask(movingMask); -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 if (threadsGiven) metric->SetNumberOfThreads( threads ); -#else - if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<StartRegistration(); -#else registration->Update(); -#endif } catch ( itk::ExceptionObject & err ) { std::cerr << "ExceptionObject caught !" << std::endl; std::cerr << err << std::endl; diff --git a/registration/clitkBLUTDIRGenericFilter.cxx b/registration/clitkBLUTDIRGenericFilter.cxx index b49d990..9f7b84f 100644 --- a/registration/clitkBLUTDIRGenericFilter.cxx +++ b/registration/clitkBLUTDIRGenericFilter.cxx @@ -396,9 +396,7 @@ namespace clitk // Crop the fixedImage to the bounding box to facilitate multi-resolution typedef itk::ExtractImageFilter ExtractImageFilterType; typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New(); -#if ITK_VERSION_MAJOR == 4 extractImageFilter->SetDirectionCollapseToSubmatrix(); -#endif extractImageFilter->SetInput(fixedImage); extractImageFilter->SetExtractionRegion(transformRegion); extractImageFilter->Update(); @@ -652,16 +650,10 @@ namespace clitk typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType; typename MetricType::Pointer metric=genericMetric->GetMetricPointer(); if (movingMask) metric->SetMovingImageMask(movingMask); - -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 if (threadsGiven) { metric->SetNumberOfThreads( threads ); if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl; } -#else - if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<StartRegistration(); -#else registration->Update(); -#endif } catch( itk::ExceptionObject & err ) { @@ -807,8 +795,6 @@ namespace clitk # else typedef itk::TransformToDisplacementFieldFilter ConvertorType; # endif -#else - typedef itk::TransformToDeformationFieldSource ConvertorType; #endif typename ConvertorType::Pointer filter= ConvertorType::New(); filter->SetNumberOfThreads(1); diff --git a/registration/clitkBSplineDeformableTransform.h b/registration/clitkBSplineDeformableTransform.h index 392247f..8f5a278 100644 --- a/registration/clitkBSplineDeformableTransform.h +++ b/registration/clitkBSplineDeformableTransform.h @@ -65,9 +65,7 @@ namespace clitk /** Standard parameters container. */ typedef typename Superclass::ParametersType ParametersType; -#if ITK_VERSION_MAJOR >= 4 typedef typename Superclass::NumberOfParametersType NumberOfParametersType; -#endif /** Standard Jacobian container. */ typedef typename Superclass::JacobianType JacobianType; @@ -257,22 +255,14 @@ namespace clitk } /** Compute the Jacobian Matrix of the transformation at one point */ -#if ITK_VERSION_MAJOR >= 4 virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const; virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const { itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() ); } -#else - virtual const JacobianType& GetJacobian(const InputPointType &point ) const; -#endif /** Return the number of parameters that completely define the Transfom */ -#if ITK_VERSION_MAJOR >= 4 virtual NumberOfParametersType GetNumberOfParameters(void) const; -#else - virtual unsigned int GetNumberOfParameters(void) const; -#endif /** Return the number of parameters per dimension */ unsigned int GetNumberOfParametersPerDimension(void) const; @@ -378,10 +368,7 @@ namespace clitk // VD Add MultipleBSplineDeformableTransform as friend to facilitate wrapping friend class MultipleBSplineDeformableTransform; -#if ITK_VERSION_MAJOR >= 4 mutable JacobianType m_SharedDataBSplineJacobian; -#endif - }; //class BSplineDeformableTransform diff --git a/registration/clitkBSplineDeformableTransform.txx b/registration/clitkBSplineDeformableTransform.txx index 4d0f50d..60a2136 100644 --- a/registration/clitkBSplineDeformableTransform.txx +++ b/registration/clitkBSplineDeformableTransform.txx @@ -31,12 +31,7 @@ namespace clitk // Constructor with default arguments template - BSplineDeformableTransform -#if ITK_VERSION_MAJOR >= 4 - ::BSplineDeformableTransform():Superclass(0) -#else - ::BSplineDeformableTransform():Superclass(OutputDimension,0) -#endif + BSplineDeformableTransform::BSplineDeformableTransform():Superclass(0) { unsigned int i; @@ -255,11 +250,7 @@ namespace clitk // Get the number of parameters template -#if ITK_VERSION_MAJOR >= 4 typename BSplineDeformableTransform::NumberOfParametersType -#else - unsigned int -#endif BSplineDeformableTransform ::GetNumberOfParameters(void) const { @@ -584,18 +575,11 @@ namespace clitk m_CoefficientImage = m_WrappedImage; //JV Wrap jacobian into OutputDimension X Vectorial images -#if ITK_VERSION_MAJOR >= 4 this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); -#else - this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); -#endif // Use memset to set the memory -#if ITK_VERSION_MAJOR >= 4 JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_SharedDataBSplineJacobian.data_block()); -#else - JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_Jacobian.data_block()); -#endif + memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType)); m_LastJacobianIndex = m_ValidRegion.GetIndex(); m_NeedResetJacobian = false; @@ -920,7 +904,6 @@ namespace clitk // JV weights are identical as for transformpoint, could be done simultaneously in metric!!!! // Compute the Jacobian in one position -#if ITK_VERSION_MAJOR >= 4 template void BSplineDeformableTransform @@ -982,76 +965,7 @@ namespace clitk jacobian = m_SharedDataBSplineJacobian; } -#else - template - const - typename BSplineDeformableTransform - ::JacobianType & - BSplineDeformableTransform - ::GetJacobian( const InputPointType & point ) const - { - // Can only compute Jacobian if parameters are set via - // SetParameters or SetParametersByValue - // if( m_InputParametersPointer == NULL ) - // { - // itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" ); - // } - - if (m_NeedResetJacobian) - ResetJacobian(); - //======================================================== - // For each dimension, copy the weight to the support region - //======================================================== - - // Check if inside mask - if(m_Mask && !(m_Mask->IsInside(point) ) ) - { - // Outside: no (deformable) displacement - return this->m_Jacobian; - } - - //Get index - this->TransformPointToContinuousIndex( point, m_Index ); - - // NOTE: if the support region does not lie totally within the grid - // we assume zero displacement and return the input point - if ( !this->InsideValidRegion( m_Index ) ) - { - return this->m_Jacobian; - } - - //Compute interpolation weights - const WeightsDataType *weights=NULL; - m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex); - m_SupportRegion.SetIndex( m_LastJacobianIndex ); - - //Reset the iterators - unsigned int j = 0; - for ( j = 0; j < OutputDimension; j++ ) - m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion); - - // For each dimension, copy the weight to the support region - while ( ! (m_Iterator[0]).IsAtEnd() ) - { - //copy weight to jacobian image - for ( j = 0; j < OutputDimension; j++ ) - { - m_ZeroVector[j]=*weights; - (m_Iterator[j]).Set( m_ZeroVector); - m_ZeroVector[j]=itk::NumericTraits::Zero; - ++(m_Iterator[j]); - } - // go to next coefficient in the support region - weights++; - } - m_NeedResetJacobian = true; - - // Return the result - return this->m_Jacobian; - - } -#endif template diff --git a/registration/clitkConvertBLUTCoeffsToVFFilter.h b/registration/clitkConvertBLUTCoeffsToVFFilter.h index 76b1578..40eb129 100644 --- a/registration/clitkConvertBLUTCoeffsToVFFilter.h +++ b/registration/clitkConvertBLUTCoeffsToVFFilter.h @@ -13,8 +13,6 @@ # else # include "itkTransformToDisplacementFieldFilter.h" # endif -#else -# include "itkTransformToDeformationFieldSource.h" #endif namespace clitk @@ -55,8 +53,6 @@ namespace clitk # else typedef itk::TransformToDisplacementFieldFilter ConvertorType; # endif -#else - typedef itk::TransformToDeformationFieldSource ConvertorType; #endif itkNewMacro(Self); diff --git a/registration/clitkConvertBLUTCoeffsToVFFilter.txx b/registration/clitkConvertBLUTCoeffsToVFFilter.txx index d10bcbb..ef63900 100644 --- a/registration/clitkConvertBLUTCoeffsToVFFilter.txx +++ b/registration/clitkConvertBLUTCoeffsToVFFilter.txx @@ -143,12 +143,7 @@ namespace clitk typedef clitk::VectorImageToImageFilter FilterType; typename FilterType::Pointer component_filter[BLUTCoefficientImageType::ImageDimension]; - -#if ITK_VERSION_MAJOR >= 4 typename ITKTransformType::CoefficientImageArray coefficient_images; -#else - typename ITKTransformType::ImagePointer coefficient_images[BLUTCoefficientImageType::ImageDimension]; -#endif for (unsigned int i=0; i < BLUTCoefficientImageType::ImageDimension; i++) { component_filter[i] = FilterType::New(); @@ -158,7 +153,6 @@ namespace clitk coefficient_images[i] = component_filter[i]->GetOutput(); } -#if ITK_VERSION_MAJOR >= 4 // RP: 16/01/2013 // ATTENTION: Apparently, there's a bug in the SetCoefficientImages function of ITK 4.x // I needed to use SetParametersByValue instead. @@ -174,9 +168,6 @@ namespace clitk m_ITKTransform->SetGridRegion(input->GetLargestPossibleRegion()); m_ITKTransform->SetGridSpacing(input->GetSpacing()); m_ITKTransform->SetParametersByValue(params); -#else - m_ITKTransform->SetCoefficientImage(coefficient_images); -#endif m_GenericTransform = m_ITKTransform; } diff --git a/registration/clitkCorrelationRatioImageToImageMetric.txx b/registration/clitkCorrelationRatioImageToImageMetric.txx index a400e51..b52c077 100644 --- a/registration/clitkCorrelationRatioImageToImageMetric.txx +++ b/registration/clitkCorrelationRatioImageToImageMetric.txx @@ -266,14 +266,8 @@ CorrelationRatioImageToImageMetric if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) { const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint ); -#if ITK_VERSION_MAJOR >= 4 TransformJacobianType jacobian; this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint , jacobian); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif - const RealType fixedValue = ti.Value(); this->m_NumberOfPixelsCounted++; @@ -389,14 +383,8 @@ CorrelationRatioImageToImageMetric if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) { const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint ); -#if ITK_VERSION_MAJOR >= 4 TransformJacobianType jacobian; this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian ); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif - const RealType fixedValue = ti.Value(); this->m_NumberOfPixelsCounted++; diff --git a/registration/clitkDeformationFieldTransform.h b/registration/clitkDeformationFieldTransform.h index 06c5288..a997d86 100644 --- a/registration/clitkDeformationFieldTransform.h +++ b/registration/clitkDeformationFieldTransform.h @@ -109,7 +109,6 @@ namespace clitk return OutputCovariantVectorType(); } -#if ITK_VERSION_MAJOR >= 4 virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const { itkExceptionMacro( << "DeformationFieldTransform doesn't declare ComputeJacobianWithRespectToParameters" ); @@ -118,13 +117,6 @@ namespace clitk { itkExceptionMacro( << "DeformationFieldTransform doesn't declare ComputeJacobianWithRespectToPosition" ); } -#else - virtual const JacobianType& GetJacobian(const InputPointType &point ) const - { - itkExceptionMacro( << "DeformationFieldTransform doesn't declare GetJacobian" ); - return this->m_Jacobian; - } -#endif protected: DeformationFieldTransform(); diff --git a/registration/clitkDeformationFieldTransform.txx b/registration/clitkDeformationFieldTransform.txx index 049f270..c93b456 100644 --- a/registration/clitkDeformationFieldTransform.txx +++ b/registration/clitkDeformationFieldTransform.txx @@ -25,12 +25,8 @@ namespace clitk // Constructor template - DeformationFieldTransform -#if ITK_VERSION_MAJOR >= 4 - ::DeformationFieldTransform():Superclass(1) -#else - ::DeformationFieldTransform():Superclass(OutputDimension,1) -#endif + DeformationFieldTransform::DeformationFieldTransform():Superclass(1) + { m_DeformationField=NULL; m_Interpolator=DefaultInterpolatorType::New(); diff --git a/registration/clitkDemonsDeformableRegistrationGenericFilter.txx b/registration/clitkDemonsDeformableRegistrationGenericFilter.txx index 6132dc5..d4b0cb2 100644 --- a/registration/clitkDemonsDeformableRegistrationGenericFilter.txx +++ b/registration/clitkDemonsDeformableRegistrationGenericFilter.txx @@ -151,11 +151,7 @@ namespace clitk //find the multiresolution filter // typedef typename RegistrationFilterType::FixedImageType InternalImageType; // typedef typename RegistrationFilterType::MovingImageType MovingImageType; -#if ITK_VERSION_MAJOR >= 4 typedef typename RegistrationFilterType::DisplacementFieldType DisplacementFieldType; -#else - typedef typename RegistrationFilterType::DeformationFieldType DisplacementFieldType; -#endif typedef clitk::MultiResolutionPDEDeformableRegistration MultiResolutionRegistrationType; typedef CommandResolutionLevelUpdate LevelObserver; @@ -537,11 +533,7 @@ namespace clitk //JV TODO // pdeFilter->SetMaximumError(m_ArgsInfo.maxError_arg); // pdeFilter->SetMaximumKernelWidth(m_ArgsInfo.maxError_arg); -#if ITK_VERSION_MAJOR >= 4 pdeFilter->SetSmoothDisplacementField(!m_ArgsInfo.fluid_flag); -#else - pdeFilter->SetSmoothDeformationField(!m_ArgsInfo.fluid_flag); -#endif pdeFilter->SetSmoothUpdateField(m_ArgsInfo.fluid_flag); pdeFilter->SetUseImageSpacing( m_ArgsInfo.spacing_flag ); @@ -607,11 +599,7 @@ namespace clitk typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType; typename WarpFilterType::Pointer warp = WarpFilterType::New(); -#if ITK_VERSION_MAJOR >= 4 warp->SetDisplacementField( deformationField ); -#else - warp->SetDeformationField( deformationField ); -#endif warp->SetInput( movingImageReader->GetOutput() ); warp->SetOutputOrigin( fixedImage->GetOrigin() ); warp->SetOutputSpacing( fixedImage->GetSpacing() ); diff --git a/registration/clitkGenericMetric.h b/registration/clitkGenericMetric.h index ad03c71..7b84008 100644 --- a/registration/clitkGenericMetric.h +++ b/registration/clitkGenericMetric.h @@ -143,11 +143,8 @@ private: typename FixedImageType::Pointer m_FixedImage; typename FixedImageMaskType::ConstPointer m_FixedImageMask; -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 FixedImagePixelType m_FixedImageSamplesIntensityThreshold; bool m_UseFixedImageSamplesIntensityThreshold; -#endif - }; } // end namespace clitk diff --git a/registration/clitkGenericMetric.txx b/registration/clitkGenericMetric.txx index cbbcaa0..852c19f 100644 --- a/registration/clitkGenericMetric.txx +++ b/registration/clitkGenericMetric.txx @@ -35,10 +35,8 @@ GenericMetric::GenericMetric() m_Maximize=false; m_Verbose=false; m_FixedImageRegionGiven=false; -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 m_FixedImageSamplesIntensityThreshold=0; m_UseFixedImageSamplesIntensityThreshold=false; -#endif m_FixedImageMask=NULL; } @@ -273,9 +271,6 @@ GenericMetric::GetMetricPointer( m_Metric->SetFixedImageRegion(m_FixedImageRegion); //m_Metric->SetFixedImageRegion(mask_region); - -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 - //============================================================================ // Set the lower intensity threshold if (m_ArgsInfo.intThreshold_given) { @@ -436,12 +431,6 @@ GenericMetric::GetMetricPointer( if (m_Verbose) std::cout<<"number of mask pixels "< 3 mutable std::ostringstream m_StopConditionDescription; -#else - mutable itk::OStringStream m_StopConditionDescription; -#endif BoundValueType m_LowerBound; BoundValueType m_UpperBound; BoundSelectionType m_BoundSelection; diff --git a/registration/clitkMatrixTransformToVFGenericFilter.h b/registration/clitkMatrixTransformToVFGenericFilter.h index cc6ae00..9aa5489 100644 --- a/registration/clitkMatrixTransformToVFGenericFilter.h +++ b/registration/clitkMatrixTransformToVFGenericFilter.h @@ -42,8 +42,6 @@ # else # include "itkTransformToDisplacementFieldFilter.h" # endif -#else -# include "itkTransformToDeformationFieldSource.h" #endif #include "itkAffineTransform.h" diff --git a/registration/clitkMatrixTransformToVFGenericFilter.txx b/registration/clitkMatrixTransformToVFGenericFilter.txx index a1241b5..8a36a88 100644 --- a/registration/clitkMatrixTransformToVFGenericFilter.txx +++ b/registration/clitkMatrixTransformToVFGenericFilter.txx @@ -84,8 +84,6 @@ namespace clitk # else typedef itk::TransformToDisplacementFieldFilter ConvertorType; # endif -#else - typedef itk::TransformToDeformationFieldSource ConvertorType; #endif typename ConvertorType::Pointer filter= ConvertorType::New(); diff --git a/registration/clitkMultiResolutionPDEDeformableRegistration.txx b/registration/clitkMultiResolutionPDEDeformableRegistration.txx index 5c6b0ae..43fd978 100644 --- a/registration/clitkMultiResolutionPDEDeformableRegistration.txx +++ b/registration/clitkMultiResolutionPDEDeformableRegistration.txx @@ -335,11 +335,7 @@ MultiResolutionPDEDeformableRegistration= 4 m_RegistrationFilter->SetInitialDisplacementField( NULL ); -#else - m_RegistrationFilter->SetInitialDeformationField( NULL ); -#endif } else { @@ -361,12 +357,7 @@ MultiResolutionPDEDeformableRegistrationGetOutput(); tempField->DisconnectPipeline(); -#if ITK_VERSION_MAJOR >= 4 m_RegistrationFilter->SetInitialDisplacementField( tempField ); -#else - m_RegistrationFilter->SetInitialDeformationField( tempField ); -#endif - } // setup registration filter and pyramids diff --git a/registration/clitkMultipleBSplineDeformableTransform.h b/registration/clitkMultipleBSplineDeformableTransform.h index 1b7909a..ae9f16f 100644 --- a/registration/clitkMultipleBSplineDeformableTransform.h +++ b/registration/clitkMultipleBSplineDeformableTransform.h @@ -59,9 +59,7 @@ namespace clitk /** Standard parameters container. */ typedef typename Superclass::ParametersType ParametersType; -#if ITK_VERSION_MAJOR >= 4 typedef typename Superclass::NumberOfParametersType NumberOfParametersType; -#endif /** Standard Jacobian container. */ typedef typename Superclass::JacobianType JacobianType; @@ -216,22 +214,14 @@ namespace clitk } /** Compute the Jacobian Matrix of the transformation at one point */ -#if ITK_VERSION_MAJOR >= 4 virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const; virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const { itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() ); } -#else - virtual const JacobianType& GetJacobian(const InputPointType &point ) const; -#endif /** Return the number of parameters that completely define the Transfom */ -#if ITK_VERSION_MAJOR >= 4 virtual NumberOfParametersType GetNumberOfParameters(void) const; -#else - virtual unsigned int GetNumberOfParameters(void) const; -#endif /** Return the number of parameters per dimension */ unsigned int GetNumberOfParametersPerDimension(void) const; @@ -278,9 +268,7 @@ namespace clitk std::vector m_parameters; mutable std::vector m_CoefficientImages; mutable int m_LastJacobian; -#if ITK_VERSION_MAJOR >= 4 mutable JacobianType m_SharedDataBSplineJacobian; -#endif void InitJacobian(); // FIXME it seems not used diff --git a/registration/clitkMultipleBSplineDeformableTransform.txx b/registration/clitkMultipleBSplineDeformableTransform.txx index 8df2706..ebcfae9 100644 --- a/registration/clitkMultipleBSplineDeformableTransform.txx +++ b/registration/clitkMultipleBSplineDeformableTransform.txx @@ -28,12 +28,7 @@ namespace clitk { // Constructor with default arguments template - MultipleBSplineDeformableTransform -#if ITK_VERSION_MAJOR >= 4 - ::MultipleBSplineDeformableTransform() : Superclass(0) -#else - ::MultipleBSplineDeformableTransform() : Superclass(OutputDimension, 0) -#endif + MultipleBSplineDeformableTransform::MultipleBSplineDeformableTransform() : Superclass(0) { m_nLabels = 1; m_labels = 0; @@ -329,11 +324,7 @@ namespace clitk #undef LOOP_ON_LABELS template -#if ITK_VERSION_MAJOR >= 4 inline typename MultipleBSplineDeformableTransform::NumberOfParametersType -#else - inline unsigned int -#endif MultipleBSplineDeformableTransform ::GetNumberOfParameters(void) const { @@ -433,7 +424,6 @@ namespace clitk return m_trans[lidx]->DeformablyTransformPoint(inputPoint); } -#if ITK_VERSION_MAJOR >= 4 template inline void MultipleBSplineDeformableTransform @@ -459,30 +449,6 @@ namespace clitk jacobian = this->m_SharedDataBSplineJacobian; } -#else - template - inline const typename MultipleBSplineDeformableTransform::JacobianType & - MultipleBSplineDeformableTransform - ::GetJacobian( const InputPointType & point ) const - { - if (m_LastJacobian != -1) - m_trans[m_LastJacobian]->ResetJacobian(); - - int lidx = 0; - if (m_labels) - lidx = m_labelInterpolator->Evaluate(point) - 1; - if (lidx == -1) - { - m_LastJacobian = lidx; - return this->m_Jacobian; - } - - m_trans[lidx]->GetJacobian(point); - m_LastJacobian = lidx; - - return this->m_Jacobian; - } -#endif template inline void @@ -497,13 +463,8 @@ namespace clitk MultipleBSplineDeformableTransform::InitJacobian() { unsigned numberOfParameters = this->GetNumberOfParameters(); -#if ITK_VERSION_MAJOR >= 4 this->m_SharedDataBSplineJacobian.set_size(OutputDimension, numberOfParameters); JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_SharedDataBSplineJacobian.data_block()); -#else - this->m_Jacobian.set_size(OutputDimension, numberOfParameters); - JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_Jacobian.data_block()); -#endif memset(jacobianDataPointer, 0, numberOfParameters * sizeof (JacobianPixelType)); unsigned tot = 0; diff --git a/registration/clitkNormalizedCorrelationImageToImageMetric.h b/registration/clitkNormalizedCorrelationImageToImageMetric.h index e6b704b..0dd921b 100644 --- a/registration/clitkNormalizedCorrelationImageToImageMetric.h +++ b/registration/clitkNormalizedCorrelationImageToImageMetric.h @@ -23,96 +23,7 @@ // This line can be removed once the optimized versions // gets integrated into the main directories. #include "itkConfigure.h" - -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "clitkOptNormalizedCorrelationImageToImageMetric.h" -#else - -#include "itkImageToImageMetric.h" -#include "itkCovariantVector.h" -#include "itkPoint.h" - - -namespace clitk -{ - -template < class TFixedImage, class TMovingImage > -class ITK_EXPORT NormalizedCorrelationImageToImageMetric : - public itk::ImageToImageMetric< TFixedImage, TMovingImage> -{ -public: - - /** Standard class typedefs. */ - typedef NormalizedCorrelationImageToImageMetric Self; - typedef itk::ImageToImageMetric Superclass; - - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(NormalizedCorrelationImageToImageMetric, itk::Object); - - - /** Types transferred from the base class */ - typedef typename Superclass::RealType RealType; - typedef typename Superclass::TransformType TransformType; - typedef typename Superclass::TransformPointer TransformPointer; - typedef typename Superclass::TransformParametersType TransformParametersType; - typedef typename Superclass::TransformJacobianType TransformJacobianType; - typedef typename Superclass::GradientPixelType GradientPixelType; - typedef typename Superclass::OutputPointType OutputPointType; - typedef typename Superclass::InputPointType InputPointType; - - typedef typename Superclass::MeasureType MeasureType; - typedef typename Superclass::DerivativeType DerivativeType; - typedef typename Superclass::FixedImageType FixedImageType; - typedef typename Superclass::MovingImageType MovingImageType; - typedef typename Superclass::FixedImageConstPointer FixedImageConstPointer; - typedef typename Superclass::MovingImageConstPointer MovingImageConstPointer; - - - /** Get the derivatives of the match measure. */ - void GetDerivative( const TransformParametersType & parameters, - DerivativeType & Derivative ) const; - - /** Get the value for single valued optimizers. */ - MeasureType GetValue( const TransformParametersType & parameters ) const; - - /** Get value and derivatives for multiple valued optimizers. */ - void GetValueAndDerivative( const TransformParametersType & parameters, - MeasureType& Value, DerivativeType& Derivative ) const; - - /** Set/Get SubtractMean boolean. If true, the sample mean is subtracted - * from the sample values in the cross-correlation formula and - * typically results in narrower valleys in the cost fucntion. - * Default value is false. */ - itkSetMacro( SubtractMean, bool ); - itkGetConstReferenceMacro( SubtractMean, bool ); - itkBooleanMacro( SubtractMean ); - -protected: - NormalizedCorrelationImageToImageMetric(); - virtual ~NormalizedCorrelationImageToImageMetric() {}; - void PrintSelf(std::ostream& os, itk::Indent indent) const; - -private: - NormalizedCorrelationImageToImageMetric(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented - - bool m_SubtractMean; - -}; - -} // end namespace clitk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "clitkNormalizedCorrelationImageToImageMetric.txx" -#endif - -#endif // opt #endif // _clitkNormalizedCorrelationImageToImageMetric.txx diff --git a/registration/clitkNormalizedCorrelationImageToImageMetric.txx b/registration/clitkNormalizedCorrelationImageToImageMetric.txx index 69aa81f..19286f6 100644 --- a/registration/clitkNormalizedCorrelationImageToImageMetric.txx +++ b/registration/clitkNormalizedCorrelationImageToImageMetric.txx @@ -24,492 +24,5 @@ // gets integrated into the main directories. #include "itkConfigure.h" -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "clitkOptNormalizedCorrelationImageToImageMetric.txx" -#else - - -#include "clitkNormalizedCorrelationImageToImageMetric.h" - -#include "itkImageRegionConstIteratorWithIndex.h" - - -namespace clitk -{ - -/* - * Constructor - */ -template -NormalizedCorrelationImageToImageMetric -::NormalizedCorrelationImageToImageMetric() -{ - m_SubtractMean = false; -} - -/* - * Get the match Measure - */ -template -typename NormalizedCorrelationImageToImageMetric::MeasureType -NormalizedCorrelationImageToImageMetric -::GetValue( const TransformParametersType & parameters ) const -{ - - FixedImageConstPointer fixedImage = this->m_FixedImage; - - if( !fixedImage ) { - itkExceptionMacro( << "Fixed image has not been assigned" ); - } - - typedef itk::ImageRegionConstIteratorWithIndex FixedIteratorType; - - FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); - - typename FixedImageType::IndexType index; - - MeasureType measure; - - this->m_NumberOfPixelsCounted = 0; - - this->SetTransformParameters( parameters ); - - typedef typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType; - - AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sf = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sm = itk::NumericTraits< AccumulateType >::Zero; - - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - sff += fixedValue * fixedValue; - smm += movingValue * movingValue; - sfm += fixedValue * movingValue; - if ( this->m_SubtractMean ) { - sf += fixedValue; - sm += movingValue; - } - this->m_NumberOfPixelsCounted++; - } - - ++ti; - } - - if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sff -= ( sf * sf / this->m_NumberOfPixelsCounted ); - smm -= ( sm * sm / this->m_NumberOfPixelsCounted ); - sfm -= ( sf * sm / this->m_NumberOfPixelsCounted ); - } - - const RealType denom = -1.0 * vcl_sqrt(sff * smm ); - - if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) { - measure = sfm / denom; - } else { - measure = itk::NumericTraits< MeasureType >::Zero; - } - - return measure; - -} - - - - - -/* - * Get the Derivative Measure - */ -template < class TFixedImage, class TMovingImage> -void -NormalizedCorrelationImageToImageMetric -::GetDerivative( const TransformParametersType & parameters, - DerivativeType & derivative ) const -{ - - 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 dimension = FixedImageType::ImageDimension; - - typedef itk::ImageRegionConstIteratorWithIndex FixedIteratorType; - - FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); - - typename FixedImageType::IndexType index; - - this->m_NumberOfPixelsCounted = 0; - - this->SetTransformParameters( parameters ); - - typedef typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType; - - AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sf = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sm = itk::NumericTraits< AccumulateType >::Zero; - - const unsigned int ParametersDimension = this->GetNumberOfParameters(); - derivative = DerivativeType( ParametersDimension ); - derivative.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeF = DerivativeType( ParametersDimension ); - derivativeF.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeM = DerivativeType( ParametersDimension ); - derivativeM.Fill( itk::NumericTraits::Zero ); - - ti.GoToBegin(); - // First compute the sums - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - sff += fixedValue * fixedValue; - smm += movingValue * movingValue; - sfm += fixedValue * movingValue; - if ( this->m_SubtractMean ) { - sf += fixedValue; - sm += movingValue; - } - this->m_NumberOfPixelsCounted++; - } - - ++ti; - } - - // Compute contributions to derivatives - ti.GoToBegin(); - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if ( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - -#if ITK_VERSION_MAJOR >= 4 - TransformJacobianType jacobian; - this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian ); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif - - // Get the gradient by NearestNeighboorInterpolation: - // which is equivalent to round up the point components. - typedef typename OutputPointType::CoordRepType CoordRepType; - typedef itk::ContinuousIndex - MovingImageContinuousIndexType; - - MovingImageContinuousIndexType tempIndex; - this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex ); - - typename MovingImageType::IndexType mappedIndex; - mappedIndex.CopyWithRound( tempIndex ); - - const GradientPixelType gradient = - this->GetGradientImage()->GetPixel( mappedIndex ); - - for(unsigned int par=0; par::Zero; - RealType sumM = itk::NumericTraits< RealType >::Zero; - for(unsigned int dim=0; dimm_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sumF -= differential * sf / this->m_NumberOfPixelsCounted; - sumM -= differential * sm / this->m_NumberOfPixelsCounted; - } - } - derivativeF[par] += sumF; - derivativeM[par] += sumM; - } - } - - ++ti; - } - - if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sff -= ( sf * sf / this->m_NumberOfPixelsCounted ); - smm -= ( sm * sm / this->m_NumberOfPixelsCounted ); - sfm -= ( sf * sm / this->m_NumberOfPixelsCounted ); - } - - const RealType denom = -1.0 * vcl_sqrt(sff * smm ); - - if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) { - for(unsigned int i=0; i::Zero; - } - } - -} - - -/* - * Get both the match Measure and theDerivative Measure - */ -template -void -NormalizedCorrelationImageToImageMetric -::GetValueAndDerivative(const TransformParametersType & parameters, - MeasureType & value, DerivativeType & derivative) const -{ - - - 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 dimension = FixedImageType::ImageDimension; - - typedef itk::ImageRegionConstIteratorWithIndex FixedIteratorType; - - FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); - - typename FixedImageType::IndexType index; - - this->m_NumberOfPixelsCounted = 0; - - this->SetTransformParameters( parameters ); - - typedef typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType; - - AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sf = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sm = itk::NumericTraits< AccumulateType >::Zero; - - const unsigned int ParametersDimension = this->GetNumberOfParameters(); - derivative = DerivativeType( ParametersDimension ); - derivative.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeF = DerivativeType( ParametersDimension ); - derivativeF.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeM = DerivativeType( ParametersDimension ); - derivativeM.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeM1 = DerivativeType( ParametersDimension ); - derivativeM1.Fill( itk::NumericTraits::Zero ); - - ti.GoToBegin(); - // First compute the sums - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - sff += fixedValue * fixedValue; - smm += movingValue * movingValue; - sfm += fixedValue * movingValue; - if ( this->m_SubtractMean ) { - sf += fixedValue; - sm += movingValue; - } - this->m_NumberOfPixelsCounted++; - } - - ++ti; - } - - - // Compute contributions to derivatives - ti.GoToBegin(); - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - -#if ITK_VERSION_MAJOR >= 4 - TransformJacobianType jacobian; - this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian ); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif - - // Get the gradient by NearestNeighboorInterpolation: - // which is equivalent to round up the point components. - typedef typename OutputPointType::CoordRepType CoordRepType; - typedef itk::ContinuousIndex - MovingImageContinuousIndexType; - - MovingImageContinuousIndexType tempIndex; - this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex ); - - typename MovingImageType::IndexType mappedIndex; - mappedIndex.CopyWithRound( tempIndex ); - - const GradientPixelType gradient = - this->GetGradientImage()->GetPixel( mappedIndex ); - - for(unsigned int par=0; par::Zero; - RealType sumM = itk::NumericTraits< RealType >::Zero; - for(unsigned int dim=0; dimm_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sumF -= differential * sf / this->m_NumberOfPixelsCounted; - sumM -= differential * sm / this->m_NumberOfPixelsCounted; - } - } - derivativeF[par] += sumF; - derivativeM[par] += sumM; - } - } - ++ti; - } - - if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sff -= ( sf * sf / this->m_NumberOfPixelsCounted ); - smm -= ( sm * sm / this->m_NumberOfPixelsCounted ); - sfm -= ( sf * sm / this->m_NumberOfPixelsCounted ); - } - - const RealType denom = -1.0 * vcl_sqrt(sff * smm ); - - if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) { - for(unsigned int i=0; i::Zero; - } - value = itk::NumericTraits< MeasureType >::Zero; - } - - - -} - -template < class TFixedImage, class TMovingImage> -void -NormalizedCorrelationImageToImageMetric -::PrintSelf(std::ostream& os, itk::Indent indent) const -{ - Superclass::PrintSelf(os, indent); - os << indent << "SubtractMean: " << m_SubtractMean << std::endl; -} - -} // end namespace itk - - -#endif // opt - #endif // _clitkNormalizedCorrelationImageToImageMetric.txx diff --git a/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h b/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h index 7a7d2b1..ad4a118 100644 --- a/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h +++ b/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h @@ -24,96 +24,7 @@ // gets integrated into the main directories. #include "itkConfigure.h" -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h" -#else - -#include "itkImageToImageMetric.h" -#include "itkCovariantVector.h" -#include "itkPoint.h" - - -namespace clitk -{ - -template < class TFixedImage, class TMovingImage > -class ITK_EXPORT NormalizedCorrelationImageToImageMetricFor3DBLUTFFD : - public itk::ImageToImageMetric< TFixedImage, TMovingImage> -{ -public: - - /** Standard class typedefs. */ - typedef NormalizedCorrelationImageToImageMetricFor3DBLUTFFD Self; - typedef itk::ImageToImageMetric Superclass; - - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(NormalizedCorrelationImageToImageMetricFor3DBLUTFFD, itk::Object); - - - /** Types transferred from the base class */ - typedef typename Superclass::RealType RealType; - typedef typename Superclass::TransformType TransformType; - typedef typename Superclass::TransformPointer TransformPointer; - typedef typename Superclass::TransformParametersType TransformParametersType; - typedef typename Superclass::TransformJacobianType TransformJacobianType; - typedef typename Superclass::GradientPixelType GradientPixelType; - typedef typename Superclass::OutputPointType OutputPointType; - typedef typename Superclass::InputPointType InputPointType; - - typedef typename Superclass::MeasureType MeasureType; - typedef typename Superclass::DerivativeType DerivativeType; - typedef typename Superclass::FixedImageType FixedImageType; - typedef typename Superclass::MovingImageType MovingImageType; - typedef typename Superclass::FixedImageConstPointer FixedImageConstPointer; - typedef typename Superclass::MovingImageConstPointer MovingImageConstPointer; - - - /** Get the derivatives of the match measure. */ - void GetDerivative( const TransformParametersType & parameters, - DerivativeType & Derivative ) const; - - /** Get the value for single valued optimizers. */ - MeasureType GetValue( const TransformParametersType & parameters ) const; - - /** Get value and derivatives for multiple valued optimizers. */ - void GetValueAndDerivative( const TransformParametersType & parameters, - MeasureType& Value, DerivativeType& Derivative ) const; - - /** Set/Get SubtractMean boolean. If true, the sample mean is subtracted - * from the sample values in the cross-correlation formula and - * typically results in narrower valleys in the cost fucntion. - * Default value is false. */ - itkSetMacro( SubtractMean, bool ); - itkGetConstReferenceMacro( SubtractMean, bool ); - itkBooleanMacro( SubtractMean ); - -protected: - NormalizedCorrelationImageToImageMetricFor3DBLUTFFD(); - virtual ~NormalizedCorrelationImageToImageMetricFor3DBLUTFFD() {}; - void PrintSelf(std::ostream& os, itk::Indent indent) const; - -private: - NormalizedCorrelationImageToImageMetricFor3DBLUTFFD(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented - - bool m_SubtractMean; - -}; - -} // end namespace clitk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx" -#endif - -#endif // opt - #endif // _clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx diff --git a/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx b/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx index 16ecec4..9dd366a 100644 --- a/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx +++ b/registration/clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx @@ -23,493 +23,5 @@ // This line can be removed once the optimized versions // gets integrated into the main directories. #include "itkConfigure.h" - -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx" -#else - - -#include "clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h" - -#include "itkImageRegionConstIteratorWithIndex.h" - - -namespace clitk -{ - -/* - * Constructor - */ -template -NormalizedCorrelationImageToImageMetricFor3DBLUTFFD -::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD() -{ - m_SubtractMean = false; -} - -/* - * Get the match Measure - */ -template -typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD::MeasureType -NormalizedCorrelationImageToImageMetricFor3DBLUTFFD -::GetValue( const TransformParametersType & parameters ) const -{ - - FixedImageConstPointer fixedImage = this->m_FixedImage; - - if( !fixedImage ) { - itkExceptionMacro( << "Fixed image has not been assigned" ); - } - - typedef itk::ImageRegionConstIteratorWithIndex FixedIteratorType; - - FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); - - typename FixedImageType::IndexType index; - - MeasureType measure; - - this->m_NumberOfPixelsCounted = 0; - - this->SetTransformParameters( parameters ); - - typedef typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType; - - AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sf = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sm = itk::NumericTraits< AccumulateType >::Zero; - - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - sff += fixedValue * fixedValue; - smm += movingValue * movingValue; - sfm += fixedValue * movingValue; - if ( this->m_SubtractMean ) { - sf += fixedValue; - sm += movingValue; - } - this->m_NumberOfPixelsCounted++; - } - - ++ti; - } - - if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sff -= ( sf * sf / this->m_NumberOfPixelsCounted ); - smm -= ( sm * sm / this->m_NumberOfPixelsCounted ); - sfm -= ( sf * sm / this->m_NumberOfPixelsCounted ); - } - - const RealType denom = -1.0 * vcl_sqrt(sff * smm ); - - if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) { - measure = sfm / denom; - } else { - measure = itk::NumericTraits< MeasureType >::Zero; - } - - return measure; - -} - - - - - -/* - * Get the Derivative Measure - */ -template < class TFixedImage, class TMovingImage> -void -NormalizedCorrelationImageToImageMetricFor3DBLUTFFD -::GetDerivative( const TransformParametersType & parameters, - DerivativeType & derivative ) const -{ - - 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 dimension = FixedImageType::ImageDimension; - - typedef itk::ImageRegionConstIteratorWithIndex FixedIteratorType; - - FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); - - typename FixedImageType::IndexType index; - - this->m_NumberOfPixelsCounted = 0; - - this->SetTransformParameters( parameters ); - - typedef typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType; - - AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sf = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sm = itk::NumericTraits< AccumulateType >::Zero; - - const unsigned int ParametersDimension = this->GetNumberOfParameters(); - derivative = DerivativeType( ParametersDimension ); - derivative.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeF = DerivativeType( ParametersDimension ); - derivativeF.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeM = DerivativeType( ParametersDimension ); - derivativeM.Fill( itk::NumericTraits::Zero ); - - ti.GoToBegin(); - // First compute the sums - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - sff += fixedValue * fixedValue; - smm += movingValue * movingValue; - sfm += fixedValue * movingValue; - if ( this->m_SubtractMean ) { - sf += fixedValue; - sm += movingValue; - } - this->m_NumberOfPixelsCounted++; - } - - ++ti; - } - - // Compute contributions to derivatives - ti.GoToBegin(); - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if ( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - -#if ITK_VERSION_MAJOR >= 4 - TransformJacobianType jacobian; - this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian ); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif - - // Get the gradient by NearestNeighboorInterpolation: - // which is equivalent to round up the point components. - typedef typename OutputPointType::CoordRepType CoordRepType; - typedef itk::ContinuousIndex - MovingImageContinuousIndexType; - - MovingImageContinuousIndexType tempIndex; - this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex ); - - typename MovingImageType::IndexType mappedIndex; - mappedIndex.CopyWithRound( tempIndex ); - - const GradientPixelType gradient = - this->GetGradientImage()->GetPixel( mappedIndex ); - - for(unsigned int par=0; par::Zero; - RealType sumM = itk::NumericTraits< RealType >::Zero; - for(unsigned int dim=0; dimm_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sumF -= differential * sf / this->m_NumberOfPixelsCounted; - sumM -= differential * sm / this->m_NumberOfPixelsCounted; - } - } - derivativeF[par] += sumF; - derivativeM[par] += sumM; - } - } - - ++ti; - } - - if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sff -= ( sf * sf / this->m_NumberOfPixelsCounted ); - smm -= ( sm * sm / this->m_NumberOfPixelsCounted ); - sfm -= ( sf * sm / this->m_NumberOfPixelsCounted ); - } - - const RealType denom = -1.0 * vcl_sqrt(sff * smm ); - - if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) { - for(unsigned int i=0; i::Zero; - } - } - -} - - -/* - * Get both the match Measure and theDerivative Measure - */ -template -void -NormalizedCorrelationImageToImageMetricFor3DBLUTFFD -::GetValueAndDerivative(const TransformParametersType & parameters, - MeasureType & value, DerivativeType & derivative) const -{ - - - 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 dimension = FixedImageType::ImageDimension; - - typedef itk::ImageRegionConstIteratorWithIndex FixedIteratorType; - - FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() ); - - typename FixedImageType::IndexType index; - - this->m_NumberOfPixelsCounted = 0; - - this->SetTransformParameters( parameters ); - - typedef typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType; - - AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sf = itk::NumericTraits< AccumulateType >::Zero; - AccumulateType sm = itk::NumericTraits< AccumulateType >::Zero; - - const unsigned int ParametersDimension = this->GetNumberOfParameters(); - derivative = DerivativeType( ParametersDimension ); - derivative.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeF = DerivativeType( ParametersDimension ); - derivativeF.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeM = DerivativeType( ParametersDimension ); - derivativeM.Fill( itk::NumericTraits::Zero ); - - DerivativeType derivativeM1 = DerivativeType( ParametersDimension ); - derivativeM1.Fill( itk::NumericTraits::Zero ); - - ti.GoToBegin(); - // First compute the sums - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - sff += fixedValue * fixedValue; - smm += movingValue * movingValue; - sfm += fixedValue * movingValue; - if ( this->m_SubtractMean ) { - sf += fixedValue; - sm += movingValue; - } - this->m_NumberOfPixelsCounted++; - } - - ++ti; - } - - - // Compute contributions to derivatives - ti.GoToBegin(); - while(!ti.IsAtEnd()) { - - index = ti.GetIndex(); - - InputPointType inputPoint; - fixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++ti; - continue; - } - - 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 RealType fixedValue = ti.Get(); - -#if ITK_VERSION_MAJOR >= 4 - TransformJacobianType jacobian; - this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian ); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif - - // Get the gradient by NearestNeighboorInterpolation: - // which is equivalent to round up the point components. - typedef typename OutputPointType::CoordRepType CoordRepType; - typedef itk::ContinuousIndex - MovingImageContinuousIndexType; - - MovingImageContinuousIndexType tempIndex; - this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex ); - - typename MovingImageType::IndexType mappedIndex; - mappedIndex.CopyWithRound( tempIndex ); - - const GradientPixelType gradient = - this->GetGradientImage()->GetPixel( mappedIndex ); - - for(unsigned int par=0; par::Zero; - RealType sumM = itk::NumericTraits< RealType >::Zero; - for(unsigned int dim=0; dimm_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sumF -= differential * sf / this->m_NumberOfPixelsCounted; - sumM -= differential * sm / this->m_NumberOfPixelsCounted; - } - } - derivativeF[par] += sumF; - derivativeM[par] += sumM; - } - } - ++ti; - } - - if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) { - sff -= ( sf * sf / this->m_NumberOfPixelsCounted ); - smm -= ( sm * sm / this->m_NumberOfPixelsCounted ); - sfm -= ( sf * sm / this->m_NumberOfPixelsCounted ); - } - - const RealType denom = -1.0 * vcl_sqrt(sff * smm ); - - if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) { - for(unsigned int i=0; i::Zero; - } - value = itk::NumericTraits< MeasureType >::Zero; - } - - - -} - -template < class TFixedImage, class TMovingImage> -void -NormalizedCorrelationImageToImageMetricFor3DBLUTFFD -::PrintSelf(std::ostream& os, itk::Indent indent) const -{ - Superclass::PrintSelf(os, indent); - os << indent << "SubtractMean: " << m_SubtractMean << std::endl; -} - -} // end namespace itk - - -#endif // opt - #endif // _clitkNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx diff --git a/registration/clitkOptNormalizedCorrelationImageToImageMetric.h b/registration/clitkOptNormalizedCorrelationImageToImageMetric.h index 83448af..242a4a1 100644 --- a/registration/clitkOptNormalizedCorrelationImageToImageMetric.h +++ b/registration/clitkOptNormalizedCorrelationImageToImageMetric.h @@ -19,11 +19,7 @@ #ifndef __clitkOptNormalizedCorrelationImageToImageMetric_h #define __clitkOptNormalizedCorrelationImageToImageMetric_h -#if ITK_VERSION_MAJOR >= 4 - #include "itkImageToImageMetric.h" -#else - #include "itkOptImageToImageMetric.h" -#endif +#include "itkImageToImageMetric.h" #include "itkCovariantVector.h" #include "itkPoint.h" #include "itkIndex.h" diff --git a/registration/clitkOptNormalizedCorrelationImageToImageMetric.txx b/registration/clitkOptNormalizedCorrelationImageToImageMetric.txx index 3c0dfab..0f8fa64 100644 --- a/registration/clitkOptNormalizedCorrelationImageToImageMetric.txx +++ b/registration/clitkOptNormalizedCorrelationImageToImageMetric.txx @@ -219,9 +219,6 @@ NormalizedCorrelationImageToImageMetric // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // MUST BE CALLED TO INITIATE PROCESSING this->GetValueMultiThreadedInitiate(); @@ -296,9 +293,6 @@ NormalizedCorrelationImageToImageMetric // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // MUST BE CALLED TO INITIATE PROCESSING this->GetValueMultiThreadedInitiate(); @@ -385,13 +379,8 @@ NormalizedCorrelationImageToImageMetric } // Jacobian should be evaluated at the unmapped (fixed image) point. -#if ITK_VERSION_MAJOR >= 4 TransformJacobianType jacobian; transform->ComputeJacobianWithRespectToParameters(fixedImagePoint, jacobian); -#else - const TransformJacobianType & jacobian = transform - ->GetJacobian( fixedImagePoint ); -#endif for(unsigned int par=0; parm_NumberOfParameters; par++) { RealType sumF = itk::NumericTraits< RealType >::Zero; @@ -431,9 +420,6 @@ NormalizedCorrelationImageToImageMetric // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif //We need the sums and the value to be calculated first value=this->ComputeSums(parameters); diff --git a/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h b/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h index 65559c2..dc7ef9f 100644 --- a/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h +++ b/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.h @@ -19,11 +19,7 @@ #ifndef __clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD_h #define __clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD_h -#if ITK_VERSION_MAJOR >= 4 - #include "itkImageToImageMetric.h" -#else - #include "itkOptImageToImageMetric.h" -#endif +#include "itkImageToImageMetric.h" #include "itkCovariantVector.h" #include "itkPoint.h" #include "itkIndex.h" diff --git a/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx b/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx index d7c7777..959bfed 100644 --- a/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx +++ b/registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx @@ -219,9 +219,6 @@ NormalizedCorrelationImageToImageMetricFor3DBLUTFFD // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // MUST BE CALLED TO INITIATE PROCESSING this->GetValueMultiThreadedInitiate(); @@ -296,9 +293,6 @@ NormalizedCorrelationImageToImageMetricFor3DBLUTFFD // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // MUST BE CALLED TO INITIATE PROCESSING this->GetValueMultiThreadedInitiate(); @@ -385,12 +379,8 @@ NormalizedCorrelationImageToImageMetricFor3DBLUTFFD } // Jacobian should be evaluated at the unmapped (fixed image) point. -#if ITK_VERSION_MAJOR >= 4 TransformJacobianType jacobian; transform->ComputeJacobianWithRespectToParameters( fixedImagePoint, jacobian ); -#else - const TransformJacobianType & jacobian = transform->GetJacobian( fixedImagePoint ); -#endif // for(unsigned int par=0; parm_NumberOfParameters; par++) // { @@ -455,9 +445,6 @@ NormalizedCorrelationImageToImageMetricFor3DBLUTFFD // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif //We need the sums and the value to be calculated first value=this->ComputeSums(parameters); diff --git a/registration/clitkPointListTransform.h b/registration/clitkPointListTransform.h index d98c525..b897f12 100644 --- a/registration/clitkPointListTransform.h +++ b/registration/clitkPointListTransform.h @@ -117,7 +117,6 @@ namespace clitk return OutputCovariantVectorType(); } -#if ITK_VERSION_MAJOR >= 4 virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const { itkExceptionMacro( << "PointListTransform doesn't declare ComputeJacobianWithRespectToParameters" ); @@ -126,13 +125,6 @@ namespace clitk { itkExceptionMacro( << "PointListTransform doesn't declare ComputeJacobianWithRespectToPosition" ); } -#else - virtual const JacobianType& GetJacobian(const InputPointType &point ) const - { - itkExceptionMacro( << "PointListTransform doesn't declare GetJacobian" ); - return this->m_Jacobian; - } -#endif protected: PointListTransform(); diff --git a/registration/clitkPointListTransform.txx b/registration/clitkPointListTransform.txx index a93a309..d827593 100644 --- a/registration/clitkPointListTransform.txx +++ b/registration/clitkPointListTransform.txx @@ -25,12 +25,7 @@ namespace clitk // Constructor template - PointListTransform -#if ITK_VERSION_MAJOR >= 4 - ::PointListTransform():Superclass(1) -#else - ::PointListTransform():Superclass(NOutputDimensions,1) -#endif + PointListTransform::PointListTransform():Superclass(1) { m_PointLists.resize(0); m_PointList.resize(1); diff --git a/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.h b/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.h index b192458..309b867 100644 --- a/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.h +++ b/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.h @@ -68,9 +68,7 @@ namespace clitk /** Standard parameters container. */ typedef typename Superclass::ParametersType ParametersType; -#if ITK_VERSION_MAJOR >= 4 typedef typename Superclass::NumberOfParametersType NumberOfParametersType; -#endif /** Standard Jacobian container. */ typedef typename Superclass::JacobianType JacobianType; @@ -266,22 +264,14 @@ namespace clitk } /** Compute the Jacobian Matrix of the transformation at one point */ -#if ITK_VERSION_MAJOR >= 4 virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const; virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const { itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() ); } -#else - virtual const JacobianType& GetJacobian(const InputPointType &point ) const; -#endif /** Return the number of parameters that completely define the Transfom */ -#if ITK_VERSION_MAJOR >= 4 virtual NumberOfParametersType GetNumberOfParameters(void) const; -#else - virtual unsigned int GetNumberOfParameters(void) const; -#endif //JV Return the padded number of parameters virtual unsigned int GetPaddedNumberOfParameters(void) const; @@ -445,9 +435,7 @@ namespace clitk // JV Shape unsigned int m_TransformShape; -#if ITK_VERSION_MAJOR >= 4 mutable JacobianType m_SharedDataBSplineJacobian; -#endif }; //class ShapedBLUTSpatioTemporalDeformableTransform diff --git a/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.txx b/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.txx index a5db85a..7cc0107 100644 --- a/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.txx +++ b/registration/clitkShapedBLUTSpatioTemporalDeformableTransform.txx @@ -31,12 +31,7 @@ namespace clitk // Constructor with default arguments template - ShapedBLUTSpatioTemporalDeformableTransform -#if ITK_VERSION_MAJOR >= 4 - ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0) -#else - ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(OutputDimension,0) -#endif + ShapedBLUTSpatioTemporalDeformableTransform::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0) { unsigned int i; @@ -383,11 +378,7 @@ namespace clitk // Get the number of parameters template -#if ITK_VERSION_MAJOR >= 4 typename ShapedBLUTSpatioTemporalDeformableTransform::NumberOfParametersType -#else - unsigned int -#endif ShapedBLUTSpatioTemporalDeformableTransform ::GetNumberOfParameters(void) const { @@ -810,19 +801,11 @@ namespace clitk //===================================== //JV Wrap jacobian into OutputDimension X Vectorial images //===================================== -#if ITK_VERSION_MAJOR >= 4 this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); -#else - this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); -#endif // Use memset to set the memory // JV four rows of three comps of parameters -#if ITK_VERSION_MAJOR >= 4 JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_SharedDataBSplineJacobian.data_block()); -#else - JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_Jacobian.data_block()); -#endif memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType)); for (unsigned int j=0; j -#if ITK_VERSION_MAJOR >= 4 void ShapedBLUTSpatioTemporalDeformableTransform ::ComputeJacobianWithRespectToParameters( const InputPointType & point, JacobianType & jacobian) const -#else - const - typename ShapedBLUTSpatioTemporalDeformableTransform - ::JacobianType & - ShapedBLUTSpatioTemporalDeformableTransform - ::GetJacobian( const InputPointType & point ) const -#endif { //======================================================== @@ -2496,12 +2471,8 @@ namespace clitk if(m_Mask && !(m_Mask->IsInside(point) ) ) { // Outside: no (deformable) displacement -#if ITK_VERSION_MAJOR >= 4 jacobian = m_SharedDataBSplineJacobian; return; -#else - return this->m_Jacobian; -#endif } // Get index @@ -2511,12 +2482,8 @@ namespace clitk // we assume zero displacement and return the input point if ( !this->InsideValidRegion( m_Index ) ) { -#if ITK_VERSION_MAJOR >= 4 jacobian = m_SharedDataBSplineJacobian; return; -#else - return this->m_Jacobian; -#endif } // Compute interpolation weights @@ -2684,11 +2651,7 @@ namespace clitk } // Return the result -#if ITK_VERSION_MAJOR >= 4 jacobian = m_SharedDataBSplineJacobian; -#else - return this->m_Jacobian; -#endif } diff --git a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h index c941b05..5342fea 100644 --- a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h +++ b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h @@ -39,497 +39,5 @@ // This line can be removed once the optimized versions // gets integrated into the main directories. #include "itkConfigure.h" - -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h" -#else - -#include "itkImageToImageMetric.h" -#include "itkCovariantVector.h" -#include "itkPoint.h" -#include "itkIndex.h" -#include "itkBSplineKernelFunction.h" -#include "itkBSplineDerivativeKernelFunction.h" -#include "itkCentralDifferenceImageFunction.h" -#include "itkBSplineInterpolateImageFunction.h" -#include "itkBSplineDeformableTransform.h" -#include "itkArray2D.h" - -namespace itk -{ - -/** \class MattesMutualInformationImageToImageMetricFor3DBLUTFFD - * \brief Computes the mutual information between two images to be - * registered using the method of Mattes et al. - * - * MattesMutualInformationImageToImageMetricFor3DBLUTFFD computes the mutual - * information between a fixed and moving image to be registered. - * - * This class is templated over the FixedImage type and the MovingImage - * type. - * - * The fixed and moving images are set via methods SetFixedImage() and - * SetMovingImage(). This metric makes use of user specified Transform and - * Interpolator. The Transform is used to map points from the fixed image to - * the moving image domain. The Interpolator is used to evaluate the image - * intensity at user specified geometric points in the moving image. - * The Transform and Interpolator are set via methods SetTransform() and - * SetInterpolator(). - * - * If a BSplineInterpolationFunction is used, this class obtain - * image derivatives from the BSpline interpolator. Otherwise, - * image derivatives are computed using central differencing. - * - * \warning This metric assumes that the moving image has already been - * connected to the interpolator outside of this class. - * - * The method GetValue() computes of the mutual information - * while method GetValueAndDerivative() computes - * both the mutual information and its derivatives with respect to the - * transform parameters. - * - * The calculations are based on the method of Mattes et al [1,2] - * where the probability density distribution are estimated using - * Parzen histograms. Since the fixed image PDF does not contribute - * to the derivatives, it does not need to be smooth. Hence, - * a zero order (box car) BSpline kernel is used - * for the fixed image intensity PDF. On the other hand, to ensure - * smoothness a third order BSpline kernel is used for the - * moving image intensity PDF. - * - * On Initialize(), the FixedImage is uniformly sampled within - * the FixedImageRegion. The number of samples used can be set - * via SetNumberOfSpatialSamples(). Typically, the number of - * spatial samples used should increase with the image size. - * - * The option UseAllPixelOn() disables the random sampling and uses - * all the pixels of the FixedImageRegion in order to estimate the - * joint intensity PDF. - * - * During each call of GetValue(), GetDerivatives(), - * GetValueAndDerivatives(), marginal and joint intensity PDF's - * values are estimated at discrete position or bins. - * The number of bins used can be set via SetNumberOfHistogramBins(). - * To handle data with arbitray magnitude and dynamic range, - * the image intensity is scale such that any contribution to the - * histogram will fall into a valid bin. - * - * One the PDF's have been contructed, the mutual information - * is obtained by doubling summing over the discrete PDF values. - * - * - * Notes: - * 1. This class returns the negative mutual information value. - * 2. This class in not thread safe due the private data structures - * used to the store the sampled points and the marginal and joint pdfs. - * - * References: - * [1] "Nonrigid multimodality image registration" - * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank - * Medical Imaging 2001: Image Processing, 2001, pp. 1609-1620. - * [2] "PET-CT Image Registration in the Chest Using Free-form Deformations" - * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank - * IEEE Transactions in Medical Imaging. Vol.22, No.1, - January 2003. pp.120-128. - * [3] "Optimization of Mutual Information for MultiResolution Image - * Registration" - * P. Thevenaz and M. Unser - * IEEE Transactions in Image Processing, 9(12) December 2000. - * - * \ingroup RegistrationMetrics - * \ingroup ThreadUnSafe - */ -template -class ITK_EXPORT MattesMutualInformationImageToImageMetricFor3DBLUTFFD : - public ImageToImageMetric< TFixedImage, TMovingImage > -{ -public: - - /** Standard class typedefs. */ - typedef MattesMutualInformationImageToImageMetricFor3DBLUTFFD Self; - typedef ImageToImageMetric< TFixedImage, TMovingImage > Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(MattesMutualInformationImageToImageMetricFor3DBLUTFFD, ImageToImageMetric); - - /** Types inherited from Superclass. */ - typedef typename Superclass::TransformType TransformType; - typedef typename Superclass::TransformPointer TransformPointer; - typedef typename Superclass::TransformJacobianType TransformJacobianType; - typedef typename Superclass::InterpolatorType InterpolatorType; - typedef typename Superclass::MeasureType MeasureType; - typedef typename Superclass::DerivativeType DerivativeType; - typedef typename Superclass::ParametersType ParametersType; - typedef typename Superclass::FixedImageType FixedImageType; - typedef typename Superclass::MovingImageType MovingImageType; - typedef typename Superclass::FixedImageConstPointer FixedImageConstPointer; - typedef typename Superclass::MovingImageConstPointer MovingImageCosntPointer; - typedef typename Superclass::InputPointType InputPointType; - typedef typename Superclass::OutputPointType OutputPointType; - - typedef typename Superclass::CoordinateRepresentationType - CoordinateRepresentationType; - - /** Index and Point typedef support. */ - typedef typename FixedImageType::IndexType FixedImageIndexType; - typedef typename FixedImageIndexType::IndexValueType FixedImageIndexValueType; - typedef typename MovingImageType::IndexType MovingImageIndexType; - typedef typename TransformType::InputPointType FixedImagePointType; - typedef typename TransformType::OutputPointType MovingImagePointType; - - /** The moving image dimension. */ - itkStaticConstMacro( MovingImageDimension, unsigned int, - MovingImageType::ImageDimension ); - - /** - * Initialize the Metric by - * (1) making sure that all the components are present and plugged - * together correctly, - * (2) uniformly select NumberOfSpatialSamples within - * the FixedImageRegion, and - * (3) allocate memory for pdf data structures. */ - virtual void Initialize(void) throw ( ExceptionObject ); - - /** Get the derivatives of the match measure. */ - void GetDerivative( const ParametersType& parameters, - DerivativeType & Derivative ) const; - - /** Get the value. */ - MeasureType GetValue( const ParametersType& parameters ) const; - - /** Get the value and derivatives for single valued optimizers. */ - void GetValueAndDerivative( const ParametersType& parameters, - MeasureType& Value, - DerivativeType& Derivative ) const; - - /** Number of spatial samples to used to compute metric */ - itkSetClampMacro( NumberOfSpatialSamples, unsigned long, - 1, NumericTraits::max() ); - itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long); - - /** Number of bins to used in the histogram. Typical value is 50. */ - itkSetClampMacro( NumberOfHistogramBins, unsigned long, - 1, NumericTraits::max() ); - itkGetConstReferenceMacro( NumberOfHistogramBins, unsigned long); - - /** Reinitialize the seed of the random number generator that selects the - * sample of pixels used for estimating the image histograms and the joint - * histogram. By nature, this metric is not deterministic, since at each run - * it may select a different set of pixels. By initializing the random number - * generator seed to the same value you can restore determinism. On the other - * hand, calling the method ReinitializeSeed() without arguments will use the - * clock from your machine in order to have a very random initialization of - * the seed. This will indeed increase the non-deterministic behavior of the - * metric. */ - void ReinitializeSeed(); - void ReinitializeSeed(int); - - /** Select whether the metric will be computed using all the pixels on the - * fixed image region, or only using a set of randomly selected pixels. */ - itkSetMacro(UseAllPixels,bool); - itkGetConstReferenceMacro(UseAllPixels,bool); - itkBooleanMacro(UseAllPixels); - - /** This variable selects the method to be used for computing the Metric - * derivatives with respect to the Transform parameters. Two modes of - * computation are available. The choice between one and the other is a - * trade-off between computation speed and memory allocations. The two modes - * are described in detail below: - * - * UseExplicitPDFDerivatives = True - * will compute the Metric derivative by first calculating the derivatives of - * each one of the Joint PDF bins with respect to each one of the Transform - * parameters and then accumulating these contributions in the final metric - * derivative array by using a bin-specific weight. The memory required for - * storing the intermediate derivatives is a 3D array of doubles with size - * equals to the product of (number of histogram bins)^2 times number of - * transform parameters. This method is well suited for Transform with a small - * number of parameters. - * - * UseExplicitPDFDerivatives = False will compute the Metric derivative by - * first computing the weights for each one of the Joint PDF bins and caching - * them into an array. Then it will revisit each one of the PDF bins for - * computing its weighted contribution to the full derivative array. In this - * method an extra 2D array is used for storing the weights of each one of - * the PDF bins. This is an array of doubles with size equals to (number of - * histogram bins)^2. This method is well suited for Transforms with a large - * number of parameters, such as, BSplineDeformableTransforms. */ - itkSetMacro(UseExplicitPDFDerivatives,bool); - itkGetConstReferenceMacro(UseExplicitPDFDerivatives,bool); - itkBooleanMacro(UseExplicitPDFDerivatives); - - /** This boolean flag is only relevant when this metric is used along - * with a BSplineDeformableTransform. The flag enables/disables the - * caching of values computed when a physical point is mapped through - * the BSplineDeformableTransform. In particular it will cache the - * values of the BSpline weights for that points, and the indexes - * indicating what BSpline-grid nodes are relevant for that specific - * point. This caching is made optional due to the fact that the - * memory arrays used for the caching can reach large sizes even for - * moderate image size problems. For example, for a 3D image of - * 256^3, using 20% of pixels, these arrays will take about 1 - * Gigabyte of RAM for storage. The ratio of computing time between - * using the cache or not using the cache can reach 1:5, meaning that - * using the caching can provide a five times speed up. It is - * therefore, interesting to enable the caching, if enough memory is - * available for it. The caching is enabled by default, in order to - * preserve backward compatibility with previous versions of ITK. */ - itkSetMacro(UseCachingOfBSplineWeights,bool); - itkGetConstReferenceMacro(UseCachingOfBSplineWeights,bool); - itkBooleanMacro(UseCachingOfBSplineWeights); - -protected: - - MattesMutualInformationImageToImageMetricFor3DBLUTFFD(); - virtual ~MattesMutualInformationImageToImageMetricFor3DBLUTFFD() {}; - void PrintSelf(std::ostream& os, Indent indent) const; - - /** \class FixedImageSpatialSample - * A fixed image spatial sample consists of the fixed domain point - * and the fixed image value at that point. */ - /// @cond - class FixedImageSpatialSample - { - public: - FixedImageSpatialSample():FixedImageValue(0.0) { - FixedImagePointValue.Fill(0.0); - } - ~FixedImageSpatialSample() {}; - - FixedImagePointType FixedImagePointValue; - double FixedImageValue; - unsigned int FixedImageParzenWindowIndex; - }; - /// @endcond - - /** FixedImageSpatialSample typedef support. */ - typedef std::vector - FixedImageSpatialSampleContainer; - - /** Container to store a set of points and fixed image values. */ - FixedImageSpatialSampleContainer m_FixedImageSamples; - - /** Uniformly select a sample set from the fixed image domain. */ - virtual void SampleFixedImageDomain( - FixedImageSpatialSampleContainer& samples); - - /** Gather all the pixels from the fixed image domain. */ - virtual void SampleFullFixedImageDomain( - FixedImageSpatialSampleContainer& samples); - - /** Transform a point from FixedImage domain to MovingImage domain. - * This function also checks if mapped point is within support region. */ - virtual void TransformPoint( unsigned int sampleNumber, - const ParametersType& parameters, - MovingImagePointType& mappedPoint, - bool& sampleWithinSupportRegion, - double& movingImageValue ) const; - -private: - - //purposely not implemented - MattesMutualInformationImageToImageMetricFor3DBLUTFFD(const Self&); - //purposely not implemented - void operator=(const Self&); - - - /** The marginal PDFs are stored as std::vector. */ - typedef float PDFValueType; - typedef std::vector MarginalPDFType; - - /** The fixed image marginal PDF. */ - mutable MarginalPDFType m_FixedImageMarginalPDF; - - /** The moving image marginal PDF. */ - mutable MarginalPDFType m_MovingImageMarginalPDF; - - /** Helper array for storing the values of the JointPDF ratios. */ - typedef double PRatioType; - typedef Array2D< PRatioType > PRatioArrayType; - mutable PRatioArrayType m_PRatioArray; - - /** Helper variable for accumulating the derivative of the metric. */ - mutable DerivativeType m_MetricDerivative; - - /** Typedef for the joint PDF and PDF derivatives are stored as ITK Images. */ - typedef Image JointPDFType; - typedef JointPDFType::IndexType JointPDFIndexType; - typedef JointPDFType::PixelType JointPDFValueType; - typedef JointPDFType::RegionType JointPDFRegionType; - typedef JointPDFType::SizeType JointPDFSizeType; - typedef Image JointPDFDerivativesType; - typedef JointPDFDerivativesType::IndexType JointPDFDerivativesIndexType; - typedef JointPDFDerivativesType::PixelType JointPDFDerivativesValueType; - typedef JointPDFDerivativesType::RegionType JointPDFDerivativesRegionType; - typedef JointPDFDerivativesType::SizeType JointPDFDerivativesSizeType; - - /** The joint PDF and PDF derivatives. */ - typename JointPDFType::Pointer m_JointPDF; - typename JointPDFDerivativesType::Pointer m_JointPDFDerivatives; - - unsigned long m_NumberOfSpatialSamples; - unsigned long m_NumberOfParameters; - - /** Variables to define the marginal and joint histograms. */ - unsigned long m_NumberOfHistogramBins; - double m_MovingImageNormalizedMin; - double m_FixedImageNormalizedMin; - double m_MovingImageTrueMin; - double m_MovingImageTrueMax; - double m_FixedImageBinSize; - double m_MovingImageBinSize; - - /** Typedefs for BSpline kernel and derivative functions. */ - typedef BSplineKernelFunction<3> CubicBSplineFunctionType; - typedef BSplineDerivativeKernelFunction<3> CubicBSplineDerivativeFunctionType; - - /** Cubic BSpline kernel for computing Parzen histograms. */ - typename CubicBSplineFunctionType::Pointer m_CubicBSplineKernel; - typename CubicBSplineDerivativeFunctionType::Pointer - m_CubicBSplineDerivativeKernel; - - /** Precompute fixed image parzen window indices. */ - virtual void ComputeFixedImageParzenWindowIndices( - FixedImageSpatialSampleContainer& samples ); - - /** - * Types and variables related to image derivative calculations. - * If a BSplineInterpolationFunction is used, this class obtain - * image derivatives from the BSpline interpolator. Otherwise, - * image derivatives are computed using central differencing. - */ - typedef CovariantVector< double, - itkGetStaticConstMacro(MovingImageDimension) > - ImageDerivativesType; - - /** Compute image derivatives at a point. */ - virtual void ComputeImageDerivatives( const MovingImagePointType& mappedPoint, - ImageDerivativesType& gradient ) const; - - /** Boolean to indicate if the interpolator BSpline. */ - bool m_InterpolatorIsBSpline; - - /** Typedefs for using BSpline interpolator. */ - typedef - BSplineInterpolateImageFunction - BSplineInterpolatorType; - - /** Pointer to BSplineInterpolator. */ - typename BSplineInterpolatorType::Pointer m_BSplineInterpolator; - - /** Typedefs for using central difference calculator. */ - typedef CentralDifferenceImageFunction - DerivativeFunctionType; - - /** Pointer to central difference calculator. */ - typename DerivativeFunctionType::Pointer m_DerivativeCalculator; - - - /** Compute PDF derivative contribution for each parameter. */ - virtual void ComputePDFDerivatives( unsigned int sampleNumber, - int movingImageParzenWindowIndex, - const ImageDerivativesType& - movingImageGradientValue, - double cubicBSplineDerivativeValue - ) const; - - /** - * Types and variables related to BSpline deformable transforms. - * If the transform is of type third order BSplineDeformableTransform, - * then we can speed up the metric derivative calculation by - * only inspecting the parameters within the support region - * of a mapped point. - */ - - /** Boolean to indicate if the transform is BSpline deformable. */ - bool m_TransformIsBSpline; - - /** The number of BSpline parameters per image dimension. */ - long m_NumParametersPerDim; - - /** - * The number of BSpline transform weights is the number of - * of parameter in the support region (per dimension ). */ - unsigned long m_NumBSplineWeights; - - /** The fixed image dimension. */ - itkStaticConstMacro( FixedImageDimension, unsigned int, - FixedImageType::ImageDimension ); - - /** - * Enum of the deformabtion field spline order. - */ - enum { DeformationSplineOrder = 3 }; - - /** - * Typedefs for the BSplineDeformableTransform. - */ - typedef BSplineDeformableTransform< - CoordinateRepresentationType, - ::itk::GetImageDimension::ImageDimension, - DeformationSplineOrder> BSplineTransformType; - typedef typename BSplineTransformType::WeightsType - BSplineTransformWeightsType; - typedef typename BSplineTransformType::ParameterIndexArrayType - BSplineTransformIndexArrayType; - - /** - * Variables used when transform is of type BSpline deformable. - */ - typename BSplineTransformType::Pointer m_BSplineTransform; - - /** - * Cache pre-transformed points, weights, indices and - * within support region flag. - */ - typedef typename BSplineTransformWeightsType::ValueType WeightsValueType; - typedef Array2D BSplineTransformWeightsArrayType; - typedef typename BSplineTransformIndexArrayType::ValueType IndexValueType; - typedef Array2D BSplineTransformIndicesArrayType; - typedef std::vector MovingImagePointArrayType; - typedef std::vector BooleanArrayType; - - BSplineTransformWeightsArrayType m_BSplineTransformWeightsArray; - BSplineTransformIndicesArrayType m_BSplineTransformIndicesArray; - MovingImagePointArrayType m_PreTransformPointsArray; - BooleanArrayType m_WithinSupportRegionArray; - - typedef FixedArray::ImageDimension> - ParametersOffsetType; - ParametersOffsetType m_ParametersOffset; - - bool m_UseAllPixels; - - virtual void PreComputeTransformValues(); - - bool m_ReseedIterator; - int m_RandomSeed; - - // Selection of explicit or implicit computation of PDF derivatives - // with respect to Transform parameters. - bool m_UseExplicitPDFDerivatives; - - // Variables needed for optionally caching values when using a BSpline transform. - bool m_UseCachingOfBSplineWeights; - mutable BSplineTransformWeightsType m_Weights; - mutable BSplineTransformIndexArrayType m_Indices; - -}; - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx" -#endif - -#endif - #endif diff --git a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx index 9391dc2..929d534 100644 --- a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx +++ b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx @@ -41,1555 +41,5 @@ // gets integrated into the main directories. #include "itkConfigure.h" -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx" -#else - - -#include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h" -#include "itkBSplineInterpolateImageFunction.h" -#include "itkCovariantVector.h" -#include "itkImageRandomConstIteratorWithIndex.h" -#include "itkImageRegionConstIterator.h" -#include "itkImageRegionIterator.h" -#include "itkImageIterator.h" -#include "vnl/vnl_math.h" -#include "itkBSplineDeformableTransform.h" - -namespace itk -{ - - -/** - * Constructor - */ -template < class TFixedImage, class TMovingImage > -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::MattesMutualInformationImageToImageMetricFor3DBLUTFFD() -{ - - m_NumberOfSpatialSamples = 500; - m_NumberOfHistogramBins = 50; - - this->SetComputeGradient(false); // don't use the default gradient for now - - m_InterpolatorIsBSpline = false; - m_TransformIsBSpline = false; - - // Initialize PDFs to NULL - m_JointPDF = NULL; - m_JointPDFDerivatives = NULL; - - m_UseExplicitPDFDerivatives = true; - - typename BSplineTransformType::Pointer transformer = - BSplineTransformType::New(); - this->SetTransform (transformer); - - typename BSplineInterpolatorType::Pointer interpolator = - BSplineInterpolatorType::New(); - this->SetInterpolator (interpolator); - - // Initialize memory - m_MovingImageNormalizedMin = 0.0; - m_FixedImageNormalizedMin = 0.0; - m_MovingImageTrueMin = 0.0; - m_MovingImageTrueMax = 0.0; - m_FixedImageBinSize = 0.0; - m_MovingImageBinSize = 0.0; - m_CubicBSplineDerivativeKernel = NULL; - m_BSplineInterpolator = NULL; - m_DerivativeCalculator = NULL; - m_NumParametersPerDim = 0; - m_NumBSplineWeights = 0; - m_BSplineTransform = NULL; - m_NumberOfParameters = 0; - m_UseAllPixels = false; - m_ReseedIterator = false; - m_RandomSeed = -1; - m_UseCachingOfBSplineWeights = true; -} - - -/** - * Print out internal information about this class - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::PrintSelf(std::ostream& os, Indent indent) const -{ - - Superclass::PrintSelf(os, indent); - - os << indent << "NumberOfSpatialSamples: "; - os << m_NumberOfSpatialSamples << std::endl; - os << indent << "NumberOfHistogramBins: "; - os << m_NumberOfHistogramBins << std::endl; - os << indent << "UseAllPixels: "; - os << m_UseAllPixels << std::endl; - - // Debugging information - os << indent << "NumberOfParameters: "; - os << m_NumberOfParameters << std::endl; - os << indent << "FixedImageNormalizedMin: "; - os << m_FixedImageNormalizedMin << std::endl; - os << indent << "MovingImageNormalizedMin: "; - os << m_MovingImageNormalizedMin << std::endl; - os << indent << "MovingImageTrueMin: "; - os << m_MovingImageTrueMin << std::endl; - os << indent << "MovingImageTrueMax: "; - os << m_MovingImageTrueMax << std::endl; - os << indent << "FixedImageBinSize: "; - os << m_FixedImageBinSize << std::endl; - os << indent << "MovingImageBinSize: "; - os << m_MovingImageBinSize << std::endl; - os << indent << "InterpolatorIsBSpline: "; - os << m_InterpolatorIsBSpline << std::endl; - os << indent << "TransformIsBSpline: "; - os << m_TransformIsBSpline << std::endl; - os << indent << "UseCachingOfBSplineWeights: "; - os << m_UseCachingOfBSplineWeights << std::endl; - os << indent << "UseExplicitPDFDerivatives: "; - os << m_UseExplicitPDFDerivatives << std::endl; - -} - - -/** - * Initialize - */ -template -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::Initialize(void) throw ( ExceptionObject ) -{ - this->Superclass::Initialize(); - - // Cache the number of transformation parameters - m_NumberOfParameters = this->m_Transform->GetNumberOfParameters(); - - /** - * Compute the minimum and maximum for the FixedImage over - * the FixedImageRegion. - * - * NB: 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( - this->m_FixedImage, this->GetFixedImageRegion() ); - - for ( fixedImageIterator.GoToBegin(); - !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) { - - double sample = static_cast( fixedImageIterator.Get() ); - - if ( sample < fixedImageMin ) { - fixedImageMin = sample; - } - - if ( sample > fixedImageMax ) { - fixedImageMax = sample; - } - } - - - /** - * Compute the minimum and maximum for the entire moving image - * in the buffer. - */ - double movingImageMin = NumericTraits::max(); - double movingImageMax = NumericTraits::NonpositiveMin(); - - typedef ImageRegionConstIterator MovingIteratorType; - MovingIteratorType movingImageIterator( - this->m_MovingImage, this->m_MovingImage->GetBufferedRegion() ); - - for ( movingImageIterator.GoToBegin(); - !movingImageIterator.IsAtEnd(); ++movingImageIterator) { - double sample = static_cast( movingImageIterator.Get() ); - - if ( sample < movingImageMin ) { - movingImageMin = sample; - } - - if ( sample > movingImageMax ) { - movingImageMax = sample; - } - } - - m_MovingImageTrueMin = movingImageMin; - m_MovingImageTrueMax = movingImageMax; - - itkDebugMacro( " FixedImageMin: " << fixedImageMin << - " FixedImageMax: " << fixedImageMax << std::endl ); - itkDebugMacro( " MovingImageMin: " << movingImageMin << - " MovingImageMax: " << movingImageMax << std::endl ); - - - /** - * Compute binsize for the histograms. - * - * The binsize for the image intensities needs to be adjusted so that - * we can avoid dealing with boundary conditions using the cubic - * spline as the Parzen window. We do this by increasing the size - * of the bins so that the joint histogram becomes "padded" at the - * borders. Because we are changing the binsize, - * we also need to shift the minimum by the padded amount in order to - * avoid minimum values filling in our padded region. - * - * Note that there can still be non-zero bin values in the padded region, - * it's just that these bins will never be a central bin for the Parzen - * window. - * - */ - const int padding = 2; // this will pad by 2 bins - - m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / - static_cast( m_NumberOfHistogramBins - 2 * padding ); - m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize - - static_cast( padding ); - - m_MovingImageBinSize = ( movingImageMax - movingImageMin ) / - static_cast( m_NumberOfHistogramBins - 2 * padding ); - m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize - - static_cast( padding ); - - - itkDebugMacro( "FixedImageNormalizedMin: " << m_FixedImageNormalizedMin ); - itkDebugMacro( "MovingImageNormalizedMin: " << m_MovingImageNormalizedMin ); - itkDebugMacro( "FixedImageBinSize: " << m_FixedImageBinSize ); - itkDebugMacro( "MovingImageBinSize; " << m_MovingImageBinSize ); - - if( m_UseAllPixels ) { - m_NumberOfSpatialSamples = - this->GetFixedImageRegion().GetNumberOfPixels(); - } - - /** - * Allocate memory for the fixed image sample container. - */ - m_FixedImageSamples.resize( m_NumberOfSpatialSamples ); - - - /** - * Allocate memory for the marginal PDF and initialize values - * to zero. The marginal PDFs are stored as std::vector. - */ - m_FixedImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 ); - m_MovingImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 ); - - /** - * Allocate memory for the joint PDF and joint PDF derivatives. - * The joint PDF and joint PDF derivatives are store as itk::Image. - */ - m_JointPDF = JointPDFType::New(); - - // Instantiate a region, index, size - JointPDFRegionType jointPDFRegion; - JointPDFIndexType jointPDFIndex; - JointPDFSizeType jointPDFSize; - - // Deallocate the memory that may have been allocated for - // previous runs of the metric. - this->m_JointPDFDerivatives = NULL; // by destroying the dynamic array - this->m_PRatioArray.SetSize( 1, 1 ); // and by allocating very small the static ones - this->m_MetricDerivative = DerivativeType( 1 ); - - // - // Now allocate memory according to the user-selected method. - // - if( this->m_UseExplicitPDFDerivatives ) { - this->m_JointPDFDerivatives = JointPDFDerivativesType::New(); - JointPDFDerivativesRegionType jointPDFDerivativesRegion; - JointPDFDerivativesIndexType jointPDFDerivativesIndex; - JointPDFDerivativesSizeType jointPDFDerivativesSize; - - // For the derivatives of the joint PDF define a region starting from {0,0,0} - // with size {m_NumberOfParameters,m_NumberOfHistogramBins, - // m_NumberOfHistogramBins}. The dimension represents transform parameters, - // fixed image parzen window index and moving image parzen window index, - // respectively. - jointPDFDerivativesIndex.Fill( 0 ); - jointPDFDerivativesSize[0] = m_NumberOfParameters; - jointPDFDerivativesSize[1] = m_NumberOfHistogramBins; - jointPDFDerivativesSize[2] = m_NumberOfHistogramBins; - - jointPDFDerivativesRegion.SetIndex( jointPDFDerivativesIndex ); - jointPDFDerivativesRegion.SetSize( jointPDFDerivativesSize ); - - // Set the regions and allocate - m_JointPDFDerivatives->SetRegions( jointPDFDerivativesRegion ); - m_JointPDFDerivatives->Allocate(); - } else { - /** Allocate memory for helper array that will contain the pRatios - * for each bin of the joint histogram. This is part of the effort - * for flattening the computation of the PDF Jacobians. - */ - this->m_PRatioArray.SetSize( this->m_NumberOfHistogramBins, this->m_NumberOfHistogramBins ); - this->m_MetricDerivative = DerivativeType( this->GetNumberOfParameters() ); - } - - // For the joint PDF define a region starting from {0,0} - // with size {m_NumberOfHistogramBins, m_NumberOfHistogramBins}. - // The dimension represents fixed image parzen window index - // and moving image parzen window index, respectively. - jointPDFIndex.Fill( 0 ); - jointPDFSize.Fill( m_NumberOfHistogramBins ); - - jointPDFRegion.SetIndex( jointPDFIndex ); - jointPDFRegion.SetSize( jointPDFSize ); - - // Set the regions and allocate - m_JointPDF->SetRegions( jointPDFRegion ); - m_JointPDF->Allocate(); - - - /** - * Setup the kernels used for the Parzen windows. - */ - m_CubicBSplineKernel = CubicBSplineFunctionType::New(); - m_CubicBSplineDerivativeKernel = CubicBSplineDerivativeFunctionType::New(); - - - if( m_UseAllPixels ) { - /** - * Take all the pixels within the fixed image region) - * to create the sample points list. - */ - this->SampleFullFixedImageDomain( m_FixedImageSamples ); - } else { - /** - * Uniformly sample the fixed image (within the fixed image region) - * to create the sample points list. - */ - this->SampleFixedImageDomain( m_FixedImageSamples ); - } - - /** - * Pre-compute the fixed image parzen window index for - * each point of the fixed image sample points list. - */ - this->ComputeFixedImageParzenWindowIndices( m_FixedImageSamples ); - - /** - * Check if the interpolator is of type BSplineInterpolateImageFunction. - * If so, we can make use of its EvaluateDerivatives method. - * Otherwise, we instantiate an external central difference - * derivative calculator. - * - * TODO: Also add it the possibility of using the default gradient - * provided by the superclass. - * - */ - m_InterpolatorIsBSpline = true; - - BSplineInterpolatorType * testPtr = dynamic_cast( - this->m_Interpolator.GetPointer() ); - if ( !testPtr ) { - m_InterpolatorIsBSpline = false; - - m_DerivativeCalculator = DerivativeFunctionType::New(); - -#ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION - m_DerivativeCalculator->UseImageDirectionOn(); -#endif - - m_DerivativeCalculator->SetInputImage( this->m_MovingImage ); - - m_BSplineInterpolator = NULL; - itkDebugMacro( "Interpolator is not BSpline" ); - } else { - m_BSplineInterpolator = testPtr; - -#ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION - m_BSplineInterpolator->UseImageDirectionOn(); -#endif - - m_DerivativeCalculator = NULL; - itkDebugMacro( "Interpolator is BSpline" ); - } - - /** - * Check if the transform is of type BSplineDeformableTransform. - * - * If so, several speed up features are implemented. - * [1] Precomputing the results of bulk transform for each sample point. - * [2] Precomputing the BSpline weights for each sample point, - * to be used later to directly compute the deformation vector - * [3] Precomputing the indices of the parameters within the - * the support region of each sample point. - */ - m_TransformIsBSpline = true; - - BSplineTransformType * testPtr2 = dynamic_cast( - this->m_Transform.GetPointer() ); - if( !testPtr2 ) { - m_TransformIsBSpline = false; - m_BSplineTransform = NULL; - itkDebugMacro( "Transform is not BSplineDeformable" ); - } else { - m_BSplineTransform = testPtr2; - m_NumParametersPerDim = - m_BSplineTransform->GetNumberOfParametersPerDimension(); - m_NumBSplineWeights = m_BSplineTransform->GetNumberOfWeights(); - itkDebugMacro( "Transform is BSplineDeformable" ); - } - - if( this->m_TransformIsBSpline ) { - // First, deallocate memory that may have been used from previous run of the Metric - this->m_BSplineTransformWeightsArray.SetSize( 1, 1 ); - this->m_BSplineTransformIndicesArray.SetSize( 1, 1 ); - this->m_PreTransformPointsArray.resize( 1 ); - this->m_WithinSupportRegionArray.resize( 1 ); - this->m_Weights.SetSize( 1 ); - this->m_Indices.SetSize( 1 ); - - - if( this->m_UseCachingOfBSplineWeights ) { - m_BSplineTransformWeightsArray.SetSize( - m_NumberOfSpatialSamples, m_NumBSplineWeights ); - m_BSplineTransformIndicesArray.SetSize( - m_NumberOfSpatialSamples, m_NumBSplineWeights ); - m_PreTransformPointsArray.resize( m_NumberOfSpatialSamples ); - m_WithinSupportRegionArray.resize( m_NumberOfSpatialSamples ); - - this->PreComputeTransformValues(); - } else { - this->m_Weights.SetSize( this->m_NumBSplineWeights ); - this->m_Indices.SetSize( this->m_NumBSplineWeights ); - } - - for ( unsigned int j = 0; j < FixedImageDimension; j++ ) { - m_ParametersOffset[j] = j * - m_BSplineTransform->GetNumberOfParametersPerDimension(); - } - } - -} - - -/** - * Uniformly sample the fixed image domain using a random walk - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples ) -{ - - // Set up a random interator within the user specified fixed image region. - typedef ImageRandomConstIteratorWithIndex RandomIterator; - RandomIterator randIter( this->m_FixedImage, this->GetFixedImageRegion() ); - - randIter.SetNumberOfSamples( m_NumberOfSpatialSamples ); - randIter.GoToBegin(); - - typename FixedImageSpatialSampleContainer::iterator iter; - typename FixedImageSpatialSampleContainer::const_iterator end=samples.end(); - - if( this->m_FixedImageMask ) { - - InputPointType inputPoint; - - iter=samples.begin(); - int count = 0; - int samples_found = 0; - int maxcount = m_NumberOfSpatialSamples * 10; - while( iter != end ) { - - if ( count > maxcount ) { -#if 0 - itkExceptionMacro( - "Drew too many samples from the mask (is it too small?): " - << maxcount << std::endl ); -#else -samples.resize(samples_found); -// this->SetNumberOfSpatialSamples(sample_found); -break; -#endif - } - count++; - - // Get sampled index - FixedImageIndexType index = randIter.GetIndex(); - // Check if the Index is inside the mask, translate index to point - this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - // If not inside the mask, ignore the point - if( !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++randIter; // jump to another random position - continue; - } - - // Get sampled fixed image value - (*iter).FixedImageValue = randIter.Get(); - // Translate index to point - (*iter).FixedImagePointValue = inputPoint; - samples_found++; - // Jump to random position - ++randIter; - ++iter; - } - } else { - for( iter=samples.begin(); iter != end; ++iter ) { - // Get sampled index - FixedImageIndexType index = randIter.GetIndex(); - // Get sampled fixed image value - (*iter).FixedImageValue = randIter.Get(); - // Translate index to point - this->m_FixedImage->TransformIndexToPhysicalPoint( index, - (*iter).FixedImagePointValue ); - // Jump to random position - ++randIter; - - } - } -} - -/** - * Sample the fixed image domain using all pixels in the Fixed image region - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::SampleFullFixedImageDomain( FixedImageSpatialSampleContainer& samples ) -{ - - // Set up a region interator within the user specified fixed image region. - typedef ImageRegionConstIteratorWithIndex RegionIterator; - RegionIterator regionIter( this->m_FixedImage, this->GetFixedImageRegion() ); - - regionIter.GoToBegin(); - - typename FixedImageSpatialSampleContainer::iterator iter; - typename FixedImageSpatialSampleContainer::const_iterator end=samples.end(); - - if( this->m_FixedImageMask ) { - InputPointType inputPoint; - - iter=samples.begin(); - unsigned long nSamplesPicked = 0; - - while( iter != end && !regionIter.IsAtEnd() ) { - // Get sampled index - FixedImageIndexType index = regionIter.GetIndex(); - // Check if the Index is inside the mask, translate index to point - this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); - - // If not inside the mask, ignore the point - if( !this->m_FixedImageMask->IsInside( inputPoint ) ) { - ++regionIter; // jump to next pixel - continue; - } - - // Get sampled fixed image value - (*iter).FixedImageValue = regionIter.Get(); - // Translate index to point - (*iter).FixedImagePointValue = inputPoint; - - ++regionIter; - ++iter; - ++nSamplesPicked; - } - - // If we picked fewer samples than the desired number, - // resize the container - if (nSamplesPicked != this->m_NumberOfSpatialSamples) { - this->m_NumberOfSpatialSamples = nSamplesPicked; - samples.resize(this->m_NumberOfSpatialSamples); - } - } else { // not restricting sample throwing to a mask - - // cannot sample more than the number of pixels in the image region - if ( this->m_NumberOfSpatialSamples - > this->GetFixedImageRegion().GetNumberOfPixels()) { - this->m_NumberOfSpatialSamples - = this->GetFixedImageRegion().GetNumberOfPixels(); - samples.resize(this->m_NumberOfSpatialSamples); - } - - for( iter=samples.begin(); iter != end; ++iter ) { - // Get sampled index - FixedImageIndexType index = regionIter.GetIndex(); - // Get sampled fixed image value - (*iter).FixedImageValue = regionIter.Get(); - // Translate index to point - this->m_FixedImage->TransformIndexToPhysicalPoint( index, - (*iter).FixedImagePointValue ); - ++regionIter; - } - } -} - -/** - * Uniformly sample the fixed image domain using a random walk - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::ComputeFixedImageParzenWindowIndices( - FixedImageSpatialSampleContainer& samples ) -{ - - typename FixedImageSpatialSampleContainer::iterator iter; - typename FixedImageSpatialSampleContainer::const_iterator end=samples.end(); - - for( iter=samples.begin(); iter != end; ++iter ) { - - // Determine parzen window arguments (see eqn 6 of Mattes paper [2]). - double windowTerm = - static_cast( (*iter).FixedImageValue ) / m_FixedImageBinSize - - m_FixedImageNormalizedMin; - unsigned int pindex = static_cast( vcl_floor(windowTerm ) ); - - // Make sure the extreme values are in valid bins - if ( pindex < 2 ) { - pindex = 2; - } else if ( pindex > (m_NumberOfHistogramBins - 3) ) { - pindex = m_NumberOfHistogramBins - 3; - } - - (*iter).FixedImageParzenWindowIndex = pindex; - - } - -} - -/** - * Get the match Measure - */ -template < class TFixedImage, class TMovingImage > -typename MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::MeasureType -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::GetValue( const ParametersType& parameters ) const -{ - - // Reset marginal pdf to all zeros. - // Assumed the size has already been set to NumberOfHistogramBins - // in Initialize(). - for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) { - m_FixedImageMarginalPDF[j] = 0.0; - m_MovingImageMarginalPDF[j] = 0.0; - } - - // Reset the joint pdfs to zero - m_JointPDF->FillBuffer( 0.0 ); - - // Set up the parameters in the transform - this->m_Transform->SetParameters( parameters ); - - - // Declare iterators for iteration over the sample container - typename FixedImageSpatialSampleContainer::const_iterator fiter; - typename FixedImageSpatialSampleContainer::const_iterator fend = - m_FixedImageSamples.end(); - - unsigned long nSamples=0; - unsigned long nFixedImageSamples=0; - - - for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) { - - // Get moving image value - MovingImagePointType mappedPoint; - bool sampleOk; - double movingImageValue; - - this->TransformPoint( nFixedImageSamples, parameters, mappedPoint, - sampleOk, movingImageValue ); - - ++nFixedImageSamples; - - if( sampleOk ) { - - ++nSamples; - - /** - * Compute this sample's contribution to the marginal and - * joint distributions. - * - */ - - // Determine parzen window arguments (see eqn 6 of Mattes paper [2]) - double movingImageParzenWindowTerm = - movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin; - unsigned int movingImageParzenWindowIndex = - static_cast( vcl_floor(movingImageParzenWindowTerm ) ); - - // Make sure the extreme values are in valid bins - if ( movingImageParzenWindowIndex < 2 ) { - movingImageParzenWindowIndex = 2; - } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) { - movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3; - } - - - // Since a zero-order BSpline (box car) kernel is used for - // the fixed image marginal pdf, we need only increment the - // fixedImageParzenWindowIndex by value of 1.0. - m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] += - static_cast( 1 ); - - /** - * The region of support of the parzen window determines which bins - * of the joint PDF are effected by the pair of image values. - * Since we are using a cubic spline for the moving image parzen - * window, four bins are affected. The fixed image parzen window is - * a zero-order spline (box car) and thus effects only one bin. - * - * The PDF is arranged so that moving image bins corresponds to the - * zero-th (column) dimension and the fixed image bins corresponds - * to the first (row) dimension. - * - */ - - // Pointer to affected bin to be updated - JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() + - ( (*fiter).FixedImageParzenWindowIndex - * m_JointPDF->GetOffsetTable()[1] ); - - // Move the pointer to the first affected bin - int pdfMovingIndex = static_cast( movingImageParzenWindowIndex ) - 1; - pdfPtr += pdfMovingIndex; - - for (; pdfMovingIndex <= static_cast( movingImageParzenWindowIndex ) - + 2; - pdfMovingIndex++, pdfPtr++ ) { - - // Update PDF for the current intensity pair - double movingImageParzenWindowArg = - static_cast( pdfMovingIndex ) - - static_cast( movingImageParzenWindowTerm ); - - *(pdfPtr) += static_cast( - m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) ); - - } //end parzen windowing for loop - - } //end if-block check sampleOk - } // end iterating over fixed image spatial sample container for loop - - itkDebugMacro( "Ratio of voxels mapping into moving image buffer: " - << nSamples << " / " << m_NumberOfSpatialSamples - << std::endl ); - - if( nSamples < m_NumberOfSpatialSamples / 16 ) { - itkExceptionMacro( "Too many samples map outside moving image buffer: " - << nSamples << " / " << m_NumberOfSpatialSamples - << std::endl ); - } - - this->m_NumberOfPixelsCounted = nSamples; - - - /** - * Normalize the PDFs, compute moving image marginal PDF - * - */ - typedef ImageRegionIterator JointPDFIteratorType; - JointPDFIteratorType jointPDFIterator ( m_JointPDF, - m_JointPDF->GetBufferedRegion() ); - - jointPDFIterator.GoToBegin(); - - // Compute joint PDF normalization factor - // (to ensure joint PDF sum adds to 1.0) - double jointPDFSum = 0.0; - - while( !jointPDFIterator.IsAtEnd() ) { - jointPDFSum += jointPDFIterator.Get(); - ++jointPDFIterator; - } - - if ( jointPDFSum == 0.0 ) { - itkExceptionMacro( "Joint PDF summed to zero" ); - } - - - // Normalize the PDF bins - jointPDFIterator.GoToEnd(); - while( !jointPDFIterator.IsAtBegin() ) { - --jointPDFIterator; - jointPDFIterator.Value() /= static_cast( jointPDFSum ); - } - - - // Normalize the fixed image marginal PDF - double fixedPDFSum = 0.0; - for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) { - fixedPDFSum += m_FixedImageMarginalPDF[bin]; - } - - if ( fixedPDFSum == 0.0 ) { - itkExceptionMacro( "Fixed image marginal PDF summed to zero" ); - } - - for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) { - m_FixedImageMarginalPDF[bin] /= static_cast( fixedPDFSum ); - } - - - // Compute moving image marginal PDF by summing over fixed image bins. - typedef ImageLinearIteratorWithIndex JointPDFLinearIterator; - JointPDFLinearIterator linearIter( m_JointPDF, - m_JointPDF->GetBufferedRegion() ); - - linearIter.SetDirection( 1 ); - linearIter.GoToBegin(); - unsigned int movingIndex1 = 0; - - while( !linearIter.IsAtEnd() ) { - - double sum = 0.0; - - while( !linearIter.IsAtEndOfLine() ) { - sum += linearIter.Get(); - ++linearIter; - } - - m_MovingImageMarginalPDF[movingIndex1] = static_cast(sum); - - linearIter.NextLine(); - ++movingIndex1; - - } - - /** - * Compute the metric by double summation over histogram. - */ - - // Setup pointer to point to the first bin - JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer(); - - // Initialze sum to zero - double sum = 0.0; - - for( unsigned int fixedIndex = 0; - fixedIndex < m_NumberOfHistogramBins; - ++fixedIndex ) { - - double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex]; - - for( unsigned int movingIndex = 0; - movingIndex < m_NumberOfHistogramBins; - ++movingIndex, jointPDFPtr++ ) { - - double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex]; - double jointPDFValue = *(jointPDFPtr); - - // check for non-zero bin contribution - if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) { - - double pRatio = vcl_log(jointPDFValue / movingImagePDFValue ); - if( fixedImagePDFValue > 1e-16) { - sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) ); - } - - } // end if-block to check non-zero bin contribution - } // end for-loop over moving index - } // end for-loop over fixed index - - return static_cast( -1.0 * sum ); - -} - - -/** - * Get the both Value and Derivative Measure - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::GetValueAndDerivative( - const ParametersType& parameters, - MeasureType& value, - DerivativeType& derivative) const -{ - - // Set output values to zero - value = NumericTraits< MeasureType >::Zero; - - if( this->m_UseExplicitPDFDerivatives ) { - m_JointPDFDerivatives->FillBuffer( 0.0 ); - derivative = DerivativeType( this->GetNumberOfParameters() ); - derivative.Fill( NumericTraits< MeasureType >::Zero ); - } else { - this->m_MetricDerivative.Fill( NumericTraits< MeasureType >::Zero ); - this->m_PRatioArray.Fill( 0.0 ); - } - - // Reset marginal pdf to all zeros. - // Assumed the size has already been set to NumberOfHistogramBins - // in Initialize(). - for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) { - m_FixedImageMarginalPDF[j] = 0.0; - m_MovingImageMarginalPDF[j] = 0.0; - } - - // Reset the joint pdfs to zero - m_JointPDF->FillBuffer( 0.0 ); - - - // Set up the parameters in the transform - this->m_Transform->SetParameters( parameters ); - - - // Declare iterators for iteration over the sample container - typename FixedImageSpatialSampleContainer::const_iterator fiter; - typename FixedImageSpatialSampleContainer::const_iterator fend = - m_FixedImageSamples.end(); - - unsigned long nSamples=0; - unsigned long nFixedImageSamples=0; - - for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) { - - // Get moving image value - MovingImagePointType mappedPoint; - bool sampleOk; - double movingImageValue; - - - this->TransformPoint( nFixedImageSamples, parameters, mappedPoint, - sampleOk, movingImageValue ); - - if( sampleOk ) { - ++nSamples; - - // Get moving image derivative at the mapped position - ImageDerivativesType movingImageGradientValue; - this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue ); - - - /** - * Compute this sample's contribution to the marginal - * and joint distributions. - * - */ - - // Determine parzen window arguments (see eqn 6 of Mattes paper [2]) - double movingImageParzenWindowTerm = - movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin; - unsigned int movingImageParzenWindowIndex = - static_cast( vcl_floor(movingImageParzenWindowTerm ) ); - - // Make sure the extreme values are in valid bins - if ( movingImageParzenWindowIndex < 2 ) { - movingImageParzenWindowIndex = 2; - } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) { - movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3; - } - - - // Since a zero-order BSpline (box car) kernel is used for - // the fixed image marginal pdf, we need only increment the - // fixedImageParzenWindowIndex by value of 1.0. - m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] += - static_cast( 1 ); - - /** - * The region of support of the parzen window determines which bins - * of the joint PDF are effected by the pair of image values. - * Since we are using a cubic spline for the moving image parzen - * window, four bins are effected. The fixed image parzen window is - * a zero-order spline (box car) and thus effects only one bin. - * - * The PDF is arranged so that moving image bins corresponds to the - * zero-th (column) dimension and the fixed image bins corresponds - * to the first (row) dimension. - * - */ - - // Pointer to affected bin to be updated - JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() + - ( (*fiter).FixedImageParzenWindowIndex * m_NumberOfHistogramBins ); - - // Move the pointer to the fist affected bin - int pdfMovingIndex = static_cast( movingImageParzenWindowIndex ) - 1; - pdfPtr += pdfMovingIndex; - - for (; pdfMovingIndex <= static_cast( movingImageParzenWindowIndex ) - + 2; - pdfMovingIndex++, pdfPtr++ ) { - // Update PDF for the current intensity pair - double movingImageParzenWindowArg = - static_cast( pdfMovingIndex ) - - static_cast(movingImageParzenWindowTerm); - - *(pdfPtr) += static_cast( - m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) ); - - if( this->m_UseExplicitPDFDerivatives ) { - // Compute the cubicBSplineDerivative for later repeated use. - double cubicBSplineDerivativeValue = - m_CubicBSplineDerivativeKernel->Evaluate( - movingImageParzenWindowArg ); - - // Compute PDF derivative contribution. - this->ComputePDFDerivatives( nFixedImageSamples, - pdfMovingIndex, - movingImageGradientValue, - cubicBSplineDerivativeValue ); - } - - } //end parzen windowing for loop - - } //end if-block check sampleOk - - ++nFixedImageSamples; - - } // end iterating over fixed image spatial sample container for loop - - itkDebugMacro( "Ratio of voxels mapping into moving image buffer: " - << nSamples << " / " << m_NumberOfSpatialSamples - << std::endl ); - - if( nSamples < m_NumberOfSpatialSamples / 16 ) { - itkExceptionMacro( "Too many samples map outside moving image buffer: " - << nSamples << " / " << m_NumberOfSpatialSamples - << std::endl ); - } - - this->m_NumberOfPixelsCounted = nSamples; - - /** - * Normalize the PDFs, compute moving image marginal PDF - * - */ - typedef ImageRegionIterator JointPDFIteratorType; - JointPDFIteratorType jointPDFIterator ( m_JointPDF, - m_JointPDF->GetBufferedRegion() ); - - jointPDFIterator.GoToBegin(); - - - // Compute joint PDF normalization factor - // (to ensure joint PDF sum adds to 1.0) - double jointPDFSum = 0.0; - - while( !jointPDFIterator.IsAtEnd() ) { - jointPDFSum += jointPDFIterator.Get(); - ++jointPDFIterator; - } - - if ( jointPDFSum == 0.0 ) { - itkExceptionMacro( "Joint PDF summed to zero" ); - } - - - // Normalize the PDF bins - jointPDFIterator.GoToEnd(); - while( !jointPDFIterator.IsAtBegin() ) { - --jointPDFIterator; - jointPDFIterator.Value() /= static_cast( jointPDFSum ); - } - - - // Normalize the fixed image marginal PDF - double fixedPDFSum = 0.0; - for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) { - fixedPDFSum += m_FixedImageMarginalPDF[bin]; - } - - if ( fixedPDFSum == 0.0 ) { - itkExceptionMacro( "Fixed image marginal PDF summed to zero" ); - } - - for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) { - m_FixedImageMarginalPDF[bin] /= static_cast( fixedPDFSum ); - } - - - // Compute moving image marginal PDF by summing over fixed image bins. - typedef ImageLinearIteratorWithIndex JointPDFLinearIterator; - JointPDFLinearIterator linearIter( - m_JointPDF, m_JointPDF->GetBufferedRegion() ); - - linearIter.SetDirection( 1 ); - linearIter.GoToBegin(); - unsigned int movingIndex1 = 0; - - while( !linearIter.IsAtEnd() ) { - - double sum = 0.0; - - while( !linearIter.IsAtEndOfLine() ) { - sum += linearIter.Get(); - ++linearIter; - } - - m_MovingImageMarginalPDF[movingIndex1] = static_cast(sum); - - linearIter.NextLine(); - ++movingIndex1; - - } - - double nFactor = 1.0 / ( m_MovingImageBinSize - * static_cast( nSamples ) ); - - if( this->m_UseExplicitPDFDerivatives ) { - // Normalize the joint PDF derivatives by the test image binsize and nSamples - typedef ImageRegionIterator - JointPDFDerivativesIteratorType; - JointPDFDerivativesIteratorType jointPDFDerivativesIterator ( - m_JointPDFDerivatives, - m_JointPDFDerivatives->GetBufferedRegion() - ); - - jointPDFDerivativesIterator.GoToBegin(); - - while( !jointPDFDerivativesIterator.IsAtEnd() ) { - jointPDFDerivativesIterator.Value() *= nFactor; - ++jointPDFDerivativesIterator; - } - } - - /** - * Compute the metric by double summation over histogram. - */ - - // Setup pointer to point to the first bin - JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer(); - - // Initialize sum to zero - double sum = 0.0; - - for( unsigned int fixedIndex = 0; - fixedIndex < m_NumberOfHistogramBins; - ++fixedIndex ) { - double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex]; - - for( unsigned int movingIndex = 0; movingIndex < m_NumberOfHistogramBins; - ++movingIndex, jointPDFPtr++ ) { - double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex]; - double jointPDFValue = *(jointPDFPtr); - - // check for non-zero bin contribution - if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) { - - double pRatio = vcl_log(jointPDFValue / movingImagePDFValue ); - - if( fixedImagePDFValue > 1e-16) { - sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) ); - } - - if( this->m_UseExplicitPDFDerivatives ) { - // move joint pdf derivative pointer to the right position - JointPDFValueType * derivPtr = m_JointPDFDerivatives->GetBufferPointer() - + ( fixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] ) - + ( movingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] ); - - for( unsigned int parameter=0; parameter < m_NumberOfParameters; ++parameter, derivPtr++ ) { - // Ref: eqn 23 of Thevenaz & Unser paper [3] - derivative[parameter] -= (*derivPtr) * pRatio; - } // end for-loop over parameters - } else { - this->m_PRatioArray[fixedIndex][movingIndex] = pRatio * nFactor; - } - } // end if-block to check non-zero bin contribution - } // end for-loop over moving index - } // end for-loop over fixed index - - if( !(this->m_UseExplicitPDFDerivatives ) ) { - // Second pass: This one is done for accumulating the contributions - // to the derivative array. - - nFixedImageSamples = 0; - - for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) { - - // Get moving image value - MovingImagePointType mappedPoint; - bool sampleOk; - double movingImageValue; - - this->TransformPoint( nFixedImageSamples, parameters, mappedPoint, - sampleOk, movingImageValue ); - - if( sampleOk ) { - // Get moving image derivative at the mapped position - ImageDerivativesType movingImageGradientValue; - this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue ); - - - /** - * Compute this sample's contribution to the marginal - * and joint distributions. - * - */ - - // Determine parzen window arguments (see eqn 6 of Mattes paper [2]). - double movingImageParzenWindowTerm = - movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin; - unsigned int movingImageParzenWindowIndex = - static_cast( vcl_floor(movingImageParzenWindowTerm ) ); - - // Make sure the extreme values are in valid bins - if ( movingImageParzenWindowIndex < 2 ) { - movingImageParzenWindowIndex = 2; - } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) { - movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3; - } - - - // Move the pointer to the fist affected bin - int pdfMovingIndex = static_cast( movingImageParzenWindowIndex ) - 1; - - for (; pdfMovingIndex <= static_cast( movingImageParzenWindowIndex ) - + 2; - pdfMovingIndex++ ) { - - // Update PDF for the current intensity pair - double movingImageParzenWindowArg = - static_cast( pdfMovingIndex ) - - static_cast(movingImageParzenWindowTerm); - - // Compute the cubicBSplineDerivative for later repeated use. - double cubicBSplineDerivativeValue = - m_CubicBSplineDerivativeKernel->Evaluate( - movingImageParzenWindowArg ); - - // Compute PDF derivative contribution. - this->ComputePDFDerivatives( nFixedImageSamples, - pdfMovingIndex, - movingImageGradientValue, - cubicBSplineDerivativeValue ); - - - } //end parzen windowing for loop - - } //end if-block check sampleOk - - ++nFixedImageSamples; - - } // end iterating over fixed image spatial sample container for loop - - - derivative = this->m_MetricDerivative; - } - - value = static_cast( -1.0 * sum ); -} - - -/** - * Get the match measure derivative - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::GetDerivative( const ParametersType& parameters, - DerivativeType & derivative ) const -{ - MeasureType value; - // call the combined version - this->GetValueAndDerivative( parameters, value, derivative ); -} - - -/** - * Compute image derivatives using a central difference function - * if we are not using a BSplineInterpolator, which includes - * derivatives. - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::ComputeImageDerivatives( - const MovingImagePointType& mappedPoint, - ImageDerivativesType& gradient ) const -{ - - if( m_InterpolatorIsBSpline ) { - // Computed moving image gradient using derivative BSpline kernel. - gradient = m_BSplineInterpolator->EvaluateDerivative( mappedPoint ); - } else { - // For all generic interpolator use central differencing. - gradient = m_DerivativeCalculator->Evaluate( mappedPoint ); - } - -} - - -/** - * Transform a point from FixedImage domain to MovingImage domain. - * This function also checks if mapped point is within support region. - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::TransformPoint( - unsigned int sampleNumber, - const ParametersType& parameters, - MovingImagePointType& mappedPoint, - bool& sampleOk, - double& movingImageValue ) const -{ - - if ( !m_TransformIsBSpline ) { - // Use generic transform to compute mapped position - mappedPoint = this->m_Transform->TransformPoint( - m_FixedImageSamples[sampleNumber].FixedImagePointValue ); - - // Check if mapped point inside image buffer - sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint ); - } else { - - if( this->m_UseCachingOfBSplineWeights ) { - // If the transform is BSplineDeformable, we can use the precomputed - // weights and indices to obtained the mapped position - // - const WeightsValueType * weights = - m_BSplineTransformWeightsArray[sampleNumber]; - const IndexValueType * indices = - m_BSplineTransformIndicesArray[sampleNumber]; - mappedPoint.Fill( 0.0 ); - - if ( m_WithinSupportRegionArray[sampleNumber] ) { - for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ ) { - for ( unsigned int j = 0; j < FixedImageDimension; j++ ) { - mappedPoint[j] += weights[k] * - parameters[ indices[k] + m_ParametersOffset[j] ]; - } - } - } - - for( unsigned int j = 0; j < FixedImageDimension; j++ ) { - mappedPoint[j] += m_PreTransformPointsArray[sampleNumber][j]; - } - - // Check if mapped point inside image buffer - sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint ); - - // Check if mapped point is within the support region of a grid point. - // This is neccessary for computing the metric gradient - sampleOk = sampleOk && m_WithinSupportRegionArray[sampleNumber]; - } else { - // If not caching values, we invoke the Transform to recompute the - // mapping of the point. - this->m_BSplineTransform->TransformPoint( - this->m_FixedImageSamples[sampleNumber].FixedImagePointValue, - mappedPoint, this->m_Weights, this->m_Indices, sampleOk); - - // Check if mapped point inside image buffer - sampleOk = sampleOk && this->m_Interpolator->IsInsideBuffer( mappedPoint ); - } - - } - - // If user provided a mask over the Moving image - if ( this->m_MovingImageMask ) { - // Check if mapped point is within the support region of the moving image - // mask - sampleOk = sampleOk && this->m_MovingImageMask->IsInside( mappedPoint ); - } - - - if ( sampleOk ) { - movingImageValue = this->m_Interpolator->Evaluate( mappedPoint ); - - if ( movingImageValue < m_MovingImageTrueMin || - movingImageValue > m_MovingImageTrueMax ) { - // need to throw out this sample as it will not fall into a valid bin - sampleOk = false; - } - } -} - - -/** - * Compute PDF derivatives contribution for each parameter - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::ComputePDFDerivatives( - unsigned int sampleNumber, - int pdfMovingIndex, - const ImageDerivativesType& movingImageGradientValue, - double cubicBSplineDerivativeValue ) const -{ - - const int pdfFixedIndex = - m_FixedImageSamples[sampleNumber].FixedImageParzenWindowIndex; - - JointPDFValueType * derivPtr = NULL; - double precomputedWeight = 0.0; - - if( this->m_UseExplicitPDFDerivatives ) { - // Update bins in the PDF derivatives for the current intensity pair - derivPtr = m_JointPDFDerivatives->GetBufferPointer() + - ( pdfFixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] ) + - ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] ); - } else { - // Recover the precomputed weight for this specific PDF bin - precomputedWeight = this->m_PRatioArray[pdfFixedIndex][pdfMovingIndex]; - } - - if( !m_TransformIsBSpline ) { - - /** - * Generic version which works for all transforms. - */ - - // Compute the transform Jacobian. - typedef typename TransformType::JacobianType JacobianType; -#if ITK_VERSION_MAJOR >= 4 - JacobianType jacobian; - this->m_Transform->ComputeJacobianWithRespectToParameters( m_FixedImageSamples[sampleNumber].FixedImagePointValue, jacobian ); -#else - const JacobianType & jacobian = - this->m_Transform->GetJacobian( m_FixedImageSamples[sampleNumber].FixedImagePointValue ); -#endif - - for ( unsigned int mu = 0; mu < m_NumberOfParameters; mu++ ) { - double innerProduct = 0.0; - for ( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) { - innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim]; - } - - const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue; - - if( this->m_UseExplicitPDFDerivatives ) { - *(derivPtr) -= derivativeContribution; - ++derivPtr; - } else { - this->m_MetricDerivative[mu] += precomputedWeight * derivativeContribution; - } - } - - } else { - const WeightsValueType * weights = NULL; - const IndexValueType * indices = NULL; - - if( this->m_UseCachingOfBSplineWeights ) { - /** - * If the transform is of type BSplineDeformableTransform, - * we can obtain a speed up by only processing the affected parameters. - * Note that these pointers are just pointing to pre-allocated rows - * of the caching arrays. There is therefore, no need to free this - * memory. - */ - weights = m_BSplineTransformWeightsArray[sampleNumber]; - indices = m_BSplineTransformIndicesArray[sampleNumber]; - } else { -#if ITK_VERSION_MAJOR >= 4 - m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition( - m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices ); -#else - m_BSplineTransform->GetJacobian( - m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices ); -#endif - } - - for( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) { - - double innerProduct; - int parameterIndex; - - for( unsigned int mu = 0; mu < m_NumBSplineWeights; mu++ ) { - - /* The array weights contains the Jacobian values in a 1-D array - * (because for each parameter the Jacobian is non-zero in only 1 of the - * possible dimensions) which is multiplied by the moving image - * gradient. */ - if( this->m_UseCachingOfBSplineWeights ) { - innerProduct = movingImageGradientValue[dim] * weights[mu]; - parameterIndex = indices[mu] + m_ParametersOffset[dim]; - } else { - innerProduct = movingImageGradientValue[dim] * this->m_Weights[mu]; - parameterIndex = this->m_Indices[mu] + this->m_ParametersOffset[dim]; - } - - const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue; - - if( this->m_UseExplicitPDFDerivatives ) { - JointPDFValueType * ptr = derivPtr + parameterIndex; - *(ptr) -= derivativeContribution; - } else { - this->m_MetricDerivative[parameterIndex] += precomputedWeight * derivativeContribution; - } - - } //end mu for loop - } //end dim for loop - - } // end if-block transform is BSpline - -} - - -// Method to reinitialize the seed of the random number generator -template < class TFixedImage, class TMovingImage > void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::ReinitializeSeed() -{ - Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed(); -} - -// Method to reinitialize the seed of the random number generator -template < class TFixedImage, class TMovingImage > void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::ReinitializeSeed(int seed) -{ - Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed( - seed); -} - - -/** - * Cache pre-transformed points, weights and indices. - * This method is only called if the flag UseCachingOfBSplineWeights is ON. - */ -template < class TFixedImage, class TMovingImage > -void -MattesMutualInformationImageToImageMetricFor3DBLUTFFD -::PreComputeTransformValues() -{ - // Create all zero dummy transform parameters - ParametersType dummyParameters( this->m_Transform->GetNumberOfParameters() ); - dummyParameters.Fill( 0.0 ); - this->m_Transform->SetParameters( dummyParameters ); - - // Cycle through each sampled fixed image point - BSplineTransformWeightsType weights( m_NumBSplineWeights ); - BSplineTransformIndexArrayType indices( m_NumBSplineWeights ); - bool valid; - MovingImagePointType mappedPoint; - - // Declare iterators for iteration over the sample container - typename FixedImageSpatialSampleContainer::const_iterator fiter; - typename FixedImageSpatialSampleContainer::const_iterator fend = - m_FixedImageSamples.end(); - unsigned long counter = 0; - - for( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter, counter++ ) { - m_BSplineTransform->TransformPoint( - m_FixedImageSamples[counter].FixedImagePointValue, - mappedPoint, weights, indices, valid ); - - for( unsigned long k = 0; k < m_NumBSplineWeights; k++ ) { - m_BSplineTransformWeightsArray[counter][k] = weights[k]; - m_BSplineTransformIndicesArray[counter][k] = indices[k]; - } - - m_PreTransformPointsArray[counter] = mappedPoint; - m_WithinSupportRegionArray[counter] = valid; - - } - -} - - -} // end namespace itk - - -#endif - #endif diff --git a/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h b/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h index f6d745d..c3be983 100644 --- a/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h +++ b/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h @@ -40,112 +40,6 @@ // gets integrated into the main directories. #include "itkConfigure.h" -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.h" -#else - -#include "itkImageToImageMetric.h" -#include "itkCovariantVector.h" -#include "itkPoint.h" - - -namespace itk -{ -/** \class MeanSquaresImageToImageMetricFor3DBLUTFFD - * \brief Computes similarity between two objects to be registered - * - * This Class is templated over the type of the fixed and moving - * images to be compared. - * - * This metric computes the sum of squared differences between pixels in - * the moving image and pixels in the fixed image. The spatial correspondance - * between both images is established through a Transform. Pixel values are - * taken from the Moving image. Their positions are mapped to the Fixed image - * and result in general in non-grid position on it. Values at these non-grid - * position of the Fixed image are interpolated using a user-selected Interpolator. - * - * \ingroup RegistrationMetrics - */ -template < class TFixedImage, class TMovingImage > -class ITK_EXPORT MeanSquaresImageToImageMetricFor3DBLUTFFD : - public ImageToImageMetric< TFixedImage, TMovingImage> -{ -public: - - /** Standard class typedefs. */ - typedef MeanSquaresImageToImageMetricFor3DBLUTFFD Self; - typedef ImageToImageMetric Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(MeanSquaresImageToImageMetricFor3DBLUTFFD, ImageToImageMetric); - - - /** Types transferred from the base class */ - typedef typename Superclass::RealType RealType; - typedef typename Superclass::TransformType TransformType; - typedef typename Superclass::TransformPointer TransformPointer; - typedef typename Superclass::TransformParametersType TransformParametersType; - typedef typename Superclass::TransformJacobianType TransformJacobianType; - typedef typename Superclass::GradientPixelType GradientPixelType; - typedef typename Superclass::GradientImageType GradientImageType; - typedef typename Superclass::InputPointType InputPointType; - typedef typename Superclass::OutputPointType OutputPointType; - - typedef typename Superclass::MeasureType MeasureType; - typedef typename Superclass::DerivativeType DerivativeType; - typedef typename Superclass::FixedImageType FixedImageType; - typedef typename Superclass::MovingImageType MovingImageType; - typedef typename Superclass::FixedImageConstPointer FixedImageConstPointer; - typedef typename Superclass::MovingImageConstPointer MovingImageConstPointer; - - - /** Get the derivatives of the match measure. */ - void GetDerivative( const TransformParametersType & parameters, - DerivativeType & derivative ) const; - - /** Get the value for single valued optimizers. */ - MeasureType GetValue( const TransformParametersType & parameters ) const; - - /** Get value and derivatives for multiple valued optimizers. */ - void GetValueAndDerivative( const TransformParametersType & parameters, - MeasureType& Value, DerivativeType& Derivative ) const; - -#ifdef ITK_USE_CONCEPT_CHECKING - /** Begin concept checking */ - itkConceptMacro(MovingPixelTypeHasNumericTraitsCheck, - (Concept::HasNumericTraits)); - itkConceptMacro(MovingRealTypeAdditivieOperatorsCheck, - (Concept::AdditiveOperators< - typename NumericTraits::RealType>)); - itkConceptMacro(MovingRealTypeMultiplyOperatorCheck, - (Concept::MultiplyOperator< - typename NumericTraits::RealType>)); - - /** End concept checking */ -#endif - -protected: - MeanSquaresImageToImageMetricFor3DBLUTFFD(); - virtual ~MeanSquaresImageToImageMetricFor3DBLUTFFD() {}; - -private: - MeanSquaresImageToImageMetricFor3DBLUTFFD(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented - -}; - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx" -#endif - -#endif - #endif diff --git a/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx b/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx index d6793e5..67e4ac4 100644 --- a/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx +++ b/registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx @@ -40,9 +40,7 @@ // gets integrated into the main directories. #include "itkConfigure.h" -#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx" -#else #include "itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h" #include "itkImageRegionConstIteratorWithIndex.h" @@ -197,13 +195,8 @@ MeanSquaresImageToImageMetricFor3DBLUTFFD if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) { const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint ); -#if ITK_VERSION_MAJOR >= 4 TransformJacobianType jacobian; this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian ); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif const RealType fixedValue = ti.Value(); this->m_NumberOfPixelsCounted++; @@ -315,13 +308,8 @@ MeanSquaresImageToImageMetricFor3DBLUTFFD if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) { const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint ); -#if ITK_VERSION_MAJOR >= 4 TransformJacobianType jacobian; this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian ); -#else - const TransformJacobianType & jacobian = - this->m_Transform->GetJacobian( inputPoint ); -#endif const RealType fixedValue = ti.Value(); this->m_NumberOfPixelsCounted++; diff --git a/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h b/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h index 7cfe51e..1eaa49a 100644 --- a/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h +++ b/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h @@ -35,11 +35,7 @@ #ifndef __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_h #define __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_h -#if ITK_VERSION_MAJOR >= 4 - #include "itkImageToImageMetric.h" -#else - #include "itkOptImageToImageMetric.h" -#endif +#include "itkImageToImageMetric.h" #include "itkCovariantVector.h" #include "itkPoint.h" #include "itkIndex.h" diff --git a/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx b/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx index f55b19a..d1b5b4c 100644 --- a/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx +++ b/registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx @@ -630,9 +630,6 @@ MattesMutualInformationImageToImageMetricFor3DBLUTFFD { // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // MUST BE CALLED TO INITIATE PROCESSING this->GetValueMultiThreadedInitiate(); @@ -941,9 +938,6 @@ MattesMutualInformationImageToImageMetricFor3DBLUTFFD // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES this->GetValueAndDerivativeMultiThreadedInitiate(); @@ -1156,13 +1150,8 @@ MattesMutualInformationImageToImageMetricFor3DBLUTFFD transform = this->m_Transform; } -#if ITK_VERSION_MAJOR >= 4 JacobianType jacobian; transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[sampleNumber].point, jacobian); -#else - const JacobianType& jacobian = - transform->GetJacobian( this->m_FixedImageSamples[sampleNumber].point ); -#endif // for ( unsigned int mu = 0; mu < this->m_NumberOfParameters; mu++ ) // { @@ -1230,15 +1219,9 @@ MattesMutualInformationImageToImageMetricFor3DBLUTFFD indicesHelper = &(this->m_BSplineTransformIndices); } -#if ITK_VERSION_MAJOR >= 4 this->m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition( this->m_FixedImageSamples[sampleNumber].point, *weightsHelper, *indicesHelper ); -#else - this->m_BSplineTransform->GetJacobian( - this->m_FixedImageSamples[sampleNumber].point, - *weightsHelper, *indicesHelper ); -#endif } for( unsigned int dim = 0; dim < Superclass::FixedImageDimension; dim++ ) { diff --git a/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.h b/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.h index acb0b1b..949246d 100644 --- a/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.h +++ b/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.h @@ -35,11 +35,7 @@ #ifndef __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_h #define __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_h -#if ITK_VERSION_MAJOR >= 4 - #include "itkImageToImageMetric.h" -#else - #include "itkOptImageToImageMetric.h" -#endif +#include "itkImageToImageMetric.h" #include "itkCovariantVector.h" #include "itkPoint.h" #include "itkIndex.h" diff --git a/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx b/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx index cf17340..dd7000d 100644 --- a/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx +++ b/registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx @@ -153,9 +153,6 @@ MeanSquaresImageToImageMetricFor3DBLUTFFD // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // MUST BE CALLED TO INITIATE PROCESSING this->GetValueMultiThreadedInitiate(); @@ -215,12 +212,8 @@ MeanSquaresImageToImageMetricFor3DBLUTFFD } // Jacobian should be evaluated at the unmapped (fixed image) point. -#if ITK_VERSION_MAJOR >= 4 TransformJacobianType jacobian; transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[fixedImageSample].point, jacobian); -#else - const TransformJacobianType & jacobian = transform ->GetJacobian( this->m_FixedImageSamples[fixedImageSample].point ); -#endif //double sum; unsigned int par, dim; for( par=0; parm_NumberOfParameters; par+=3) { @@ -257,9 +250,6 @@ MeanSquaresImageToImageMetricFor3DBLUTFFD // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); -#if ITK_VERSION_MAJOR < 4 - this->m_Parameters = parameters; -#endif // Reset the joint pdfs to zero memset( m_ThreaderMSE, diff --git a/segmentation/clitkAnatomicalFeatureDatabase.cxx b/segmentation/clitkAnatomicalFeatureDatabase.cxx index 0f5bfdd..13fb269 100644 --- a/segmentation/clitkAnatomicalFeatureDatabase.cxx +++ b/segmentation/clitkAnatomicalFeatureDatabase.cxx @@ -106,11 +106,7 @@ void clitk::AnatomicalFeatureDatabase::Load() //-------------------------------------------------------------------- void clitk::AnatomicalFeatureDatabase::SetPoint3D(std::string tag, PointType3D & p) { -#if ITK_VERSION_MAJOR > 3 std::ostringstream value; -#else - ::itk::OStringStream value; -#endif value << p[0] << " " << p[1] << " " << p[2]; m_MapOfTag[tag] = value.str(); } diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index b76001f..596604b 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -668,9 +668,7 @@ SearchForTracheaSeed2(int numberOfSlices) opening->Update(); typename SlicerFilterType::Pointer slicer = SlicerFilterType::New(); -#if ITK_VERSION_MAJOR >= 4 slicer->SetDirectionCollapseToIdentity(); -#endif slicer->SetInput(opening->GetOutput()); // label result @@ -752,11 +750,7 @@ SearchForTracheaSeed2(int numberOfSlices) for (unsigned int j = 0; j < nlables; j++) { shape = label_map->GetNthLabelObject(j); if (shape->Size() > 150 && shape->Size() <= max_size) { -#if ITK_VERSION_MAJOR < 4 - double e = 1 - 1/shape->GetBinaryElongation(); -#else double e = 1 - 1/shape->GetElongation(); -#endif //double area = 1 - r->Area() ; if (e < max_elongation) { nshapes++; diff --git a/segmentation/clitkFillMaskFilter.txx b/segmentation/clitkFillMaskFilter.txx index 7d2369c..cfe366a 100644 --- a/segmentation/clitkFillMaskFilter.txx +++ b/segmentation/clitkFillMaskFilter.txx @@ -141,9 +141,7 @@ GenerateData() start2D[m_Directions[i]]=sliceIndex; desiredRegion.SetIndex( start2D ); extractFilter->SetExtractionRegion( desiredRegion ); -#if ITK_VERSION_MAJOR == 4 extractFilter->SetDirectionCollapseToSubmatrix(); -#endif extractFilter->Update( ); typename ImageSliceType::Pointer slice= extractFilter->GetOutput(); diff --git a/segmentation/clitkMotionMaskGenericFilter.txx b/segmentation/clitkMotionMaskGenericFilter.txx index c6476c2..fd21b46 100644 --- a/segmentation/clitkMotionMaskGenericFilter.txx +++ b/segmentation/clitkMotionMaskGenericFilter.txx @@ -541,11 +541,8 @@ MotionMaskGenericFilter::InitializeEllips( typename itk::Vector= 4 typename LabelType::RegionType lung_bbox = label->GetBoundingBox(); -#else - typename LabelType::RegionType lung_bbox = label->GetRegion(); -#endif + axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2; axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2; axes[2] = 1.25*fabs(lung_centroid[2] - center[2]); diff --git a/segmentation/clitkRegionGrowingGenericFilter.txx b/segmentation/clitkRegionGrowingGenericFilter.txx index e8f3653..681b98e 100644 --- a/segmentation/clitkRegionGrowingGenericFilter.txx +++ b/segmentation/clitkRegionGrowingGenericFilter.txx @@ -110,24 +110,8 @@ void clitk::RegionGrowingGenericFilter::UpdateWithInputImageType() IteratorType it(ball.GetRadius(), input, input->GetLargestPossibleRegion()); -#if ITK_VERSION_MAJOR < 4 - typename BallType::ConstIterator nit; - unsigned idx = 0; - for (nit = ball.Begin(); nit != ball.End(); ++nit, ++idx) - { - if (*nit) - { - it.ActivateOffset(it.GetOffset(idx)); - } - else - { - it.DeactivateOffset(it.GetOffset(idx)); - } - } -#else it.CreateActiveListFromNeighborhood(ball); it.NeedToUseBoundaryConditionOff(); -#endif it.SetLocation(seeds[0]); for (typename IteratorType::ConstIterator i = it.Begin(); !i.IsAtEnd(); ++i) diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index fbab5f5..e0981f3 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -373,20 +373,15 @@ if(CLITK_BUILD_TOOLS) #========================================================= - if(ITK_VERSION_MAJOR VERSION_LESS 4) - message("clitkDice is not compatible with ITK<4. It will not be built.") - message("clitkDicomRTPlan2Gate is not compatible with GDCM<2 (ITK<4). It will not be built.") - else(ITK_VERSION_MAJOR VERSION_LESS 4) - WRAP_GGO(clitkDice_GGO_C clitkDice.ggo) - add_executable(clitkDice clitkDice.cxx ${clitkDice_GGO_C}) - target_link_libraries(clitkDice clitkSegmentationGgoLib clitkCommon ${ITK_LIBRARIES} ) - set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDice) - - WRAP_GGO(clitkDicomRTPlan2Gate_GGO_C clitkDicomRTPlan2Gate.ggo) - add_executable(clitkDicomRTPlan2Gate clitkDicomRTPlan2Gate.cxx clitkDicomRTPlan2Gate_ggo.c) - target_link_libraries(clitkDicomRTPlan2Gate clitkCommon) - set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDicomRTPlan2Gate) - endif(ITK_VERSION_MAJOR VERSION_LESS 4) + WRAP_GGO(clitkDice_GGO_C clitkDice.ggo) + add_executable(clitkDice clitkDice.cxx ${clitkDice_GGO_C}) + target_link_libraries(clitkDice clitkSegmentationGgoLib clitkCommon ${ITK_LIBRARIES} ) + set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDice) + + WRAP_GGO(clitkDicomRTPlan2Gate_GGO_C clitkDicomRTPlan2Gate.ggo) + add_executable(clitkDicomRTPlan2Gate clitkDicomRTPlan2Gate.cxx clitkDicomRTPlan2Gate_ggo.c) + target_link_libraries(clitkDicomRTPlan2Gate clitkCommon) + set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDicomRTPlan2Gate) #========================================================= diff --git a/tools/clitkImageToVectorImageGenericFilter.txx b/tools/clitkImageToVectorImageGenericFilter.txx index 8bbeab8..f159dc7 100644 --- a/tools/clitkImageToVectorImageGenericFilter.txx +++ b/tools/clitkImageToVectorImageGenericFilter.txx @@ -26,11 +26,7 @@ * @brief * ===================================================*/ -#if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2) -# include "itkCompose3DVectorImageFilter.h" -#else # include "itkComposeImageFilter.h" -#endif namespace clitk { @@ -82,11 +78,7 @@ namespace clitk typedef itk::Image, Dimension> OutputImageType; // Filter -#if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2) - typedef itk::Compose3DVectorImageFilter ComposeFilterType; -#else typedef itk::ComposeImageFilter ComposeFilterType; -#endif typename ComposeFilterType::Pointer composeFilter=ComposeFilterType::New(); // Read the inputs diff --git a/tools/clitkInvertVFGenericFilter.h b/tools/clitkInvertVFGenericFilter.h index 012524c..7765e56 100644 --- a/tools/clitkInvertVFGenericFilter.h +++ b/tools/clitkInvertVFGenericFilter.h @@ -36,11 +36,7 @@ //itk include #include "itkLightObject.h" -#if ITK_VERSION_MAJOR >= 4 - #include "itkInverseDisplacementFieldImageFilter.h" -#else - #include "itkInverseDeformationFieldImageFilter.h" -#endif +#include "itkInverseDisplacementFieldImageFilter.h" namespace clitk { diff --git a/tools/clitkInvertVFGenericFilter.txx b/tools/clitkInvertVFGenericFilter.txx index 812ed0c..a18a641 100644 --- a/tools/clitkInvertVFGenericFilter.txx +++ b/tools/clitkInvertVFGenericFilter.txx @@ -205,11 +205,7 @@ InvertVFGenericFilter::UpdateWithDimAndPixelType() case 2: { // Create the InverseDeformationFieldFilter -#if ITK_VERSION_MAJOR >= 4 typedef itk::InverseDisplacementFieldImageFilter FilterType; -#else - typedef itk::InverseDeformationFieldImageFilter FilterType; -#endif typename FilterType::Pointer filter =FilterType::New(); filter->SetInput(input); filter->SetOutputOrigin(input->GetOrigin()); diff --git a/tools/clitkJacobianImageGenericFilter.h b/tools/clitkJacobianImageGenericFilter.h index 73e9b64..1ddc301 100644 --- a/tools/clitkJacobianImageGenericFilter.h +++ b/tools/clitkJacobianImageGenericFilter.h @@ -35,11 +35,7 @@ //itk include #include "itkLightObject.h" -#if ITK_VERSION_MAJOR >= 4 - #include "itkInverseDisplacementFieldImageFilter.h" -#else - #include "itkInverseDeformationFieldImageFilter.h" -#endif +#include "itkInverseDisplacementFieldImageFilter.h" namespace clitk { diff --git a/tools/clitkSplitImageGenericFilter.cxx b/tools/clitkSplitImageGenericFilter.cxx index 72f712c..4eb00b9 100644 --- a/tools/clitkSplitImageGenericFilter.cxx +++ b/tools/clitkSplitImageGenericFilter.cxx @@ -105,9 +105,7 @@ void clitk::SplitImageGenericFilter::UpdateWithInputImageType() size[mSplitDimension]=0; typename ImageType::RegionType extracted_region; extracted_region.SetSize(size); -#if ITK_VERSION_MAJOR >= 4 filter->SetDirectionCollapseToIdentity(); -#endif filter->SetExtractionRegion(extracted_region); filter->Update(); diff --git a/tools/clitkVectorImageToImageFilter.h b/tools/clitkVectorImageToImageFilter.h index e2971f4..f2288f0 100644 --- a/tools/clitkVectorImageToImageFilter.h +++ b/tools/clitkVectorImageToImageFilter.h @@ -86,11 +86,7 @@ namespace clitk //---------------------------------------- // Update //---------------------------------------- -#if ITK_VERSION_MAJOR >= 4 void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId ); -#else - void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, int threadId ); -#endif //---------------------------------------- // Data members diff --git a/tools/clitkVectorImageToImageFilter.txx b/tools/clitkVectorImageToImageFilter.txx index 539e8dc..361f8dd 100644 --- a/tools/clitkVectorImageToImageFilter.txx +++ b/tools/clitkVectorImageToImageFilter.txx @@ -46,12 +46,7 @@ namespace clitk // Generate Data //------------------------------------------------------------------- template - void -#if ITK_VERSION_MAJOR >= 4 - VectorImageToImageFilter::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId) -#else - VectorImageToImageFilter::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, int threadId) -#endif + void VectorImageToImageFilter::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId) { // Iterators typename OutputImageType::Pointer output=this->GetOutput(); diff --git a/tools/clitkWarpImageGenericFilter.txx b/tools/clitkWarpImageGenericFilter.txx index cc3b17d..d1d1f71 100644 --- a/tools/clitkWarpImageGenericFilter.txx +++ b/tools/clitkWarpImageGenericFilter.txx @@ -205,11 +205,7 @@ WarpImageGenericFilter::UpdateWithDimAndPixelType() //Backward mapping typedef itk::WarpImageFilter BackwardWarpFilterType; typename BackwardWarpFilterType::Pointer backwardWarpFilter= BackwardWarpFilterType::New(); -#if ITK_VERSION_MAJOR >= 4 backwardWarpFilter->SetDisplacementField( deformationField ); -#else - backwardWarpFilter->SetDeformationField( deformationField ); -#endif backwardWarpFilter->SetEdgePaddingValue( static_cast(m_ArgsInfo.pad_arg) ); backwardWarpFilter->SetOutputSpacing( deformationField->GetSpacing() ); backwardWarpFilter->SetOutputOrigin( input->GetOrigin() ); diff --git a/vv/vvImageWarp.cxx b/vv/vvImageWarp.cxx index ff65803..2890ec8 100644 --- a/vv/vvImageWarp.cxx +++ b/vv/vvImageWarp.cxx @@ -101,11 +101,7 @@ void vvImageWarp::Update_WithDimAndPixelType() jacobian_filter->SetUseImageSpacingOn(); vf_connector->SetInput(mVF->GetVTKImages()[num]); warp_filter->SetInput(input[num]); -#if ITK_VERSION_MAJOR >= 4 warp_filter->SetDisplacementField(vf_connector->GetOutput()); -#else - warp_filter->SetDeformationField(vf_connector->GetOutput()); -#endif jacobian_filter->SetInput(vf_connector->GetOutput()); warp_filter->SetOutputSpacing(input[num]->GetSpacing()); warp_filter->SetOutputOrigin(input[num]->GetOrigin()); diff --git a/vv/vvMidPosition.cxx b/vv/vvMidPosition.cxx index b4582e2..2290d87 100644 --- a/vv/vvMidPosition.cxx +++ b/vv/vvMidPosition.cxx @@ -108,11 +108,7 @@ vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,in typename FilterType::Pointer warp_filter = FilterType::New(); warp_filter->SetInput(input); -#if ITK_VERSION_MAJOR >= 4 warp_filter->SetDisplacementField(resampler->GetOutput()); -#else - warp_filter->SetDeformationField(resampler->GetOutput()); -#endif warp_filter->SetOutputSpacing(input->GetSpacing()); warp_filter->SetOutputOrigin(input->GetOrigin()); warp_filter->SetOutputSize(input->GetLargestPossibleRegion().GetSize());