[ITK Community] [Insight-users] Bspline artefact multi resolution - origin offset - Bug?
Luis Ibanez
luis.ibanez at kitware.com
Fri Dec 13 09:04:57 EST 2013
Emma,
It is likely that your settings for the Grid of the
BSplineDeformableTransform are not fitting well to the FixedImage.
You may want to use the helper class:
http://www.itk.org/Doxygen/html/classitk_1_1BSplineTransformInitializer.html
that will take care of doing this properly, so you don't have to adjust
those parameters by hand.
More details are available in the article in the Insight Journal
http://www.insight-journal.org/browse/publication/216
Thanks
Luis
On Wed, Dec 11, 2013 at 5:26 PM, Emma Saunders <emmasaunders123 at gmail.com>wrote:
> 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;
> }
>
>
>
>
> _____________________________________
> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20131213/602522ab/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