#=========================================================
### 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()
+
#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<itk::GDCMImageIOFactory *>(*it))
itk::PNGImageIOFactory::UnRegisterFactory(*it);
break;
}
-#endif
#if CLITK_PRIVATE_FEATURES
clitk::UsfImageIOFactory::RegisterOneFactory();
clitk::USVoxImageIOFactory::RegisterOneFactory();
#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();
rtk::ImagXImageIOFactory::RegisterOneFactory();
rtk::XRadImageIOFactory::RegisterOneFactory();
clitk::EsrfHstImageIOFactory::RegisterOneFactory();
-#if ITK_VERSION_MAJOR >= 4
itk::GDCMImageIOFactory::RegisterOneFactory();
itk::PNGImageIOFactory::RegisterOneFactory();
-#endif
} ////
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();
filter->SetExtractionRegion(extractedRegion);
filter->SetInput(reader->GetOutput());
filter->ReleaseDataFlagOn();
-#if ITK_VERSION_MAJOR == 4
filter->SetDirectionCollapseToSubmatrix();
-#endif
try {
mImage->AddItkImage<SlicedImageType>(filter->GetOutput());
mImage->ComputeScalarRangeBase<InputPixelType, VImageDimension-1>(filter->GetOutput());
* This filter is implemented using the propagation algorithm
*/
-#if ITK_VERSION_MAJOR == 4
template <class TInputImage, class TOutputImage, class TtNorm=Functor::Minimum<
typename TOutputImage::PixelType,
typename TOutputImage::PixelType,
typename TOutputImage::PixelType> >
-#else
- template <class TInputImage, class TOutputImage, class TtNorm=Function::Minimum<
- typename TOutputImage::PixelType,
- typename TOutputImage::PixelType,
- typename TOutputImage::PixelType> >
-#endif
class ITK_EXPORT RelativePositionPropImageFilter :
public ImageToImageFilter< TInputImage, TOutputImage >
{
#include <itkBinaryErodeImageFilter.h>
#include <itkBinaryBallStructuringElement.h>
#include <itkAddImageFilter.h>
-#if ITK_VERSION_MAJOR >= 4
- #include <itkDivideImageFilter.h>
-#else
- #include <itkDivideByConstantImageFilter.h>
-#endif
+#include <itkDivideImageFilter.h>
// itk [Bloch et al]
#include "RelativePositionPropImageFilter.h"
// Divide by the number of relpos
if (GetNumberOfAngles() != 1) {
-#if ITK_VERSION_MAJOR >= 4
typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
typename DivideFilter::Pointer divideFilter = DivideFilter::New();
divideFilter->SetConstant2(GetNumberOfAngles());
-#else
- typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
- typename DivideFilter::Pointer divideFilter = DivideFilter::New();
- divideFilter->SetConstant(GetNumberOfAngles());
-#endif
divideFilter->SetInput(m_FuzzyMap);
divideFilter->Update();
m_FuzzyMap = divideFilter->GetOutput();
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
//-----------------------------------------------------------------------
template <class InputImageType, class OutputImageType>
void
- BackProjectImageFilter<InputImageType, OutputImageType>
-#if ITK_VERSION_MAJOR >= 4
- ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId )
-#else
- ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, int threadId )
-#endif
+ BackProjectImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId )
{
//Projection pointer
InputImageConstPointer inputPtr=this->GetInput();
//========================================================================================
//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;
//=========================================================================================================================
//update the output for the outputRegionForThread
-#if ITK_VERSION_MAJOR >= 4
template<class InputImageType, class OutputImageType>
void ComposeVFFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
-#else
- template<class InputImageType, class OutputImageType>
- void ComposeVFFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
-#endif
{
//Get pointer to the output
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;
*/
template <class TInputImage, class TOutputImage>
void
-ExtractImageFilter<TInputImage,TOutputImage>
-#if ITK_VERSION_MAJOR >= 4
-::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
-#else
-::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId)
-#endif
+ExtractImageFilter<TInputImage,TOutputImage>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
itkDebugMacro(<<"Actually executing");
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
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());
}
//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;
//=========================================================================================================================
//update the output for the outputRegionForThread
template<class InputImageType, class OutputImageType>
-#if ITK_VERSION_MAJOR >= 4
void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
-#else
-void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
-#endif
{
// std::cout << "HelperClass1::ThreadedGenerateData - IN " << threadId << std::endl;
//Get pointer to the input
//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;
//=========================================================================================================================
//update the output for the outputRegionForThread
-#if ITK_VERSION_MAJOR >= 4
template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
-#else
-template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
-#endif
{
// std::cout << "HelperClass2::ThreadedGenerateData - IN " << threadId << std::endl;
#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
typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
typedef itk::BinaryBallStructuringElement<InternalPixelType,InputImageDimension > KernelType;
typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType , KernelType> ConditionalBinaryDilateImageFilterType;
-#if ITK_VERSION_MAJOR >= 4
typedef itk::Testing::ComparisonImageFilter<InternalImageType, InternalImageType> DifferenceImageFilterType;
-#else
- typedef itk::DifferenceImageFilter<InternalImageType, InternalImageType> DifferenceImageFilterType;
-#endif
typedef itk::CastImageFilter<InternalImageType, OutputImageType> OutputCastImageFilterType;
typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType> SetBackgroundImageFilterType;
*
* \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
template <class TInputImage1, class TInputImage2, class TOutputImage, class TFunction >
void
FlexibleBinaryFunctorImageFilter<TInputImage1, TInputImage2, TOutputImage, TFunction>
-::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;
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;
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)..."<<std::endl;
-#endif
//============================================================================
// Initialize using image moments.
if (m_Verbose) std::cout << "Starting the registration now..." << std::endl;
try {
-#if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2)
- registration->StartRegistration();
-#else
registration->Update();
-#endif
} catch ( itk::ExceptionObject & err ) {
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
// Crop the fixedImage to the bounding box to facilitate multi-resolution
typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
-#if ITK_VERSION_MAJOR == 4
extractImageFilter->SetDirectionCollapseToSubmatrix();
-#endif
extractImageFilter->SetInput(fixedImage);
extractImageFilter->SetExtractionRegion(transformRegion);
extractImageFilter->Update();
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)..."<<std::endl;
-#endif
-
//=======================================================
// Optimizer
try
{
-#if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2)
- registration->StartRegistration();
-#else
registration->Update();
-#endif
}
catch( itk::ExceptionObject & err )
{
# else
typedef itk::TransformToDisplacementFieldFilter<DisplacementFieldType, double> ConvertorType;
# endif
-#else
- typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
#endif
typename ConvertorType::Pointer filter= ConvertorType::New();
filter->SetNumberOfThreads(1);
/** 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;
}
/** 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;
// VD Add MultipleBSplineDeformableTransform as friend to facilitate wrapping
friend class MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>;
-#if ITK_VERSION_MAJOR >= 4
mutable JacobianType m_SharedDataBSplineJacobian;
-#endif
-
}; //class BSplineDeformableTransform
// Constructor with default arguments
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
- ::BSplineDeformableTransform():Superclass(0)
-#else
- ::BSplineDeformableTransform():Superclass(OutputDimension,0)
-#endif
+ BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::BSplineDeformableTransform():Superclass(0)
{
unsigned int i;
// Get the number of parameters
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
-#else
- unsigned int
-#endif
BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::GetNumberOfParameters(void) const
{
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<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
-#else
- JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
-#endif
+
memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
m_LastJacobianIndex = m_ValidRegion.GetIndex();
m_NeedResetJacobian = false;
// 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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
void
BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
jacobian = m_SharedDataBSplineJacobian;
}
-#else
- template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- const
- typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::JacobianType &
- BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::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<JacobianValueType>::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
# else
# include "itkTransformToDisplacementFieldFilter.h"
# endif
-#else
-# include "itkTransformToDeformationFieldSource.h"
#endif
namespace clitk
# else
typedef itk::TransformToDisplacementFieldFilter<OutputImageType, double> ConvertorType;
# endif
-#else
- typedef itk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
#endif
itkNewMacro(Self);
typedef clitk::VectorImageToImageFilter<BLUTCoefficientImageType, typename ITKTransformType::ImageType> 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();
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.
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;
}
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++;
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++;
return OutputCovariantVectorType();
}
-#if ITK_VERSION_MAJOR >= 4
virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const
{
itkExceptionMacro( << "DeformationFieldTransform doesn't declare ComputeJacobianWithRespectToParameters" );
{
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();
// Constructor
template<class TScalarType, unsigned int InputDimension, unsigned int OutputDimension, unsigned int SpaceDimension>
- DeformationFieldTransform<TScalarType, InputDimension, OutputDimension, SpaceDimension>
-#if ITK_VERSION_MAJOR >= 4
- ::DeformationFieldTransform():Superclass(1)
-#else
- ::DeformationFieldTransform():Superclass(OutputDimension,1)
-#endif
+ DeformationFieldTransform<TScalarType, InputDimension, OutputDimension, SpaceDimension>::DeformationFieldTransform():Superclass(1)
+
{
m_DeformationField=NULL;
m_Interpolator=DefaultInterpolatorType::New();
//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<FixedImageType, MovingImageType, DisplacementFieldType> MultiResolutionRegistrationType;
typedef CommandResolutionLevelUpdate<MultiResolutionRegistrationType> LevelObserver;
//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 );
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() );
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
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;
}
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) {
if (m_Verbose) std::cout<<"number of mask pixels "<<totalNumberOfMaskPixels<<std::endl;
}
-
-#else
- if (m_Verbose) std::cout<<"Not setting the fixed image intensity threshold or the fraction of samples to use (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
-
-
-#endif
//============================================================================
//return the pointer
return m_Metric;
bool m_OptimizerInitialized;
InternalOptimizerType * m_VnlOptimizer;
-#if ITK_VERSION_MAJOR > 3
mutable std::ostringstream m_StopConditionDescription;
-#else
- mutable itk::OStringStream m_StopConditionDescription;
-#endif
BoundValueType m_LowerBound;
BoundValueType m_UpperBound;
BoundSelectionType m_BoundSelection;
# else
# include "itkTransformToDisplacementFieldFilter.h"
# endif
-#else
-# include "itkTransformToDeformationFieldSource.h"
#endif
#include "itkAffineTransform.h"
# else
typedef itk::TransformToDisplacementFieldFilter<OutputImageType, double> ConvertorType;
# endif
-#else
- typedef itk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
#endif
typename ConvertorType::Pointer filter= ConvertorType::New();
if( tempField.IsNull() )
{
-#if ITK_VERSION_MAJOR >= 4
m_RegistrationFilter->SetInitialDisplacementField( NULL );
-#else
- m_RegistrationFilter->SetInitialDeformationField( NULL );
-#endif
}
else
{
tempField = m_FieldExpander->GetOutput();
tempField->DisconnectPipeline();
-#if ITK_VERSION_MAJOR >= 4
m_RegistrationFilter->SetInitialDisplacementField( tempField );
-#else
- m_RegistrationFilter->SetInitialDeformationField( tempField );
-#endif
-
}
// setup registration filter and pyramids
/** 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;
}
/** 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;
std::vector<ParametersType> m_parameters;
mutable std::vector<CoefficientImagePointer> m_CoefficientImages;
mutable int m_LastJacobian;
-#if ITK_VERSION_MAJOR >= 4
mutable JacobianType m_SharedDataBSplineJacobian;
-#endif
void InitJacobian();
// FIXME it seems not used
{
// Constructor with default arguments
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
- ::MultipleBSplineDeformableTransform() : Superclass(0)
-#else
- ::MultipleBSplineDeformableTransform() : Superclass(OutputDimension, 0)
-#endif
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::MultipleBSplineDeformableTransform() : Superclass(0)
{
m_nLabels = 1;
m_labels = 0;
#undef LOOP_ON_LABELS
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
-#else
- inline unsigned int
-#endif
MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::GetNumberOfParameters(void) const
{
return m_trans[lidx]->DeformablyTransformPoint(inputPoint);
}
-#if ITK_VERSION_MAJOR >= 4
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
inline void
MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
jacobian = this->m_SharedDataBSplineJacobian;
}
-#else
- template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::JacobianType &
- MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
inline void
MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>::InitJacobian()
{
unsigned numberOfParameters = this->GetNumberOfParameters();
-#if ITK_VERSION_MAJOR >= 4
this->m_SharedDataBSplineJacobian.set_size(OutputDimension, numberOfParameters);
JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
-#else
- this->m_Jacobian.set_size(OutputDimension, numberOfParameters);
- JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
-#endif
memset(jacobianDataPointer, 0, numberOfParameters * sizeof (JacobianPixelType));
unsigned tot = 0;
// 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<TFixedImage, TMovingImage > Superclass;
-
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> 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
// 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 <class TFixedImage, class TMovingImage>
-NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
-::NormalizedCorrelationImageToImageMetric()
-{
- m_SubtractMean = false;
-}
-
-/*
- * Get the match Measure
- */
-template <class TFixedImage, class TMovingImage>
-typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
-NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
-::GetValue( const TransformParametersType & parameters ) const
-{
-
- FixedImageConstPointer fixedImage = this->m_FixedImage;
-
- if( !fixedImage ) {
- itkExceptionMacro( << "Fixed image has not been assigned" );
- }
-
- typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
-
- FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
-
- typename FixedImageType::IndexType index;
-
- MeasureType measure;
-
- 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<TFixedImage,TMovingImage>
-::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<FixedImageType> 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<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeF = DerivativeType( ParametersDimension );
- derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeM = DerivativeType( ParametersDimension );
- derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::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<CoordRepType,MovingImageType::ImageDimension>
- 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<ParametersDimension; par++) {
- RealType sumF = itk::NumericTraits< RealType >::Zero;
- RealType sumM = itk::NumericTraits< RealType >::Zero;
- for(unsigned int dim=0; dim<dimension; dim++) {
- const RealType differential = jacobian( dim, par ) * gradient[dim];
- sumF += fixedValue * differential;
- sumM += movingValue * differential;
- if ( this->m_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<ParametersDimension; i++) {
- derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
- }
- } else {
- for(unsigned int i=0; i<ParametersDimension; i++) {
- derivative[i] = itk::NumericTraits< MeasureType >::Zero;
- }
- }
-
-}
-
-
-/*
- * Get both the match Measure and theDerivative Measure
- */
-template <class TFixedImage, class TMovingImage>
-void
-NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
-::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<FixedImageType> 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<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeF = DerivativeType( ParametersDimension );
- derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeM = DerivativeType( ParametersDimension );
- derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeM1 = DerivativeType( ParametersDimension );
- derivativeM1.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::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<CoordRepType,MovingImageType::ImageDimension>
- 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<ParametersDimension; par++) {
- RealType sumF = itk::NumericTraits< RealType >::Zero;
- RealType sumM = itk::NumericTraits< RealType >::Zero;
- for(unsigned int dim=0; dim<dimension; dim++) {
- const RealType differential = jacobian( dim, par ) * gradient[dim];
- sumF += fixedValue * differential;
- sumM += movingValue * differential;
- if ( this->m_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<ParametersDimension; i++) {
- derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
- }
- value = sfm / denom;
- } else {
- for(unsigned int i=0; i<ParametersDimension; i++) {
- derivative[i] = itk::NumericTraits< MeasureType >::Zero;
- }
- value = itk::NumericTraits< MeasureType >::Zero;
- }
-
-
-
-}
-
-template < class TFixedImage, class TMovingImage>
-void
-NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
-::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
// 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<TFixedImage, TMovingImage > Superclass;
-
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> 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
// 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 <class TFixedImage, class TMovingImage>
-NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD()
-{
- m_SubtractMean = false;
-}
-
-/*
- * Get the match Measure
- */
-template <class TFixedImage, class TMovingImage>
-typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>::MeasureType
-NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::GetValue( const TransformParametersType & parameters ) const
-{
-
- FixedImageConstPointer fixedImage = this->m_FixedImage;
-
- if( !fixedImage ) {
- itkExceptionMacro( << "Fixed image has not been assigned" );
- }
-
- typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
-
- FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
-
- typename FixedImageType::IndexType index;
-
- MeasureType measure;
-
- 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<TFixedImage,TMovingImage>
-::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<FixedImageType> 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<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeF = DerivativeType( ParametersDimension );
- derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeM = DerivativeType( ParametersDimension );
- derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::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<CoordRepType,MovingImageType::ImageDimension>
- 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<ParametersDimension; par++) {
- RealType sumF = itk::NumericTraits< RealType >::Zero;
- RealType sumM = itk::NumericTraits< RealType >::Zero;
- for(unsigned int dim=0; dim<dimension; dim++) {
- const RealType differential = jacobian( dim, par ) * gradient[dim];
- sumF += fixedValue * differential;
- sumM += movingValue * differential;
- if ( this->m_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<ParametersDimension; i++) {
- derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
- }
- } else {
- for(unsigned int i=0; i<ParametersDimension; i++) {
- derivative[i] = itk::NumericTraits< MeasureType >::Zero;
- }
- }
-
-}
-
-
-/*
- * Get both the match Measure and theDerivative Measure
- */
-template <class TFixedImage, class TMovingImage>
-void
-NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::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<FixedImageType> 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<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeF = DerivativeType( ParametersDimension );
- derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeM = DerivativeType( ParametersDimension );
- derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
-
- DerivativeType derivativeM1 = DerivativeType( ParametersDimension );
- derivativeM1.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::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<CoordRepType,MovingImageType::ImageDimension>
- 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<ParametersDimension; par++) {
- RealType sumF = itk::NumericTraits< RealType >::Zero;
- RealType sumM = itk::NumericTraits< RealType >::Zero;
- for(unsigned int dim=0; dim<dimension; dim++) {
- const RealType differential = jacobian( dim, par ) * gradient[dim];
- sumF += fixedValue * differential;
- sumM += movingValue * differential;
- if ( this->m_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<ParametersDimension; i++) {
- derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
- }
- value = sfm / denom;
- } else {
- for(unsigned int i=0; i<ParametersDimension; i++) {
- derivative[i] = itk::NumericTraits< MeasureType >::Zero;
- }
- value = itk::NumericTraits< MeasureType >::Zero;
- }
-
-
-
-}
-
-template < class TFixedImage, class TMovingImage>
-void
-NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::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
#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"
// 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();
// 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();
}
// 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; par<this->m_NumberOfParameters; par++) {
RealType sumF = itk::NumericTraits< RealType >::Zero;
// 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);
#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"
// 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();
// 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();
}
// 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; par<this->m_NumberOfParameters; par++)
// {
// 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);
return OutputCovariantVectorType();
}
-#if ITK_VERSION_MAJOR >= 4
virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const
{
itkExceptionMacro( << "PointListTransform doesn't declare ComputeJacobianWithRespectToParameters" );
{
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();
// Constructor
template<class TScalarType, unsigned int NDimensions, unsigned int NOutputDimensions>
- PointListTransform<TScalarType, NDimensions, NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
- ::PointListTransform():Superclass(1)
-#else
- ::PointListTransform():Superclass(NOutputDimensions,1)
-#endif
+ PointListTransform<TScalarType, NDimensions, NOutputDimensions>::PointListTransform():Superclass(1)
{
m_PointLists.resize(0);
m_PointList.resize(1);
/** 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;
}
/** 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;
// JV Shape
unsigned int m_TransformShape;
-#if ITK_VERSION_MAJOR >= 4
mutable JacobianType m_SharedDataBSplineJacobian;
-#endif
}; //class ShapedBLUTSpatioTemporalDeformableTransform
// Constructor with default arguments
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
- ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0)
-#else
- ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(OutputDimension,0)
-#endif
+ ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0)
{
unsigned int i;
// Get the number of parameters
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
-#else
- unsigned int
-#endif
ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::GetNumberOfParameters(void) const
{
//=====================================
//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<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
-#else
- JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
-#endif
memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
for (unsigned int j=0; j<OutputDimension; j++)
// JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
// Compute the Jacobian in one position
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
-#if ITK_VERSION_MAJOR >= 4
void
ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::ComputeJacobianWithRespectToParameters( const InputPointType & point, JacobianType & jacobian) const
-#else
- const
- typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::JacobianType &
- ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::GetJacobian( const InputPointType & point ) const
-#endif
{
//========================================================
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
// 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
}
// Return the result
-#if ITK_VERSION_MAJOR >= 4
jacobian = m_SharedDataBSplineJacobian;
-#else
- return this->m_Jacobian;
-#endif
}
// 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 TFixedImage,class TMovingImage >
-class ITK_EXPORT MattesMutualInformationImageToImageMetricFor3DBLUTFFD :
- public ImageToImageMetric< TFixedImage, TMovingImage >
-{
-public:
-
- /** Standard class typedefs. */
- typedef MattesMutualInformationImageToImageMetricFor3DBLUTFFD Self;
- typedef ImageToImageMetric< TFixedImage, TMovingImage > Superclass;
- typedef SmartPointer<Self> Pointer;
- typedef SmartPointer<const Self> 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<unsigned long>::max() );
- itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long);
-
- /** Number of bins to used in the histogram. Typical value is 50. */
- itkSetClampMacro( NumberOfHistogramBins, unsigned long,
- 1, NumericTraits<unsigned long>::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<FixedImageSpatialSample>
- 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<PDFValueType> 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<PDFValueType,2> JointPDFType;
- typedef JointPDFType::IndexType JointPDFIndexType;
- typedef JointPDFType::PixelType JointPDFValueType;
- typedef JointPDFType::RegionType JointPDFRegionType;
- typedef JointPDFType::SizeType JointPDFSizeType;
- typedef Image<PDFValueType,3> 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<MovingImageType,
- CoordinateRepresentationType>
- BSplineInterpolatorType;
-
- /** Pointer to BSplineInterpolator. */
- typename BSplineInterpolatorType::Pointer m_BSplineInterpolator;
-
- /** Typedefs for using central difference calculator. */
- typedef CentralDifferenceImageFunction<MovingImageType,
- CoordinateRepresentationType>
- 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<FixedImageType>::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<WeightsValueType> BSplineTransformWeightsArrayType;
- typedef typename BSplineTransformIndexArrayType::ValueType IndexValueType;
- typedef Array2D<IndexValueType> BSplineTransformIndicesArrayType;
- typedef std::vector<MovingImagePointType> MovingImagePointArrayType;
- typedef std::vector<bool> BooleanArrayType;
-
- BSplineTransformWeightsArrayType m_BSplineTransformWeightsArray;
- BSplineTransformIndicesArrayType m_BSplineTransformIndicesArray;
- MovingImagePointArrayType m_PreTransformPointsArray;
- BooleanArrayType m_WithinSupportRegionArray;
-
- typedef FixedArray<unsigned long,
- ::itk::GetImageDimension<FixedImageType>::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
// 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<TFixedImage,TMovingImage>
-::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<TFixedImage,TMovingImage>
-::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 <class TFixedImage, class TMovingImage>
-void
-MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::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<double>::max();
- double fixedImageMax = NumericTraits<double>::NonpositiveMin();
-
- typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
- FixedIteratorType fixedImageIterator(
- this->m_FixedImage, this->GetFixedImageRegion() );
-
- for ( fixedImageIterator.GoToBegin();
- !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
-
- double sample = static_cast<double>( fixedImageIterator.Get() );
-
- if ( sample < fixedImageMin ) {
- fixedImageMin = sample;
- }
-
- if ( sample > fixedImageMax ) {
- fixedImageMax = sample;
- }
- }
-
-
- /**
- * Compute the minimum and maximum for the entire moving image
- * in the buffer.
- */
- double movingImageMin = NumericTraits<double>::max();
- double movingImageMax = NumericTraits<double>::NonpositiveMin();
-
- typedef ImageRegionConstIterator<MovingImageType> MovingIteratorType;
- MovingIteratorType movingImageIterator(
- this->m_MovingImage, this->m_MovingImage->GetBufferedRegion() );
-
- for ( movingImageIterator.GoToBegin();
- !movingImageIterator.IsAtEnd(); ++movingImageIterator) {
- double sample = static_cast<double>( 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<double>( m_NumberOfHistogramBins - 2 * padding );
- m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize -
- static_cast<double>( padding );
-
- m_MovingImageBinSize = ( movingImageMax - movingImageMin ) /
- static_cast<double>( m_NumberOfHistogramBins - 2 * padding );
- m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize -
- static_cast<double>( 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<BSplineInterpolatorType *>(
- 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<BSplineTransformType *>(
- 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<TFixedImage,TMovingImage>
-::SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples )
-{
-
- // Set up a random interator within the user specified fixed image region.
- typedef ImageRandomConstIteratorWithIndex<FixedImageType> 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<TFixedImage,TMovingImage>
-::SampleFullFixedImageDomain( FixedImageSpatialSampleContainer& samples )
-{
-
- // Set up a region interator within the user specified fixed image region.
- typedef ImageRegionConstIteratorWithIndex<FixedImageType> 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<TFixedImage,TMovingImage>
-::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<double>( (*iter).FixedImageValue ) / m_FixedImageBinSize -
- m_FixedImageNormalizedMin;
- unsigned int pindex = static_cast<unsigned int>( 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<TFixedImage,TMovingImage>
-::MeasureType
-MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::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<unsigned int>( 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<PDFValueType>( 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<int>( movingImageParzenWindowIndex ) - 1;
- pdfPtr += pdfMovingIndex;
-
- for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
- + 2;
- pdfMovingIndex++, pdfPtr++ ) {
-
- // Update PDF for the current intensity pair
- double movingImageParzenWindowArg =
- static_cast<double>( pdfMovingIndex ) -
- static_cast<double>( movingImageParzenWindowTerm );
-
- *(pdfPtr) += static_cast<PDFValueType>(
- 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<JointPDFType> 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<PDFValueType>( 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<PDFValueType>( fixedPDFSum );
- }
-
-
- // Compute moving image marginal PDF by summing over fixed image bins.
- typedef ImageLinearIteratorWithIndex<JointPDFType> 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<PDFValueType>(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<MeasureType>( -1.0 * sum );
-
-}
-
-
-/**
- * Get the both Value and Derivative Measure
- */
-template < class TFixedImage, class TMovingImage >
-void
-MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::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<unsigned int>( 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<PDFValueType>( 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<int>( movingImageParzenWindowIndex ) - 1;
- pdfPtr += pdfMovingIndex;
-
- for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
- + 2;
- pdfMovingIndex++, pdfPtr++ ) {
- // Update PDF for the current intensity pair
- double movingImageParzenWindowArg =
- static_cast<double>( pdfMovingIndex ) -
- static_cast<double>(movingImageParzenWindowTerm);
-
- *(pdfPtr) += static_cast<PDFValueType>(
- 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<JointPDFType> 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<PDFValueType>( 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<PDFValueType>( fixedPDFSum );
- }
-
-
- // Compute moving image marginal PDF by summing over fixed image bins.
- typedef ImageLinearIteratorWithIndex<JointPDFType> 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<PDFValueType>(sum);
-
- linearIter.NextLine();
- ++movingIndex1;
-
- }
-
- double nFactor = 1.0 / ( m_MovingImageBinSize
- * static_cast<double>( nSamples ) );
-
- if( this->m_UseExplicitPDFDerivatives ) {
- // Normalize the joint PDF derivatives by the test image binsize and nSamples
- typedef ImageRegionIterator<JointPDFDerivativesType>
- 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<unsigned int>( 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<int>( movingImageParzenWindowIndex ) - 1;
-
- for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
- + 2;
- pdfMovingIndex++ ) {
-
- // Update PDF for the current intensity pair
- double movingImageParzenWindowArg =
- static_cast<double>( pdfMovingIndex ) -
- static_cast<double>(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<MeasureType>( -1.0 * sum );
-}
-
-
-/**
- * Get the match measure derivative
- */
-template < class TFixedImage, class TMovingImage >
-void
-MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::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<TFixedImage,TMovingImage>
-::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<TFixedImage,TMovingImage>
-::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<TFixedImage,TMovingImage>
-::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<TFixedImage,TMovingImage>
-::ReinitializeSeed()
-{
- Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed();
-}
-
-// Method to reinitialize the seed of the random number generator
-template < class TFixedImage, class TMovingImage > void
-MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
-::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<TFixedImage,TMovingImage>
-::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
// 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<TFixedImage, TMovingImage > Superclass;
- typedef SmartPointer<Self> Pointer;
- typedef SmartPointer<const Self> 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<typename TMovingImage::PixelType>));
- itkConceptMacro(MovingRealTypeAdditivieOperatorsCheck,
- (Concept::AdditiveOperators<
- typename NumericTraits<typename TMovingImage::PixelType>::RealType>));
- itkConceptMacro(MovingRealTypeMultiplyOperatorCheck,
- (Concept::MultiplyOperator<
- typename NumericTraits<typename TMovingImage::PixelType>::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
// 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"
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++;
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++;
#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"
{
// 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();
// 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();
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++ )
// {
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++ ) {
#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"
// 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();
}
// 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; par<this->m_NumberOfParameters; par+=3) {
// 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,
//--------------------------------------------------------------------
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();
}
opening->Update();
typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
-#if ITK_VERSION_MAJOR >= 4
slicer->SetDirectionCollapseToIdentity();
-#endif
slicer->SetInput(opening->GetOutput());
// label result
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++;
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();
// try to guess ideal ellipse axes from the lungs' bounding box and centroid
// with some hard-coded "magic" constants...
-#if ITK_VERSION_MAJOR >= 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]);
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)
#=========================================================
- 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)
#=========================================================
* @brief
*
===================================================*/
-#if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2)
-# include "itkCompose3DVectorImageFilter.h"
-#else
# include "itkComposeImageFilter.h"
-#endif
namespace clitk
{
typedef itk::Image<itk::Vector<PixelType,3>, Dimension> OutputImageType;
// Filter
-#if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2)
- typedef itk::Compose3DVectorImageFilter<InputImageType,OutputImageType> ComposeFilterType;
-#else
typedef itk::ComposeImageFilter<InputImageType,OutputImageType> ComposeFilterType;
-#endif
typename ComposeFilterType::Pointer composeFilter=ComposeFilterType::New();
// Read the inputs
//itk include
#include "itkLightObject.h"
-#if ITK_VERSION_MAJOR >= 4
- #include "itkInverseDisplacementFieldImageFilter.h"
-#else
- #include "itkInverseDeformationFieldImageFilter.h"
-#endif
+#include "itkInverseDisplacementFieldImageFilter.h"
namespace clitk
{
case 2: {
// Create the InverseDeformationFieldFilter
-#if ITK_VERSION_MAJOR >= 4
typedef itk::InverseDisplacementFieldImageFilter<InputImageType,OutputImageType> FilterType;
-#else
- typedef itk::InverseDeformationFieldImageFilter<InputImageType,OutputImageType> FilterType;
-#endif
typename FilterType::Pointer filter =FilterType::New();
filter->SetInput(input);
filter->SetOutputOrigin(input->GetOrigin());
//itk include
#include "itkLightObject.h"
-#if ITK_VERSION_MAJOR >= 4
- #include "itkInverseDisplacementFieldImageFilter.h"
-#else
- #include "itkInverseDeformationFieldImageFilter.h"
-#endif
+#include "itkInverseDisplacementFieldImageFilter.h"
namespace clitk
{
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();
//----------------------------------------
// 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
// Generate Data
//-------------------------------------------------------------------
template<class InputImageType, class OutputImageType>
- void
-#if ITK_VERSION_MAJOR >= 4
- VectorImageToImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId)
-#else
- VectorImageToImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, int threadId)
-#endif
+ void VectorImageToImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId)
{
// Iterators
typename OutputImageType::Pointer output=this->GetOutput();
//Backward mapping
typedef itk::WarpImageFilter<InputImageType, InputImageType, DeformationFieldType> BackwardWarpFilterType;
typename BackwardWarpFilterType::Pointer backwardWarpFilter= BackwardWarpFilterType::New();
-#if ITK_VERSION_MAJOR >= 4
backwardWarpFilter->SetDisplacementField( deformationField );
-#else
- backwardWarpFilter->SetDeformationField( deformationField );
-#endif
backwardWarpFilter->SetEdgePaddingValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
backwardWarpFilter->SetOutputSpacing( deformationField->GetSpacing() );
backwardWarpFilter->SetOutputOrigin( input->GetOrigin() );
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());
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());
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff