Hi, I am using the following code for registering two sets of points. my input fixed and moving images are in .mat file but code is giving me the error. does anyone know the solution ?
// Last used for diaphragm registration for testing the data with artifact (004...)
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImageRegistrationMethod.h"
#include "itkMattesMutualInformationImageToImageMetric.h" //Mattes Mutual information as the metric
#include "itkLinearInterpolateImageFunction.h" // Linear interpolation
#include "itkImage.h"
#include "itkTimeProbesCollectorBase.h"
#ifdef ITK_USE_REVIEW
#include "itkMemoryProbesCollectorBase.h"
#define itkProbesCreate() \
itk::TimeProbesCollectorBase chronometer; \
itk::MemoryProbesCollectorBase memorymeter
#define itkProbesStart( text ) memorymeter.Start( text ); chronometer.Start( text )
#define itkProbesStop( text ) chronometer.Stop( text ); memorymeter.Stop( text )
#define itkProbesReport( stream ) chronometer.Report( stream ); memorymeter.Report( stream )
#else
#define itkProbesCreate() \
itk::TimeProbesCollectorBase chronometer
#define itkProbesStart( text ) chronometer.Start( text )
#define itkProbesStop( text ) chronometer.Stop( text )
#define itkProbesReport( stream ) chronometer.Report( stream )
#endif
// initializer for the affine/rigid registration
#include "itkCenteredTransformInitializer.h"
// the rigid transform
#include "itkVersorRigid3DTransform.h"
// the affine transform
#include "itkAffineTransform.h"
// the BSpline transform
#include "itkBSplineDeformableTransform.h"
// Gradient Descent
#include "itkRegularStepGradientDescentOptimizer.h"
// Resample
#include "itkBSplineResampleImageFunction.h"
#include "itkBSplineDecompositionImageFilter.h"
// read write filters
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
////////////////////////////////////////////////////////////////////////
/////////////////////////// //////////////////////////////////
//#include <itkNearestNeighborInterpolateImageFunction.h>
//************************************************
//#include "itkSquaredDifferenceImageFilter.h"
#include "itkAbsoluteValueDifferenceImageFilter.h"
//************************************************
#include "itkTransformFileReader.h"
char * metricFilename;
// The following section of code implements a Command observer
// used to monitor the evolution of the registration process.
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef const OptimizerType * OptimizerPointer;
void Execute(itk::Object *caller, const itk::EventObject & event)
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object * object, const itk::EventObject & event)
{
OptimizerPointer optimizer =
dynamic_cast< OptimizerPointer >( object );
if( !(itk::IterationEvent().CheckEvent( &event )) )
{
return;
}
//***************************************************
std::ofstream metricFile;
metricFile.open( metricFilename,std::ios_base::app );
metricFile << optimizer->GetCurrentIteration() << " ";
metricFile << optimizer->GetValue() << " ";
metricFile << std::endl;
metricFile.close();
//****************************************************
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << std::endl;
}
};
// The Registration program starts here
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
// Check the number of input arguments to make sure they are correct
//******************************************************************************
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile outputImagefile ";
std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
std::cerr << " [useExplicitPDFderivatives ] [useCachingBSplineWeights ] ";
std::cerr << " [filenameForFinalTransformParameters] ";
std::cerr << " [numberOfGridNodesInsideImageInOneDimensionCoarse_0, _1 _2 ] (3 Param)";
std::cerr << " [numberOfGridNodesInsideImageInOneDimensionFine_0, _1 _2] (3 Param)";
std::cerr << " [maximumStepLength] [maximumNumberOfIterations]";
std::cerr << "[filenameForReportOfTimeAndMemoryTakenByTheRegistration]";
std::cerr << "[filenameForAffineTransform Parameters] [filenameMetricValues]";
std::cerr << " [deformationField] ";
std::cerr << std::endl;
return EXIT_FAILURE;
//**********************************************************************************
}
//***************************
metricFilename = "MetricValues.txt";
if( argc > 19 )
{
metricFilename = argv[19];
}
//***************************
const unsigned int ImageDimension = 3;
typedef signed short PixelType;
// Defining the fixed and moving image types
typedef itk::Image< PixelType, ImageDimension > FixedImageType;
typedef itk::Image< PixelType, ImageDimension > MovingImageType;
//////////////////////////////// /////////////////////
// Defining the image to be deformed. Labelmap containing Emph
//typedef itk::Image< PixelType, ImageDimension > DeformingImageType;
const unsigned int SpaceDimension = ImageDimension;
// Cubic Spline, order = 3.
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
// the rigid, the affine, and the BSpline instantiation
typedef itk::VersorRigid3DTransform< double > RigidTransformType;
typedef itk::AffineTransform< double, SpaceDimension > AffineTransformType;
typedef itk::BSplineDeformableTransform<
CoordinateRepType,
SpaceDimension,
SplineOrder > DeformableTransformType;
// rigid/affine initializer
typedef itk::CenteredTransformInitializer< RigidTransformType,
FixedImageType,
MovingImageType
> TransformInitializerType;
// Optimizer
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
//Metric
typedef itk::MattesMutualInformationImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
typedef itk:: LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
///////////////////////// /////////////////////////////////
//typedef itk:: NearestNeighborInterpolateImageFunction<
// DeformingImageType,
// double > NNInterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > RegistrationType;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
//////////////////////////// //////////////////////////////
//NNInterpolatorType::Pointer nninterpolator = NNInterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
// Connect to the registration method
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
// Auxiliary identity transform.
typedef itk::IdentityTransform<double,SpaceDimension> IdentityTransformType;
IdentityTransformType::Pointer identityTransform = IdentityTransformType::New();
//
// Read the Fixed and Moving images.
//
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
////////////////////////// ////////////////////////////
// typedef itk::ImageFileReader< DeformingImageType > DeformingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
////////////////////////////// ////////////////////////////////////////////
//DeformingImageReaderType::Pointer deformingImageReader = DeformingImageReaderType::New();
// The first argument will be the fixed and the second will be the moving image
fixedImageReader->SetFileName( argv[1] ); // EI
movingImageReader->SetFileName( argv[2] ); // 50
///////////////////////////// /////////////////////////////
// deformingImageReader->SetFileName("Hu028-LA-00.nrrd");
try
{
fixedImageReader->Update();
movingImageReader->Update();
//////////////////////////////// ////////////////////
//deformingImageReader->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
// Connecting the fixed and moving images to the registration pipeline
registration->SetFixedImage( fixedImage );
registration->SetMovingImage( movingImageReader->GetOutput() );
itkProbesCreate();
// Set Mattes Mutual information parameter
//*************************************************
metric->SetNumberOfHistogramBins( 50 * 4);
//*************************************************
FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
const unsigned int numberOfPixels = fixedRegion.GetNumberOfPixels();
// Initialize the transform
metric->ReinitializeSeed( 76926294 );
if( argc > 6 )
{
// Define whether to calculate the metric derivative by explicitly
// computing the derivatives of the joint PDF with respect to the Transform
// parameters, or doing it by progressively accumulating contributions from
// each bin in the joint PDF.
metric->SetUseExplicitPDFDerivatives( atoi( argv[6] ) );
}
//if( argc > 7 )
// {
// Define whether to cache the BSpline weights and indexes corresponding to
// each one of the samples used to compute the metric. Enabling caching will
// make the algorithm run faster but it will have a cost on the amount of memory
// that needs to be allocated. This option is only relevant when using the
// BSplineDeformableTransform.
bool b1 = false;
metric->SetUseCachingOfBSplineWeights (b1);
// }
//
// Initialize a rigid transform by using Image Intensity Moments (center of mass)
//
TransformInitializerType::Pointer initializer = TransformInitializerType::New();
RigidTransformType::Pointer rigidTransform = RigidTransformType::New();
initializer->SetTransform( rigidTransform );
initializer->SetFixedImage( fixedImageReader->GetOutput() );
initializer->SetMovingImage( movingImageReader->GetOutput() );
initializer->MomentsOn();
std::cout << "Starting Rigid Transform Initialization " << std::endl;
itkProbesStart( "Rigid Initialization");
initializer->InitializeTransform();
itkProbesStop( "Rigid Initialization");
std::cout << "Rigid Transform Initialization completed" << std::endl;
std::cout << std::endl;
registration->SetFixedImageRegion( fixedRegion );
registration->SetInitialTransformParameters( rigidTransform->GetParameters() );
registration->SetTransform( rigidTransform );
// set the scales
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales( rigidTransform->GetNumberOfParameters() );
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales( optimizerScales );
optimizer->SetMaximumStepLength( 0.2000 );
optimizer->SetMinimumStepLength( 0.0001 );
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
//show Initial Metric
//**********************************************
optimizer->SetNumberOfIterations( 1 );// was 1 H03 high res=3
metric->SetNumberOfSpatialSamples( 100000L ); // (890763) for bspline fine
std::cout << "Metric Initial value: " << std::endl;
try
{
registration->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
//********************************************
optimizer->SetNumberOfIterations( 200 );
//**************************************************************************
metric->SetNumberOfSpatialSamples( 10000L);
//********************************************************************
// doing rigid registration here:
std::cout << "Starting Rigid Registration " << std::endl;
try
{
itkProbesStart( "Rigid Registration" );
registration->Update();
itkProbesStop( "Rigid Registration" );
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// Rigid Completed
std::cout << "Rigid Registration completed" << std::endl;
std::cout << std::endl;
rigidTransform->SetParameters( registration->GetLastTransformParameters() );
std::cout << "Writing transform parameter file ...";
std::ofstream parametersFile;
parametersFile.open( "Hu03parameters6107.txt" );
parametersFile << registration->GetLastTransformParameters() << std::endl;
parametersFile.close();
std::cout << " Done!" << std::endl;
FixedImageType::Pointer im = fixedImageReader->GetOutput();
FixedImageType::RegionType region = im->GetLargestPossibleRegion();
typedef itk::ImageRegionIteratorWithIndex<FixedImageType> ImageIterator;
ImageIterator it( im, region);
typedef itk::Vector< float, ImageDimension > VectorType;
typedef itk::Image< VectorType, ImageDimension > DeformationFieldType;
DeformationFieldType::Pointer field = DeformationFieldType::New();
field->SetRegions( fixedRegion );
field->SetOrigin( fixedImage->GetOrigin() );
field->SetSpacing( fixedImage->GetSpacing() );
field->SetDirection( fixedImage->GetDirection() );
field->Allocate();
typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
FieldIterator fi( field, fixedRegion );
RigidTransformType::InputPointType fixedPoint;
RigidTransformType::OutputPointType movingPoint;
DeformationFieldType::IndexType index;
VectorType displacement;
std::ofstream DispFile;
DispFile.open( "Lung3-16-16-Crop61.txt" );
it.GoToBegin();
fi.GoToBegin();
while( ! fi.IsAtEnd() )
{
index = fi.GetIndex();
field->TransformIndexToPhysicalPoint( index, fixedPoint );
movingPoint = rigidTransform->TransformPoint( fixedPoint );
displacement = movingPoint - fixedPoint;
fi.Set( displacement );
if (it.Value() ==1)
{
// std::cout << displacement;
DispFile << displacement << std::endl;
}
++fi;
++it;
}
DispFile.close();
//////////////////////////////////////////////////////////////////////////////////////////////
///// Perform Affine Registration
AffineTransformType::Pointer affineTransform = AffineTransformType::New();
affineTransform->SetCenter( rigidTransform->GetCenter() );
affineTransform->SetTranslation( rigidTransform->GetTranslation() );
affineTransform->SetMatrix( rigidTransform->GetMatrix() );
registration->SetTransform( affineTransform );
registration->SetInitialTransformParameters( affineTransform->GetParameters() );
// set scales
optimizerScales = OptimizerScalesType( affineTransform->GetNumberOfParameters() );
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = 1.0;
optimizerScales[4] = 1.0;
optimizerScales[5] = 1.0;
optimizerScales[6] = 1.0;
optimizerScales[7] = 1.0;
optimizerScales[8] = 1.0;
optimizerScales[9] = translationScale;
optimizerScales[10] = translationScale;
optimizerScales[11] = translationScale;
optimizer->SetScales( optimizerScales );
optimizer->SetMaximumStepLength( 0.2000 );
optimizer->SetMinimumStepLength( 0.0001 );
optimizer->SetNumberOfIterations( 100 ); //50 to 100 for high res
//***********************************************************************
metric->SetNumberOfSpatialSamples( 50000L);
//***********************************************************************
std::cout << "Starting Affine Registration " << std::endl;
try
{
itkProbesStart( "Affine Registration" );
registration->Update();
itkProbesStop( "Affine Registration" );
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// affine completes here
std::cout << "Affine Registration completed" << std::endl;
std::cout << std::endl;
affineTransform->SetParameters( registration->GetLastTransformParameters() );
if( argc > 18 )
{
std::cout << "Writing Affine Transform Parameters ...";
std::ofstream affineFile;
affineFile.open( argv[18] );
affineFile << affineTransform->GetCenter()<<std::endl;
affineFile << affineTransform->GetParameters();
affineFile << std::endl;
affineFile.close();
std::cout << " Done!" << std::endl;
}
// Perform Deformable Registration
//
DeformableTransformType::Pointer bsplineTransformCoarse = DeformableTransformType::New();
//*********************************************************
//unsigned int numberOfGridNodesInOneDimensionCoarse = 5;
unsigned int numberOfGridNodesInOneDimensionCoarse[3];
// different in each dimension
numberOfGridNodesInOneDimensionCoarse[0] = 20;
numberOfGridNodesInOneDimensionCoarse[1] = 20;
numberOfGridNodesInOneDimensionCoarse[2] = 20;
// check if user defined...
if( argc > 9 )
{
numberOfGridNodesInOneDimensionCoarse[0] = atoi( argv[9] );
}
if( argc > 10 )
{
numberOfGridNodesInOneDimensionCoarse[1] = atoi( argv[10] );
}
if( argc > 11 )
{
numberOfGridNodesInOneDimensionCoarse[2] = atoi( argv[11] );
}
//************************************************************
typedef DeformableTransformType::RegionType RegionType;
RegionType bsplineRegion;
RegionType::SizeType gridSizeOnImage;
RegionType::SizeType gridBorderSize;
RegionType::SizeType totalGridSize;
//*******************************************************************
//gridSizeOnImage.Fill( numberOfGridNodesInOneDimensionCoarse );
gridSizeOnImage.SetElement(0, numberOfGridNodesInOneDimensionCoarse[0]);
gridSizeOnImage.SetElement(1, numberOfGridNodesInOneDimensionCoarse[1]);
gridSizeOnImage.SetElement(2, numberOfGridNodesInOneDimensionCoarse[2]);
//*******************************************************************
gridBorderSize.Fill( SplineOrder ); // Border for spline order = 3 ( 1 lower, 2 upper )
totalGridSize = gridSizeOnImage + gridBorderSize;
bsplineRegion.SetSize( totalGridSize );
typedef DeformableTransformType::SpacingType SpacingType;
SpacingType spacing = fixedImage->GetSpacing();
typedef DeformableTransformType::OriginType OriginType;
OriginType origin = fixedImage->GetOrigin();;
FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
for(unsigned int r=0; r<ImageDimension; r++)
{
spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /
static_cast<double>(gridSizeOnImage[r] - 1);
}
FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
SpacingType gridOriginOffset = gridDirection * spacing;
OriginType gridOrigin = origin - gridOriginOffset;
bsplineTransformCoarse->SetGridSpacing( spacing );
bsplineTransformCoarse->SetGridOrigin( gridOrigin );
bsplineTransformCoarse->SetGridRegion( bsplineRegion );
bsplineTransformCoarse->SetGridDirection( gridDirection );
bsplineTransformCoarse->SetBulkTransform( affineTransform );
// set parameters of BSpline here
typedef DeformableTransformType::ParametersType ParametersType;
unsigned int numberOfBSplineParameters = bsplineTransformCoarse->GetNumberOfParameters();
optimizerScales = OptimizerScalesType( numberOfBSplineParameters );
optimizerScales.Fill( 1.0 );
optimizer->SetScales( optimizerScales );
ParametersType initialDeformableTransformParameters( numberOfBSplineParameters );
initialDeformableTransformParameters.Fill( 0.0 );
bsplineTransformCoarse->SetParameters( initialDeformableTransformParameters );
registration->SetInitialTransformParameters( bsplineTransformCoarse->GetParameters() );
registration->SetTransform( bsplineTransformCoarse );
optimizer->SetMaximumStepLength( 10.0 );
optimizer->SetMinimumStepLength( 0.01 );
optimizer->SetRelaxationFactor( 0.7 );
optimizer->SetNumberOfIterations( 70 );// 50 to 100 for high res
// Optionally, get the step length from the command line arguments
if( argc > 15 )
{
optimizer->SetMaximumStepLength( atof( argv[15] ) );
}
// Optionally, get the number of iterations from the command line arguments
if( argc > 16 )
{
optimizer->SetNumberOfIterations( atoi( argv[16] ) );
}
//**********************************************************************
metric->SetNumberOfSpatialSamples( numberOfBSplineParameters * 100);
//**********************************************************************
// Starting from a coarse grid
std::cout << std::endl << "Starting Deformable Registration Coarse Grid" << std::endl;
try
{
itkProbesStart( "Deformable Registration Coarse" );
registration->Update();
itkProbesStop( "Deformable Registration Coarse" );
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
std::cout << "Deformable Registration Coarse Grid completed" << std::endl;
std::cout << std::endl;
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
bsplineTransformCoarse->SetParameters( finalParameters );
DeformableTransformType::Pointer bsplineTransformFine = DeformableTransformType::New();
//*********************************************************
//unsigned int numberOfGridNodesInOneDimensionFine = 20;
unsigned int numberOfGridNodesInOneDimensionFine[3];
// Parya used 50, 50, 50!
numberOfGridNodesInOneDimensionFine[0] = 50;
numberOfGridNodesInOneDimensionFine[1] = 50;
numberOfGridNodesInOneDimensionFine[2] = 35;
if( argc > 12 )
{
numberOfGridNodesInOneDimensionFine[0] = atoi( argv[12] );
}
if( argc > 13 )
{
numberOfGridNodesInOneDimensionFine[1] = atoi( argv[13] );
}
if( argc > 14 )
{
numberOfGridNodesInOneDimensionFine[2] = atoi( argv[14] );
}
//************************************************************
RegionType::SizeType gridHighSizeOnImage;
//*******************************************************************
//gridHighSizeOnImage.Fill( numberOfGridNodesInOneDimensionFine );
gridHighSizeOnImage.SetElement(0, numberOfGridNodesInOneDimensionFine[0]);
gridHighSizeOnImage.SetElement(1, numberOfGridNodesInOneDimensionFine[1]);
gridHighSizeOnImage.SetElement(2, numberOfGridNodesInOneDimensionFine[2]);
//*******************************************************************
totalGridSize = gridHighSizeOnImage + gridBorderSize; // --> total number of parameters = ((20+3)^3)*3 = 36501
bsplineRegion.SetSize( totalGridSize );
SpacingType spacingHigh = fixedImage->GetSpacing();
OriginType originHigh = fixedImage->GetOrigin();;
for(unsigned int rh=0; rh<ImageDimension; rh++)
{
spacingHigh[rh] *= static_cast<double>(fixedImageSize[rh] - 1) /
static_cast<double>(gridHighSizeOnImage[rh] - 1);
originHigh[rh] -= spacingHigh[rh];
}
SpacingType gridOriginOffsetHigh = gridDirection * spacingHigh;
OriginType gridOriginHigh = origin - gridOriginOffsetHigh;
bsplineTransformFine->SetGridSpacing( spacingHigh );
bsplineTransformFine->SetGridOrigin( gridOriginHigh );
bsplineTransformFine->SetGridRegion( bsplineRegion );
bsplineTransformFine->SetGridDirection( gridDirection );
bsplineTransformFine->SetBulkTransform( affineTransform );
numberOfBSplineParameters = bsplineTransformFine->GetNumberOfParameters();
ParametersType parametersHigh( numberOfBSplineParameters );
parametersHigh.Fill( 0.0 );
unsigned int counter = 0;
for ( unsigned int k = 0; k < SpaceDimension; k++ )
{
typedef DeformableTransformType::ImageType ParametersImageType;
typedef itk::ResampleImageFilter<ParametersImageType,ParametersImageType> ResamplerType;
ResamplerType::Pointer upsampler = ResamplerType::New();
typedef itk::BSplineResampleImageFunction<ParametersImageType,double> FunctionType;
FunctionType::Pointer function = FunctionType::New();
upsampler->SetInput( bsplineTransformCoarse->GetCoefficientImages()[k] );
upsampler->SetInterpolator( function );
upsampler->SetTransform( identityTransform );
upsampler->SetSize( bsplineTransformFine->GetGridRegion().GetSize() );
upsampler->SetOutputSpacing( bsplineTransformFine->GetGridSpacing() );
upsampler->SetOutputOrigin( bsplineTransformFine->GetGridOrigin() );
typedef itk::BSplineDecompositionImageFilter<ParametersImageType,ParametersImageType>
DecompositionType;
DecompositionType::Pointer decomposition = DecompositionType::New();
decomposition->SetSplineOrder( SplineOrder );
decomposition->SetInput( upsampler->GetOutput() );
decomposition->Update();
ParametersImageType::Pointer newCoefficients = decomposition->GetOutput();
// copy the coefficients into the parameter array
typedef itk::ImageRegionIterator<ParametersImageType> Iterator;
Iterator it( newCoefficients, bsplineTransformFine->GetGridRegion() );
while ( !it.IsAtEnd() )
{
parametersHigh[ counter++ ] = it.Get();
++it;
}
}
optimizerScales = OptimizerScalesType( numberOfBSplineParameters );
optimizerScales.Fill( 1.0 );
optimizer->SetScales( optimizerScales );
bsplineTransformFine->SetParameters( parametersHigh );
//**************************************************************
// Define whether to cache the BSpline weights and indexes corresponding to
// each one of the samples used to compute the metric. Enabling caching will
// make the algorithm run faster but it will have a cost on the amount of memory
// that needs to be allocated. This option is only relevant when using the
// BSplineDeformableTransform.
metric->SetUseCachingOfBSplineWeights( atoi( "0" ) );
//***************************************************************
std::cout << "Starting Registration with high resolution transform" << std::endl;
registration->SetInitialTransformParameters( bsplineTransformFine->GetParameters() );
registration->SetTransform( bsplineTransformFine );
const unsigned long numberOfSamples =
static_cast<unsigned long>(
vcl_sqrt( static_cast<double>( numberOfBSplineParameters ) *
static_cast<double>( numberOfPixels ) ) );
//*******************************************************
metric->SetNumberOfSpatialSamples( numberOfSamples );
//*******************************************************
try
{
itkProbesStart( "Deformable Registration Fine" );
registration->Update();
itkProbesStop( "Deformable Registration Fine" );
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// Software Guide : EndCodeSnippet
std::cout << "Deformable Registration Fine Grid completed" << std::endl;
std::cout << std::endl;
// Report the time and memory taken by the registration
itkProbesReport( std::cout );
finalParameters = registration->GetLastTransformParameters();
bsplineTransformFine->SetParameters( finalParameters );
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( bsplineTransformFine);
resample->SetInput( movingImageReader->GetOutput() );
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetOutputDirection( fixedImage->GetDirection() );
//****************************************************
resample->SetInterpolator(interpolator);
//*****************************************************
//**************************************************************************
resample->SetDefaultPixelValue( -1000 ); // I changed it based on approximation on CT images
//***************************************************************************
typedef signed short OutputPixelType;
typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
////////////////////////////////////// /////////////////////////////////////
// typedef itk::ResampleImageFilter<
// DeformingImageType,
// FixedImageType > DefResampleFilterType;
// DefResampleFilterType::Pointer defresample = ResampleFilterType::New();
// defresample->SetTransform(bsplineTransformFine);
// defresample->SetInput(deformingImageReader->GetOutput());
// defresample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
// defresample->SetOutputOrigin(fixedImage->GetOrigin());
// defresample->SetOutputSpacing(fixedImage->GetSpacing());
// defresample->SetOutputDirection(fixedImage->GetDirection());
//****************************************************
// defresample->SetInterpolator(nninterpolator);
//*****************************************************
//**************************************************************************
// defresample->SetDefaultPixelValue(0); // I changed it for the emphysema Labelmap!
//***************************************************************************
// typedef signed short OutputPixelType;
// typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
// typedef itk::CastImageFilter<
// FixedImageType,
// OutputImageType > CastFilterType;
// typedef itk::ImageFileWriter< OutputImageType > DefWriterType;
///////////////////////////////////////////// /////////////////////////
typedef signed short DiffPixelType;
typedef itk::Image< DiffPixelType, ImageDimension > DiffImageType;
//*****************************************************************
/*typedef itk::SquaredDifferenceImageFilter<
FixedImageType,
FixedImageType,
DiffImageType > DifferenceFilterType;*/
typedef itk::AbsoluteValueDifferenceImageFilter<
FixedImageType,
FixedImageType,
DiffImageType > DifferenceFilterType;
//********************************************************************
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
typedef itk::ImageFileWriter< DiffImageType > DiffWriterType;
DiffWriterType::Pointer writer2 = DiffWriterType::New();
writer2->SetInput( difference->GetOutput() );
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[3] ); //
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
std::cout << "Writing resampled moving image...";
try
{
writer->Update();
}
catch (itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
std::cout << " Done!" << std::endl;
/////////////////////////////// ///////////////////////////////////
//////////////////////Write the resampled Emph def here/////////////////////////
//WriterType::Pointer defwriter = DefWriterType::New();
// CastFilterType::Pointer caster = CastFilterType::New();
// defwriter->SetFileName("Hu028-LA-00-defto-50.nrrd"); //
// caster->SetInput(defresample->GetOutput());
// defwriter->SetInput(caster->GetOutput());
// std::cout << "Writing resampled emphysema deforming image...";
// try
// {
// defwriter->Update();
// }
// catch (itk::ExceptionObject & err)
// {
// std::cerr << "ExceptionObject caught !" << std::endl;
// std::cerr << err << std::endl;
// return EXIT_FAILURE;
// }
// std::cout << " Done!" << std::endl;
/////////////////////////////////////////////////////
////////////////////////////////////////////////////
//// Compute the difference image between the
//// fixed and resampled moving image.
if( argc > 4 )
{
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resample->GetOutput() );
writer2->SetFileName( argv[4] );
std::cout << "Writing difference image after registration...";
try
{
writer2->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
std::cout << " Done!" << std::endl;
}
//// Compute the difference image between the
//// fixed and moving image before registration.
if( argc > 5 )
{
writer2->SetFileName( argv[5] );
difference->SetInput1( fixedImageReader->GetOutput() );
resample->SetTransform( identityTransform );
std::cout << "Writing difference image before registration...";
try
{
writer2->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
std::cout << " Done!" << std::endl;
}
// Generate the explicit deformation field resulting from
// the registration.
if( argc > 6 )
{
typedef itk::Vector< float, ImageDimension > VectorType;
typedef itk::Image< VectorType, ImageDimension > DeformationFieldType;
DeformationFieldType::Pointer field = DeformationFieldType::New();
field->SetRegions( fixedRegion );
field->SetOrigin( fixedImage->GetOrigin() );
field->SetSpacing( fixedImage->GetSpacing() );
field->SetDirection( fixedImage->GetDirection() );
field->Allocate();
typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
FieldIterator fi( field, fixedRegion );
fi.GoToBegin();
DeformableTransformType::InputPointType fixedPoint;
DeformableTransformType::OutputPointType movingPoint;
DeformationFieldType::IndexType index;
VectorType displacement;
while( ! fi.IsAtEnd() )
{
index = fi.GetIndex();
field->TransformIndexToPhysicalPoint( index, fixedPoint );
movingPoint = bsplineTransformFine->TransformPoint( fixedPoint );
displacement = movingPoint - fixedPoint;
fi.Set( displacement );
++fi;
}
typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
fieldWriter->SetInput( field );
fieldWriter->SetFileName( argv[6] );
std::cout << "Writing deformation field ...";
try
{
fieldWriter->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
std::cout << " Done!" << std::endl;
}
// Optionally, save the transform parameters in a file
if( argc > 7 )
{
std::cout << "Writing transform parameter file ...";
std::ofstream parametersFile;
parametersFile.open( argv[7] );
parametersFile << finalParameters << std::endl;
parametersFile.close();
std::cout << " Done!" << std::endl;
}
// Optionally, save the time and memory taken by the registration
if( argc > 17 )
{
std::cout << "Writing the time and memory taken by the registration ...";
std::ofstream timeFile;
timeFile.open( argv[17] );
itkProbesReport( timeFile );
timeFile << std::endl;
timeFile.close();
std::cout << " Done!" << std::endl;
}
return EXIT_SUCCESS;
}
#undef itkProbesCreate
#undef itkProbesStart
#undef itkProbesStop
#undef itkProbesReport