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