[ITK Community] [Insight-users] Bspline artefact multi resolution - origin offset - Bug?
Emma Saunders
emmasaunders123 at gmail.com
Wed Dec 11 17:26:16 EST 2013
Hi everyone
I appear to be getting a step artefact at the boundaries of my registered
image, using a multi-resolution Bspline approach. The distance of the
artefact co-indices with the origin of the original image. I can't seem to
find why this occurs. I set the Bspline grid wrt the original fixed image
so unsure why this is happening.
I would really appreciate any advice:
Below is the code if anyone is interested:
Many Thanks
Emma
#include "itkMedianImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkImageRegistrationMethod.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkAffineTransform.h"
#include "itkBSplineDeformableTransform.h" //version 3
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkMemoryProbesCollectorBase.h"
#include <itkHistogramMatchingImageFilter.h>
#include "itkBSplineResampleImageFunction.h"
#include "itkLBFGSBOptimizer.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLBFGSOptimizer.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkBSplineResampleImageFunction.h"
#include "itkIdentityTransform.h"
#include "itkBSplineDecompositionImageFilter.h"
#include "itkCommand.h"
#include <fstream>
#include "itkCommand.h"
unsigned int Counter = 0;
unsigned int Countold = 0;
//This is for the affine registration
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::cout << optimizer->GetValue() << " " <<std::endl;
}
};
//This is for B spline
class BCommandIterationUpdate : public itk::Command
{
public:
typedef BCommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<BCommandIterationUpdate> Pointer;
itkNewMacro( BCommandIterationUpdate );
protected:
BCommandIterationUpdate() {};
public:
typedef itk::LBFGSBOptimizer OptimizerType;
typedef const OptimizerType * OptimizerPointer;
public:
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)
{
const OptimizerPointer optimizer =
dynamic_cast< OptimizerPointer >( object );
if( !(itk::IterationEvent().CheckEvent( &event )) )
{
return;
}
Countold=Counter;
Counter = optimizer->GetCurrentIteration();
if (Counter!=Countold)
{
std::cout << optimizer->GetValue() << " " <<std::endl;
}
}
};
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile " << std::endl;
std::cerr << " outputImagefile" << std::endl;
std::cerr << " outputDeformationField " << std::endl;
return EXIT_FAILURE;
}
// READ AND CREATE THE IMAGES AND DISPLACEMENT FIELD
const unsigned int Dimension = 3;
typedef float PixelType;
typedef itk::Vector< double, 3 > VectorPixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
typedef itk::Image< PixelType, Dimension > OutputImageType;
typedef itk::Image< VectorPixelType, 3 > DisplacementFieldType;
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
typedef itk::ImageFileWriter< OutputImageType > ImageWriterType;
typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType;
typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType;
FieldWriterType::Pointer FieldWriter = FieldWriterType::New();
FixedImageReaderType::Pointer FixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer MovingImageReader =
MovingImageReaderType::New();
ImageWriterType::Pointer ImageWriter = ImageWriterType::New();
FixedImageReader->SetFileName( argv[1] );
MovingImageReader->SetFileName( argv[2] );
ImageWriter->SetFileName( argv[3] );
FieldWriter->SetFileName( argv[4] );
FixedImageType::Pointer FixedImage = FixedImageReader->GetOutput();
MovingImageType::Pointer MovedImage = MovingImageReader->GetOutput();
OutputImageType::Pointer OutImage = OutputImageType::New();
DisplacementFieldType::Pointer DeformField = DisplacementFieldType::New();
std::cout << "First Fixed Image origin is " << FixedImage->GetOrigin() <<
std::endl ;
FixedImage->Update();
MovedImage->Update();
//Median and Gaussian Filter
typedef itk::MedianImageFilter<FixedImageType, FixedImageType > FilterType;
// Create and setup a median filter
FilterType::Pointer medianFilterFixed = FilterType::New();
FilterType::InputSizeType radius;
radius.Fill(1);
medianFilterFixed->SetRadius(radius);
medianFilterFixed->SetInput( FixedImage );
medianFilterFixed->Update();
FilterType::Pointer medianFilterMoved = FilterType::New();
medianFilterMoved->SetRadius(radius);
medianFilterMoved->SetInput( MovedImage );
medianFilterMoved->Update();
typedef itk::DiscreteGaussianImageFilter<
FixedImageType,
FixedImageType
> GaussianFilterType;
GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
fixedSmoother->SetVariance(1.0 );
movingSmoother->SetVariance(1.0 );
fixedSmoother->SetInput( medianFilterFixed->GetOutput() );
movingSmoother->SetInput ( medianFilterMoved->GetOutput() );
fixedSmoother->Update();
movingSmoother->Update();
std::cout<< "The fixed origin is " <<
fixedSmoother->GetOutput()->GetOrigin() <<std::endl;
std::cout << "Fixed Image origin is " << FixedImage->GetOrigin() <<
std::endl ;
//Define the transformations to be used
//PERFROM AN INITIAL AFFINE
typedef itk::AffineTransform<
double,
Dimension > AffineTransformType;
//This will be followed by Bspline
const unsigned int SpaceDimension = Dimension;
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
typedef itk::BSplineDeformableTransform<
CoordinateRepType,
SpaceDimension,
SplineOrder > DeformableTransformType;
//version 3
//Now for the Optimizer, Metric and Interpolatror and Registration
//which are connected to the registration
//typedef itk::LBFGSOptimizer bOptimizerType;
typedef itk::LBFGSBOptimizer bOptimizerType;
bOptimizerType::Pointer boptimizer = bOptimizerType::New();
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef itk::MattesMutualInformationImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
typedef itk:: LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > RegistrationType;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
// Setup the metric parameters
metric->ReinitializeSeed( 76926294 );
metric->SetUseCachingOfBSplineWeights(1);
const unsigned long numberOfImagePixels =
FixedImage->GetLargestPossibleRegion().GetNumberOfPixels();
const unsigned long numberOfSpatialSamplesAffine =
static_cast<unsigned long>( numberOfImagePixels * 1 );
metric->SetNumberOfSpatialSamples( numberOfSpatialSamplesAffine );
metric->SetNumberOfHistogramBins( 40 ); // Reduce the number of bins if a
deformable registration fails. If the number of bins is too large, the
estimated PDFs will be a field of impulses and will inhibit reliable
registration estimation.
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
//create the Affine transform
AffineTransformType::Pointer Affinetransform =
AffineTransformType::New();
registration->SetTransform( Affinetransform );
//input the images
registration->SetFixedImage(fixedSmoother->GetOutput());
registration->SetMovingImage(movingSmoother->GetOutput());
//input the region
registration->SetFixedImageRegion(FixedImage->GetBufferedRegion() );
//set the initial parameters
typedef RegistrationType::ParametersType AffineParametersType;
AffineParametersType initialParameters(
Affinetransform->GetNumberOfParameters() );
Affinetransform->SetIdentity();
registration->SetInitialTransformParameters(
Affinetransform->GetParameters() );
//set scale values for optimizer to search sensibly and set up optimizer
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales(
Affinetransform->GetNumberOfParameters() );
FixedImageType:: SpacingType fixedspacing = FixedImage->GetSpacing();
FixedImageType::RegionType region=FixedImage->GetLargestPossibleRegion();
const unsigned int numberOfPixels = region.GetNumberOfPixels();
FixedImageType:: SizeType fixedsize =region.GetSize();
double translationScalea = 1.0 /( 10.0 * fixedspacing[0] * fixedsize[0]
);
double translationScaleb = 1.0 /( 10.0 * fixedspacing[1] * fixedsize[1]
);
double translationScalec = 1.0 /( 10.0 * fixedspacing[2] * fixedsize[2]
);
//for 3D use below
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] = translationScalea;
optimizerScales[10] = translationScaleb;
optimizerScales[11] = translationScalec;
//step length was 0.1 larger step length quicker converge but may get
erreaic values in the metric as it decreases
unsigned int maxNumberOfIterations = 60;
optimizer->SetMaximumStepLength( 0.05 );
optimizer->SetMinimumStepLength( 0.001 ); //was 0.00001
optimizer->SetNumberOfIterations( maxNumberOfIterations );
optimizer->MaximizeOff();
//Perform Affine Registration put observer on the Affine optimizer
CommandIterationUpdate::Pointer observera = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observera );
try
{
registration->Update();
std::cout << " " <<std::endl;
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
std::cout << "Affine Registration completed" << std::endl;
std::cout << std::endl;
optimizer->RemoveAllObservers();
optimizer->RemoveAllObservers();
//Obtain the last Affine transform which will be used as a bulk transform
in the Bspline
Affinetransform->SetParameters(registration->GetLastTransformParameters()
);
//////////////////////////////////////////////////////////////////////////////////////////////
///NOW FOR THE Bspline Registration
// We construct two transform objects, each one will be configured for a
resolution level.
// Notice than in this multi-resolution scheme we are not modifying the
// resolution of the image, but rather the flexibility of the deformable
// transform itself.
//SO must define a transform and a grid for each level
DeformableTransformType::Pointer BtransformLow =
DeformableTransformType::New();
typedef DeformableTransformType::RegionType RegionType;
RegionType bsplineRegionlow;
RegionType::SizeType gridBorderSize;
RegionType::SizeType gridSizeOnImagelow;
RegionType::SizeType totalGridSizelow;
//.//////////////////////////////////////////////////////////////////////
///GRID SIZE FOR LOW RESOLUTION
/////////////////////////////////////////////////////////////////////////
gridSizeOnImagelow[0]=12; //lat coronal, sup-inf //want 5 2 7 10 4 15
gridSizeOnImagelow[1]=4; // want
gridSizeOnImagelow[2]=12;
gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2
upper )
totalGridSizelow = gridSizeOnImagelow + gridBorderSize;
bsplineRegionlow.SetSize( totalGridSizelow );
typedef DeformableTransformType::SpacingType SpacingTypelow;
SpacingTypelow spacinglow = FixedImage->GetSpacing();
typedef DeformableTransformType::OriginType OriginTypelow;
OriginTypelow originlow = FixedImage->GetOrigin();
FixedImageType::RegionType FixedRegionlow =
FixedImage->GetBufferedRegion();
FixedImageType::SpacingType fixedSpacinglow =
FixedImage->GetSpacing();
FixedImageType::SizeType FixedImageSizelow = FixedRegionlow.GetSize();
for(unsigned int r=0; r<Dimension; r++)
{
spacinglow[r] = floor( fixedSpacinglow[r] * (FixedImageSizelow[r] - 1)
/ (gridSizeOnImagelow[r] - 1) ); //new way
}
FixedImageType::DirectionType gridDirection = FixedImage->GetDirection();
SpacingTypelow gridOriginOffsetlow = gridDirection * spacinglow;
OriginTypelow gridOriginlow = originlow - gridOriginOffsetlow;
BtransformLow->SetGridSpacing( spacinglow );
BtransformLow->SetGridOrigin( gridOriginlow );
BtransformLow->SetGridRegion( bsplineRegionlow ); //i.e region of gridsize
BtransformLow->SetGridDirection( gridDirection );
std::cout << "grid region for low Bspline is " <<
BtransformLow->GetGridRegion() <<std::endl;
//////////////////////////////////////////////////////////////////////////////
///GRID FOR MEDIUM
////////////////////////////////////////////////////////////////////////////////
DeformableTransformType::Pointer BtransformMed=
DeformableTransformType::New();
RegionType bsplineRegionMed;
RegionType::SizeType gridSizeOnImageMed;
RegionType::SizeType totalGridSizeMed;
//.//////////////////////////////////////////////////////////////////////
///GRID SIZE FOR MEDIUM RESOLUTION
/////////////////////////////////////////////////////////////////////////
gridSizeOnImageMed[0]=24; //lat coronal, sup-inf
gridSizeOnImageMed[1]=8; // want
gridSizeOnImageMed[2]=24;
totalGridSizeMed = gridSizeOnImageMed + gridBorderSize;
bsplineRegionMed.SetSize( totalGridSizeMed );
typedef DeformableTransformType::SpacingType SpacingTypeMed;
SpacingTypeMed spacingMed = FixedImage->GetSpacing();
typedef DeformableTransformType::OriginType OriginTypeMed;
OriginTypeMed originMed = FixedImage->GetOrigin();
FixedImageType::RegionType FixedRegionMed =
FixedImage->GetBufferedRegion();
FixedImageType::SpacingType fixedSpacingMed =
FixedImage->GetSpacing();
FixedImageType::SizeType FixedImageSizeMed = FixedRegionMed.GetSize();
for(unsigned int r=0; r<Dimension; r++)
{
spacingMed[r] = floor( fixedSpacingMed[r] * (FixedImageSizeMed[r] - 1)
/ (gridSizeOnImageMed[r] - 1) ); //new way
}
SpacingTypeMed gridOriginOffsetMed = gridDirection * spacingMed;
OriginTypeMed gridOriginMed = originMed - gridOriginOffsetMed;
BtransformMed->SetGridSpacing( spacingMed );
BtransformMed->SetGridOrigin( gridOriginMed );
BtransformMed->SetGridRegion( bsplineRegionMed ); //i.e region of gridsize
BtransformMed->SetGridDirection( gridDirection );
std::cout << "grid region for Med Bspline is " <<
BtransformMed->GetGridRegion() <<std::endl;
//////////////////////////////////////////////////////////////////////////////
///GRID FOR HIGH
////////////////////////////////////////////////////////////////////////////////
//Now the Higher Transform
DeformableTransformType::Pointer BtransformHigh =
DeformableTransformType::New();
RegionType bsplineRegionhigh;
RegionType::SizeType gridSizeOnImagehigh;
RegionType::SizeType totalGridSizehigh;
//.//////////////////////////////////////////////////////////////////////
///GRID SIZE FOR HIGH RESOLUTION
/////////////////////////////////////////////////////////////////////////
gridSizeOnImagehigh[0]=50; //lat coronal, sup-inf
gridSizeOnImagehigh[1]=16; // want
gridSizeOnImagehigh[2]=50;
totalGridSizehigh = gridSizeOnImagehigh + gridBorderSize;
bsplineRegionhigh.SetSize( totalGridSizehigh );
typedef DeformableTransformType::SpacingType SpacingTypehigh;
SpacingTypehigh spacinghigh = FixedImage->GetSpacing();
typedef DeformableTransformType::OriginType OriginTypehigh;
OriginTypehigh originhigh = FixedImage->GetOrigin();
FixedImageType::RegionType FixedRegionhigh =
FixedImage->GetBufferedRegion();
FixedImageType::SpacingType fixedSpacinghigh =
FixedImage->GetSpacing();
FixedImageType::SizeType FixedImageSizehigh = FixedRegionhigh.GetSize();
for(unsigned int r=0; r<Dimension; r++)
{
spacinghigh[r] = floor( fixedSpacinghigh[r] * (FixedImageSizehigh[r] -
1) / (gridSizeOnImagehigh[r] - 1) ); //new way
}
SpacingTypehigh gridOriginOffsethigh = gridDirection * spacinghigh;
OriginTypehigh gridOriginhigh = originhigh - gridOriginOffsethigh;
BtransformHigh->SetGridSpacing( spacinghigh );
BtransformHigh->SetGridOrigin( gridOriginhigh );
BtransformHigh->SetGridRegion( bsplineRegionhigh ); //i.e region of
gridsize
BtransformHigh->SetGridDirection( gridDirection );
std::cout << "grid region for high Bspline is " <<
BtransformHigh->GetGridRegion() <<std::endl;
///////////////////////////////////////////////////////////////////////
///METRIC NUMBER OF SAMPLES
//////////////////////////////////////////////////////////////////////
//Low
const unsigned int numberOfBSplineParametersLow =
BtransformLow->GetNumberOfParameters();
typedef DeformableTransformType::ParametersType ParametersTypeLow;
ParametersTypeLow parametersLow( numberOfBSplineParametersLow );
//Medium
const unsigned int numberOfBSplineParametersMed =
BtransformMed->GetNumberOfParameters();
typedef DeformableTransformType::ParametersType ParametersTypeMed;
ParametersTypeMed parametersMed( numberOfBSplineParametersMed );
const unsigned int numberOfBSplineParametersHigh =
BtransformHigh->GetNumberOfParameters();
typedef DeformableTransformType::ParametersType ParametersTypeHigh;
ParametersTypeHigh parametersHigh( numberOfBSplineParametersHigh );
const unsigned long numberOfSamplesBsplineHigh =
static_cast<unsigned long>(
vcl_sqrt( static_cast<double>( numberOfBSplineParametersHigh ) *
static_cast<double>( numberOfPixels ) ) );
///////////////////////////////////////////////////////////////////////
///LOW PARAMETERS
//////////////////////////////////////////////////////////////////////
typedef DeformableTransformType::ParametersType ParametersTypelow;
const unsigned int numberOfBSplineParameterslow =
BtransformLow->GetNumberOfParameters();
ParametersTypelow parameterslow( numberOfBSplineParameterslow );
parameterslow.Fill( 0.0 );
BtransformLow->SetBulkTransform( Affinetransform );
//////////////////////////////////////////////////////////////////////////////
///OPTIMIZER FOR LOW RESOLUTION
////////////////////////////////////////////////////////////////////////////////
// if get IFLAG= -1. LINE SEARCH FAILED
// increase step lentgth but keep other parameters the same
/////////////////////////////////////////////////////////////////////
// SEE HERE FOR OPTIMIZER ADVICE
//http://www.itk.org/pipermail/insight-users/2011-March/040286.html
///FOR the LBFGSBOptimizer got LOW
bOptimizerType::BoundSelectionType boundSelectlow(
BtransformLow->GetNumberOfParameters() );
bOptimizerType::BoundValueType upperBoundlow(
BtransformLow->GetNumberOfParameters() );
bOptimizerType::BoundValueType lowerBoundlow(
BtransformLow->GetNumberOfParameters() );
boundSelectlow.Fill( 0 );
upperBoundlow.Fill( 0.0 );
lowerBoundlow.Fill( 0.0 );
boptimizer->SetBoundSelection( boundSelectlow );
boptimizer->SetUpperBound( upperBoundlow );
boptimizer->SetLowerBound( lowerBoundlow );
boptimizer->SetCostFunctionConvergenceFactor( 1.e1 );
boptimizer->SetProjectedGradientTolerance( 1e-6 );
boptimizer->SetMaximumNumberOfIterations( 200 ); //was 200
boptimizer->SetMaximumNumberOfEvaluations( 300 ); //
boptimizer->SetMaximumNumberOfCorrections( 5 );
///////////////////////////////////////////////////////////////
metric->SetNumberOfSpatialSamples( numberOfPixels*1 ); //
RegistrationType::Pointer bregistration = RegistrationType::New();
bregistration->SetNumberOfThreads(1);
bregistration->SetInitialTransformParameters(
BtransformLow->GetParameters() );
bregistration->SetTransform(BtransformLow);
bregistration->SetFixedImageRegion(FixedImage->GetBufferedRegion() );
bregistration->SetMetric( metric );
bregistration->SetOptimizer( boptimizer );
bregistration->SetInterpolator( interpolator );
bregistration->SetFixedImage(fixedSmoother->GetOutput());
bregistration->SetMovingImage(movingSmoother->GetOutput());
bregistration->SetNumberOfThreads(1);
//////////////////////////////////////////////////////////////////////
///REGISTRATION FOR LOW
///////////////////////////////////////////////////////////////////////
BCommandIterationUpdate::Pointer observerl =
BCommandIterationUpdate::New();
boptimizer->AddObserver( itk::IterationEvent(), observerl );
try
{
bregistration->Update();
std::cout << " " <<std::endl;
std::cout << "Optimizer stop condition: "
<< bregistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
boptimizer->RemoveAllObservers();
/////////////////////////////////////////////////////////////////////
///MEDIUM PARAMETERS
////////////////////////////////////////////////////////////////////
typedef DeformableTransformType::ParametersType ParametersTypeMed;
parametersMed.Fill( 0.0 );
// Now we need to initialize the BSpline coefficients of the higher
resolution
// transform. This is done by first computing the actual deformation field
// at the higher resolution from the lower resolution BSpline
coefficients.
// Then a BSpline decomposition is done to obtain the BSpline coefficient
of
// the higher resolution transform.
unsigned int counterMed = 0;
for ( unsigned int k = 0; k < SpaceDimension; k++ )
{
typedef DeformableTransformType::ImageType ParametersImageTypeMed;
typedef
itk::ResampleImageFilter<ParametersImageTypeMed,ParametersImageTypeMed>
ResamplerTypeMed;
ResamplerTypeMed::Pointer upsamplerMed = ResamplerTypeMed::New();
typedef
itk::BSplineResampleImageFunction<ParametersImageTypeMed,double>
FunctionTypeMed;
FunctionTypeMed::Pointer functionMed = FunctionTypeMed::New();
typedef itk::IdentityTransform<double,SpaceDimension>
IdentityTransformTypeMed;
IdentityTransformTypeMed::Pointer identityMed =
IdentityTransformTypeMed::New();
upsamplerMed->SetInput( BtransformLow->GetCoefficientImages()[k] );
upsamplerMed->SetInterpolator( functionMed );
upsamplerMed->SetTransform( identityMed );
upsamplerMed->SetSize( BtransformMed->GetGridRegion().GetSize() );
upsamplerMed->SetOutputSpacing( BtransformMed->GetGridSpacing() );
upsamplerMed->SetOutputOrigin( BtransformMed->GetGridOrigin() );
upsamplerMed->SetOutputDirection( FixedImage->GetDirection() );
upsamplerMed->Update();
typedef
itk::BSplineDecompositionImageFilter<ParametersImageTypeMed,ParametersImageTypeMed>
DecompositionTypeMed;
DecompositionTypeMed::Pointer decompositionMed =
DecompositionTypeMed::New();
decompositionMed->SetSplineOrder( SplineOrder );
decompositionMed->SetInput( upsamplerMed->GetOutput() );
decompositionMed->Update();
ParametersImageTypeMed::Pointer newCoefficientsMed =
decompositionMed->GetOutput();
// copy the coefficients into the parameter array
typedef itk::ImageRegionIterator<ParametersImageTypeMed> IteratorMed;
IteratorMed it( newCoefficientsMed, BtransformMed->GetGridRegion() );
while ( !it.IsAtEnd() )
{
parametersMed[ counterMed++ ] = it.Get();
++it;
}
}
BtransformMed->SetParameters( parametersMed );
////////////////////////////////////////////////////////////////
///FOR the LBFGSBOptimizer MEDIUM
bOptimizerType::BoundSelectionType boundSelectMed(
BtransformMed->GetNumberOfParameters() );
bOptimizerType::BoundValueType upperBoundMed(
BtransformMed->GetNumberOfParameters() );
bOptimizerType::BoundValueType lowerBoundMed(
BtransformMed->GetNumberOfParameters() );
boundSelectMed.Fill( 0 );
upperBoundMed.Fill( 0.0 );
lowerBoundMed.Fill( 0.0 );
boptimizer->SetBoundSelection( boundSelectMed );
boptimizer->SetUpperBound( upperBoundMed );
boptimizer->SetLowerBound( lowerBoundMed );
boptimizer->SetCostFunctionConvergenceFactor( 1.e12 ); //was 1.e7
//set/Get the CostFunctionConvergenceFactor. Algorithm terminates
// when the reduction in cost function is less than factor * epsmcj
// where epsmch is the machine precision.
// Typical values for factor: 1e+12 for low accuracy;
// 1e+7 for moderate accuracy and 1e+1 for extremely high accuracy.
boptimizer->SetProjectedGradientTolerance( 1e-6 );
//Gradient tolerance determines a lower threshold which cuases the
optimization to stop if its
//value reached, recommended values for gradient tolerance are 1e-6 to
1e-10, higher values might result in premature elimination of the
optimization process
boptimizer->SetMaximumNumberOfIterations( 200); //was 200
boptimizer->SetMaximumNumberOfEvaluations( 300 );
boptimizer->SetMaximumNumberOfCorrections( 5 );
//M Number of correction- number of function/gradient pairs used to build
model
//This parameter is passed to the function which is used to create
optimizer object.
// Increase of number of corrections decreases number of function
evaluations
//and iterations needed to converge, but increases computational overhead
associated with iteration.
// On well-conditioned problems can be as small as 3-10.
//If function changes rapidly in some directions and slowly in other
ones,
//then you can try increasing M in order to increase convergenc
///////////////////////////////////////////////////////////////////////
///REGISTRATION FOR MEDIUM
/////////////////////////////////////////////////////////////////////
//For time probe
metric->SetNumberOfSpatialSamples( numberOfPixels*1 ); //
bregistration->SetInitialTransformParameters(
BtransformMed->GetParameters() );
bregistration->SetTransform( BtransformMed );
itk::TimeProbesCollectorBase chronometerm;
itk::MemoryProbesCollectorBase memorymeterm;
BCommandIterationUpdate::Pointer observerm =
BCommandIterationUpdate::New();
boptimizer->AddObserver( itk::IterationEvent(), observerm );
try
{
bregistration->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
boptimizer->RemoveAllObservers();
BtransformMed->SetParameters( bregistration->GetLastTransformParameters()
);
//NOW HIGH RESOLUTION
/////////////////////////////////////////////////////////////////////
///HIGH PARAMETERS
////////////////////////////////////////////////////////////////////
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();
typedef itk::IdentityTransform<double,SpaceDimension>
IdentityTransformType;
IdentityTransformType::Pointer identity = IdentityTransformType::New();
upsampler->SetInput( BtransformMed->GetCoefficientImages()[k] );
upsampler->SetInterpolator( function );
upsampler->SetTransform( identity );
upsampler->SetSize( BtransformHigh->GetGridRegion().GetSize() );
upsampler->SetOutputSpacing( BtransformHigh->GetGridSpacing() );
upsampler->SetOutputOrigin( BtransformHigh->GetGridOrigin() );
upsampler->SetOutputDirection( FixedImage->GetDirection() );
upsampler->Update();
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, BtransformHigh->GetGridRegion() );
while ( !it.IsAtEnd() )
{
parametersHigh[ counter++ ] = it.Get();
++it;
}
}
BtransformHigh->SetParameters( parametersHigh );
/////////////////////////////////////////////////////////////////////////////////////////
///FOR the LBFGSBOptimizer HIGH
bOptimizerType::BoundSelectionType boundSelecthigh(
BtransformHigh->GetNumberOfParameters() );
bOptimizerType::BoundValueType upperBoundhigh(
BtransformHigh->GetNumberOfParameters() );
bOptimizerType::BoundValueType lowerBoundhigh(
BtransformHigh->GetNumberOfParameters() );
boundSelecthigh.Fill( 0 );
upperBoundhigh.Fill( 0.0 );
lowerBoundhigh.Fill( 0.0 );
boptimizer->SetBoundSelection( boundSelecthigh );
boptimizer->SetUpperBound( upperBoundhigh );
boptimizer->SetLowerBound( lowerBoundhigh );
boptimizer->SetCostFunctionConvergenceFactor( 1.e12 );
boptimizer->SetProjectedGradientTolerance( 1e-6 );//was 1e-6
boptimizer->SetMaximumNumberOfIterations( 1000 ); //was 200
boptimizer->SetMaximumNumberOfEvaluations( 300 );
boptimizer->SetMaximumNumberOfCorrections( 5 );
///////////////////////////////////////////////////////////////////////
///REGISTRATION FOR HIGH
/////////////////////////////////////////////////////////////////////
//For time probe
metric->SetNumberOfSpatialSamples( numberOfPixels*1 ); //
bregistration->SetInitialTransformParameters(
BtransformHigh->GetParameters() );
bregistration->SetTransform( BtransformHigh );
itk::TimeProbesCollectorBase chronometerh;
itk::MemoryProbesCollectorBase memorymeterh;
BCommandIterationUpdate::Pointer observerh =
BCommandIterationUpdate::New();
boptimizer->AddObserver( itk::IterationEvent(), observerh );
try
{
bregistration->Update();
std::cout << " " <<std::endl;
std::cout << "Optimizer stop condition: "
<< bregistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// Finally we use the last transform parameters in order to resample the
image.
//
boptimizer->RemoveAllObservers();
BtransformHigh->SetParameters( bregistration->GetLastTransformParameters()
);
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( BtransformHigh );
resample->SetInput( MovedImage );
resample->SetSize( FixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( FixedImage->GetOrigin() );
resample->SetOutputSpacing( FixedImage->GetSpacing() );
resample->SetOutputDirection( FixedImage->GetDirection() );
resample->SetDefaultPixelValue( 0 );
resample->Update();
ImageWriter->SetInput(resample->GetOutput());
try
{
ImageWriter->Update();
}
catch( itk::ExceptionObject & excp )
{
std::cerr << "Exception thrown " << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
FixedImageType::RegionType FixedRegion = FixedImage->GetBufferedRegion();
DeformField->SetRegions( FixedRegion );
DeformField->SetOrigin( FixedImage->GetOrigin() );
DeformField->SetSpacing( FixedImage->GetSpacing() );
FixedImageType:: DirectionType Fixeddirection =FixedImage->GetDirection();
DeformField->SetDirection(Fixeddirection);
DeformField->Allocate();
DeformField->Update();
//WRITE THE DEFORMATION FIELD
typedef itk::ImageRegionIterator< DisplacementFieldType > FieldIterator;
FieldIterator fi( DeformField, FixedRegion );
fi.GoToBegin();
DeformableTransformType::InputPointType fixedPointa;
DeformableTransformType::OutputPointType movingPointa;
DisplacementFieldType::IndexType indexa;
VectorPixelType displacement;
while( ! fi.IsAtEnd() )
{
indexa = fi.GetIndex();
DeformField->TransformIndexToPhysicalPoint( indexa, fixedPointa );
movingPointa = BtransformHigh->TransformPoint( fixedPointa );
displacement = movingPointa - fixedPointa;
fi.Set( displacement );
++fi;
}
FieldWriter->SetInput( DeformField );
//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;
}
return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20131211/d0d69676/attachment.html>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users
More information about the Community
mailing list