<div>Dear All,</div>
<div> </div>
<div>I tried to apply different optimizers for doing registration of the sample images in the BrainWeb.</div>
<div>I tried using both Powell and RegularStepGradientDescent Optimizer.</div>
<div>Both give me different results. The algorithm was modified from the sample programs; attached below.</div>
<div>Powell's result is as expected, though if I changed some parameters (i.e. step size, optimizer scales, multiresolution level ) I also got different results. Even, if I tried to give initial transformation close to the answer, it ended up to give wrong result.
</div>
<div> </div>
<div>My questions:</div>
<div>1. What's wrong with the code using the RegularStepGradientDescent Optimizer below ? Why couldn't it give the correct result ?</div>
<div>2. What's the significance of the parameters in the codes ?</div>
<div>3. Can I say that registration is dependant on the parameters that we specify ? So, it's like a guessing problem ? When would only can get result if we guess it correctly ?</div>
<div> </div>
<div>Could anybody answer me please.</div>
<div>Thank you in advance.</div>
<div>Regards.</div>
<div> </div>
<div>Dilla</div>
<div> </div>
<div>Codes:</div>
<div> </div>
<div>
<p>#if defined(_MSC_VER)<br>#pragma warning ( disable : 4786 )<br>#endif</p>
<p>#include "itkMultiResolutionImageRegistrationMethod.h"<br>#include "itkMultiResolutionPyramidImageFilter.h"<br>#include "itkMattesMutualInformationImageToImageMetric.h"<br>//#include "
itkLinearInterpolateImageFunction.h"<br>#include "itkBSplineInterpolateImageFunction.h"<br>//#include "itkPowellOptimizer.h"<br>#include "itkRegularStepGradientDescentOptimizer.h"<br>#include "
itkNormalVariateGenerator.h"</p>
<p>#include "itkImage.h"</p>
<p>#include "itkTimeProbesCollectorBase.h"</p>
<p>#include "itkEuler3DTransform.h"<br>#include "itkCenteredTransformInitializer.h"</p>
<p>#include "itkImageFileReader.h"<br>#include "itkImageFileWriter.h"</p>
<p>#include "itkResampleImageFilter.h"<br>#include "itkSubtractImageFilter.h"<br>#include "itkRescaleIntensityImageFilter.h"</p>
<p>#include "itkCommand.h"</p>
<p>class CommandIterationUpdate : public itk::Command <br>{<br>public:<br> typedef CommandIterationUpdate Self;<br> typedef itk::Command Superclass;<br> typedef itk::SmartPointer<Self> Pointer;<br>
itkNewMacro( Self );<br>protected:<br> CommandIterationUpdate() {};<br>public:<br> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br> typedef const OptimizerType * OptimizerPointer;</p>
<p> void Execute(itk::Object *caller, const itk::EventObject & event)<br> {<br> Execute( (const itk::Object *)caller, event);<br> }</p>
<p> void Execute(const itk::Object * object, const itk::EventObject & event)<br> {<br> OptimizerPointer optimizer = <br> dynamic_cast< OptimizerPointer >( object );<br> if( ! itk::IterationEvent().CheckEvent( &event ) )
<br> {<br> return;<br> }<br> std::cout << optimizer->GetCurrentIteration() << " ";<br> std::cout << optimizer->GetValue() << " ";<br> std::cout << optimizer->GetCurrentPosition() << std::endl;
<br> }<br>};</p>
<p>template <typename TRegistration><br>class RegistrationInterfaceCommand : public itk::Command <br>{<br>public:<br> typedef RegistrationInterfaceCommand Self;<br> typedef itk::Command Superclass;
<br> typedef itk::SmartPointer<Self> Pointer;<br> itkNewMacro( Self );<br>protected:<br> RegistrationInterfaceCommand() {};<br>public:<br> typedef TRegistration RegistrationType;
<br> typedef RegistrationType * RegistrationPointer;<br> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<br> typedef OptimizerType * OptimizerPointer;
<br> void Execute(itk::Object * object, const itk::EventObject & event)<br> {<br> if(itk::IterationEvent().CheckEvent( &event ))<br> {</p>
<p> RegistrationPointer registration =<br> dynamic_cast<RegistrationPointer>( object );<br> OptimizerPointer optimizer = dynamic_cast< OptimizerPointer >( <br> registration->GetOptimizer() );
</p>
<p> if ( registration->GetCurrentLevel() == 0 )<br> {<br> optimizer->SetMaximumStepLength( 0.1 ); <br> optimizer->SetMinimumStepLength( 0.001 );<br> }<br> else<br> {<br> double currentLevel = (double)registration->GetCurrentLevel();
<br> double totalLevel = (double)registration->GetNumberOfLevels();<br> optimizer->SetMaximumStepLength( 0.1 * currentLevel/totalLevel);<br> optimizer->SetMinimumStepLength( 0.001 * currentLevel/totalLevel);
<br> }<br> std::cout << "Interation: " << registration->GetCurrentLevel() << std::endl;<br> }<br> }<br> void Execute(const itk::Object * , const itk::EventObject & )<br> { return; }
<br>};</p>
<p>int main( )<br>{<br>/* if( argc < 4 )<br> {<br> std::cerr << "Missing Parameters " << std::endl;<br> std::cerr << "Usage: " << argv[0];<br> std::cerr << " fixedImageFile movingImageFile ";
<br> std::cerr << " outputImagefile [differenceAfterRegistration] ";<br> std::cerr << " [differenceBeforeRegistration] ";<br> std::cerr << " [initialStepLength] "<< std::endl;
<br> return EXIT_FAILURE;<br> }<br> */ <br> const unsigned int Dimension = 3;<br> typedef unsigned char PixelType;</p>
<p> typedef itk::Image< PixelType, Dimension > FixedImageType;<br> typedef itk::Image< PixelType, Dimension > MovingImageType;<br> typedef float InternalPixelType;<br> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
</p>
<p> typedef itk::Euler3DTransform< double > TransformType;</p>
<p> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;</p>
<p> typedef itk::MattesMutualInformationImageToImageMetric< <br> InternalImageType, <br> InternalImageType > MetricType;</p>
<p> typedef itk:: LinearInterpolateImageFunction< <br> InternalImageType,<br> double > InterpolatorType;</p>
<p><br>/* typedef itk::BSplineInterpolateImageFunction< <br> InternalImageType, double > InterpolatorType;<br>*/ </p>
<p> typedef itk::MultiResolutionImageRegistrationMethod< <br> InternalImageType, <br> InternalImageType > RegistrationType; </p>
<p> typedef itk::MultiResolutionPyramidImageFilter<<br> InternalImageType,<br> InternalImageType > FixedImagePyramidType;<br> typedef itk::MultiResolutionPyramidImageFilter<
<br> InternalImageType,<br> InternalImageType > MovingImagePyramidType;</p>
<p> MetricType::Pointer metric = MetricType::New();<br> OptimizerType::Pointer optimizer = OptimizerType::New();<br> InterpolatorType::Pointer interpolator = InterpolatorType::New();<br> RegistrationType::Pointer registration = RegistrationType::New();
<br> <br> registration->SetMetric( metric );<br> registration->SetOptimizer( optimizer );<br> registration->SetInterpolator( interpolator );</p>
<p><br> TransformType::Pointer transform = TransformType::New();<br> registration->SetTransform( transform );</p>
<p> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;<br> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;</p>
<p> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();<br> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();</p>
<p> fixedImageReader->SetFileName("brainweb1e1a10f20.mha");<br> movingImageReader->SetFileName( "brainweb1e1a10f20Rot10Tx15.mha" );</p>
<p> fixedImageReader->Update();<br> FixedImageType::SizeType size;<br> size = fixedImageReader->GetOutput()->GetLargestPossibleRegion().GetSize();</p>
<p> unsigned int numberOfBins = 30;<br> double percentOfSamples = 0.2;<br> unsigned int numberOfSamples = (unsigned int)(percentOfSamples*size[0]*size[1]*size[2]);<br> metric->SetNumberOfSpatialSamples(numberOfSamples);
<br> metric->SetNumberOfHistogramBins(numberOfBins);</p>
<p> typedef itk::CastImageFilter< <br> FixedImageType, InternalImageType > FixedCastFilterType;<br> typedef itk::CastImageFilter< <br> MovingImageType, InternalImageType > MovingCastFilterType;
</p>
<p> FixedCastFilterType::Pointer fixedCaster = FixedCastFilterType::New();<br> MovingCastFilterType::Pointer movingCaster = MovingCastFilterType::New();</p>
<p> fixedCaster->SetInput(fixedImageReader->GetOutput());<br> movingCaster->SetInput(movingImageReader->GetOutput());<br> <br> FixedImagePyramidType::Pointer fixedImagePyramid = <br> FixedImagePyramidType::New();
<br> MovingImagePyramidType::Pointer movingImagePyramid =<br> MovingImagePyramidType::New();<br> registration->SetFixedImagePyramid(fixedImagePyramid);<br> registration->SetMovingImagePyramid(movingImagePyramid);
</p>
<p> registration->SetFixedImage( fixedCaster->GetOutput() );<br> registration->SetMovingImage( movingCaster->GetOutput() );</p>
<p>////////////////////////////////////////////</p>
<p> InternalImageType::IndexType regStart;<br> regStart[0] = 0;<br> regStart[1] = 0;<br> regStart[2] = 0;</p>
<p> InternalImageType::SizeType regSize;<br> regSize[0] = size[0];<br> regSize[1] = size[1];<br> regSize[2] = size[2];</p>
<p> InternalImageType::RegionType region;<br> region.SetIndex(regStart);<br> region.SetSize(regSize);</p>
<p> registration->SetFixedImageRegion(region);</p>
<p>/////////////////////////////////////////////</p>
<p> fixedImageReader->Update();</p>
<p> movingImageReader->Update();</p>
<p><br> typedef FixedImageType::SpacingType SpacingType;<br> typedef FixedImageType::PointType OriginType;<br> typedef FixedImageType::RegionType RegionType;<br> typedef FixedImageType::SizeType SizeType;
</p>
<p> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();</p>
<p> const SpacingType fixedSpacing = fixedImage->GetSpacing();<br> OriginType fixedOrigin = fixedImage->GetOrigin();<br> const RegionType fixedRegion = fixedImage->GetLargestPossibleRegion(); <br> const SizeType fixedSize =
fixedRegion.GetSize();<br> </p>
<p> TransformType::InputPointType centerFixed;<br> <br> centerFixed[0] = regStart[0] + fixedSpacing[0] * regSize[0] / 2.0 ;<br> centerFixed[1] = regStart[1] + fixedSpacing[1] * regSize[1] / 2.0 ;<br> centerFixed[2] = regStart[2] + fixedSpacing[2] * regSize[2] /
2.0 ;<br> // moving image<br> MovingImageType::Pointer movingImage = movingImageReader->GetOutput();</p>
<p> const SpacingType movingSpacing = movingImage->GetSpacing();<br> const OriginType movingOrigin = movingImage->GetOrigin();<br> const RegionType movingRegion = movingImage->GetLargestPossibleRegion();<br>
const SizeType movingSize = movingRegion.GetSize();<br> <br> TransformType::InputPointType centerMoving;<br> <br> centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] / 2.0;<br> centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] /
2.0;<br> centerMoving[2] = movingOrigin[2] + movingSpacing[2] * movingSize[2] / 2.0;</p>
<p> transform->SetCenter( centerFixed );<br> transform->SetTranslation( centerMoving - centerFixed );</p>
<p> // transform->SetIdentity( );<br> <br> registration->SetInitialTransformParameters( transform->GetParameters() );</p>
<p> typedef OptimizerType::ScalesType OptimizerScalesType;<br> OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );<br> const double translationScale = 1.0;</p>
<p> optimizerScales[0] = 1.0;//rotation<br> optimizerScales[1] = 1.0;//rotation<br> optimizerScales[2] = 1.0;//rotation<br> optimizerScales[3] = translationScale;<br> optimizerScales[4] = translationScale;<br> optimizerScales[5] = translationScale;
</p>
<p> optimizer->SetScales( optimizerScales );</p>
<p>/* double initialStepLength = 5;</p>
<p><br> // Power's Optimizer<br> optimizer->SetMaximumIteration(200);<br> optimizer->SetStepLength(2.0);<br> optimizer->SetStepTolerance(0.01);<br>// optimizer->SetMaximize(false);<br>*/<br> // Create the Command observer and register it with the optimizer.
<br> //<br> optimizer->SetNumberOfIterations( 200 );<br>//optimizer->SetMaximumStepLength( 10.00 ); <br>// optimizer->SetMinimumStepLength( 2 );<br> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
<br> optimizer->AddObserver( itk::IterationEvent(), observer );<br> // Create the command observer for registration<br> //<br> typedef RegistrationInterfaceCommand<RegistrationType> CommandType;<br> CommandType::Pointer command = CommandType::New();
<br> registration->AddObserver( itk::IterationEvent(), command );</p>
<p> // if ( argc > 6 ) {<br>// registration->SetNumberOfLevels(std::atoi(argv[6]));<br> // } else {<br> registration->SetNumberOfLevels( 5 );<br> // }<br> <br> itk::TimeProbesCollectorBase timer;<br> try <br>
{ <br> timer.Start("registration");<br> registration->StartRegistration(); <br> timer.Stop("registration");<br> } <br> catch( itk::ExceptionObject & err ) <br> { <br> std::cerr << "ExceptionObject caught !" << std::endl;
<br> std::cerr << err << std::endl; <br> return EXIT_FAILURE;<br> } <br> <br> timer.Report();</p>
<p> OptimizerType::ParametersType finalParameters = <br> registration->GetLastTransformParameters();</p>
<p> const double finalAngleX = finalParameters[0];<br> const double finalAngleY = finalParameters[1];<br> const double finalAngleZ = finalParameters[2];<br> const double finalTranslationX = finalParameters[3];
<br> const double finalTranslationY = finalParameters[4];<br> const double finalTranslationZ = finalParameters[5];</p>
<p> const unsigned int numberOfIterations = optimizer->GetCurrentIteration();</p>
<p> const double bestValue = optimizer->GetValue();</p>
<p><br> // Print out results<br> //<br> const double finalAngleXInDegrees = finalAngleX * 45.0 / atan(1.0);<br> const double finalAngleYInDegrees = finalAngleY * 45.0 / atan(1.0);<br> const double finalAngleZInDegrees = finalAngleZ *
45.0 / atan(1.0);</p>
<p> std::cout << "Result = " << std::endl;</p>
<p> std::cout << " Angle X (degrees) = " << finalAngleXInDegrees << std::endl;<br> std::cout << " Angle Y (degrees) = " << finalAngleYInDegrees << std::endl;
<br> std::cout << " Angle Z (degrees) = " << finalAngleZInDegrees << std::endl;<br> std::cout << " Translation X = " << finalTranslationX << std::endl;<br> std::cout << " Translation Y = " << finalTranslationY << std::endl;
<br> std::cout << " Translation Z = " << finalTranslationZ << std::endl;<br> std::cout << " Iterations = " << numberOfIterations << std::endl;<br> std::cout << " Metric value = " << bestValue << std::endl;
</p>
<p> return EXIT_SUCCESS;<br>}</p></div>