[Insight-users] Mutual Information rotation & translation and
the metric value (with the code)
Karthik Krishnan
Karthik.Krishnan at kitware.com
Tue Jul 26 16:31:17 EDT 2005
Reduce your learning rate. The idea is to watch the metric trace and
reduce the learning rate if it oscillates too much. You may also want to
use the RegularStepGradientDescent instead of the GradientDescent. This
halves the step size at every iteration and attempts to provide stability.
Also, there is little sense in optimizing both the translation and the
center. The Centered2DTransform has 5 parameters, [angle , Cx, Cy, Tx,
Ty ]. Really, you could represent rotation about any center Cx, Cy with
an equivalent transform that is rotated about the origin [0, 0], the
value of the translation calculated accordingly. Not that there is
anything wrong in optimizing both the center and the translation, but it
just increases the work of the optimizer.
For the given image using scales of 10 for rotation, and the your chosen
scales for translation and 100000 (a large number for the center), and a
learning rate of 0.02, and 500 spatial samples for the metric, I get a
pretty good registration with your code.
HTH
karthik
Seniha Esen Yuksel wrote:
> Hi Lydia and other listeners,
>
> I am still trying to get a mutual information registration code that
> both translates and rotates the images using centered2DTransform. So I
> added the optimizerScales to my previous code and played with the
> translationScale parameters, but no matter I did I am still getting an
> exception after going to big translation values (such as 391 in y
> direction) or minus metric values. Initially (ImageRegistration5.cxx),
> for these two images ( BrainProtonDensitySliceBorder20.png and
> BrainProtonDensitySliceR10X13Y17.png from ITK Data), the
> translationScale is given as 1000, which didn't work. So I read
> _http://www.itk.org/pipermail/insight-users/2004-July/009558.html_,
> and calculated the translationScale to be 1/ ( (221^2 + 257^2) ) *10
> = 3390, which didn't work either. I also tried couple of other values
> but I still ge! t the exception after some iterations. Do you think
> there is another problem in my code that I need to change?
>
> I am attaching the code and your help will be very much appreciated,
>
> Thanks
> Seniha
>
>
> ------------------------------------------------------------------------
> Start your day with Yahoo! - make it your home page
> <http://us.rd.yahoo.com/evt=34442/*http://www.yahoo.com/r/hs>
>
>------------------------------------------------------------------------
>
>/*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: ImageRegistration2.cxx,v $
> Language: C++
> Date: $Date: 2004/12/28 14:42:49 $
> Version: $Revision: 1.33 $
>
> Copyright (c) Insight Software Consortium. All rights reserved.
> See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notices for more information.
>
>=========================================================================*/
>
>
>#include "itkImageRegistrationMethod.h"
>#include "itkCenteredRigid2DTransform.h"
>#include "itkMutualInformationImageToImageMetric.h"
>#include "itkLinearInterpolateImageFunction.h"
>#include "itkGradientDescentOptimizer.h"
>#include "itkImage.h"
>#include "itkNormalizeImageFilter.h"
>#include "itkDiscreteGaussianImageFilter.h"
>
>#include "itkImageFileReader.h"
>#include "itkImageFileWriter.h"
>
>#include "itkResampleImageFilter.h"
>#include "itkCastImageFilter.h"
>
>
>// The following section of code implements a Command observer
>// that will monitor the evolution of the registration process.
>//
>#include "itkCommand.h"
>class CommandIterationUpdate : public itk::Command
>{
>public:
> typedef CommandIterationUpdate Self;
> typedef itk::Command Superclass;
> typedef itk::SmartPointer<Self> Pointer;
> itkNewMacro( Self );
>protected:
> CommandIterationUpdate() {};
>public:
> typedef itk::GradientDescentOptimizer 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( typeid( event ) != typeid( itk::IterationEvent ) )
> {
> return;
> }
> std::cout << optimizer->GetCurrentIteration() << " ";
> std::cout << optimizer->GetValue() << " ";
> std::cout << optimizer->GetCurrentPosition() << 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::cerr << " outputImagefile "<< std::endl;
> return 1;
> }
>*/
>
> const unsigned int Dimension = 2;
> typedef unsigned short PixelType;
>
> typedef itk::Image< PixelType, Dimension > FixedImageType;
> typedef itk::Image< PixelType, Dimension > MovingImageType;
>
> typedef float InternalPixelType;
> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
>
> typedef itk::CenteredRigid2DTransform< double > TransformType;
>
> typedef itk::GradientDescentOptimizer OptimizerType;
> typedef itk::LinearInterpolateImageFunction<
> InternalImageType,
> double > InterpolatorType;
> typedef itk::ImageRegistrationMethod<
> InternalImageType,
> InternalImageType > RegistrationType;
>
> typedef itk::MutualInformationImageToImageMetric<
> InternalImageType,
> InternalImageType > MetricType;
>
>
> TransformType::Pointer transform = TransformType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
>
> registration->SetOptimizer( optimizer );
> registration->SetTransform( transform );
> registration->SetInterpolator( interpolator );
>
>
> MetricType::Pointer metric = MetricType::New();
> registration->SetMetric( metric );
>
> metric->SetFixedImageStandardDeviation( 0.4 );
> metric->SetMovingImageStandardDeviation( 0.4 );
> metric->SetNumberOfSpatialSamples( 50 );
>
>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>
> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
>
> //fixedImageReader->SetFileName( argv[1] );
>//movingImageReader->SetFileName( argv[2] );
> //fixedImageReader->SetFileName( "original.png" );
> //movingImageReader->SetFileName( "rotate.png");
>
> //fixedImageReader->SetFileName( "kidney_original_noise_phantom.bmp" );
>// movingImageReader->SetFileName( "kidney_rotated_noise_phantom.bmp");
>
> fixedImageReader->SetFileName( "BrainProtonDensitySliceBorder20.png" );
> movingImageReader->SetFileName( "BrainProtonDensitySliceR10X13Y17.png" );
>
>
>
>
> typedef itk::NormalizeImageFilter<
> FixedImageType,
> InternalImageType
> > FixedNormalizeFilterType;
>
> typedef itk::NormalizeImageFilter<
> MovingImageType,
> InternalImageType
> > MovingNormalizeFilterType;
>
> FixedNormalizeFilterType::Pointer fixedNormalizer =
> FixedNormalizeFilterType::New();
>
> MovingNormalizeFilterType::Pointer movingNormalizer =
> MovingNormalizeFilterType::New();
>
> typedef itk::DiscreteGaussianImageFilter<
> InternalImageType,
> InternalImageType
> > GaussianFilterType;
>
> GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
> GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
>
> fixedSmoother->SetVariance( 2.0 );
> movingSmoother->SetVariance( 2.0 );
>
>
> fixedNormalizer->SetInput( fixedImageReader->GetOutput() );
> movingNormalizer->SetInput( movingImageReader->GetOutput() );
>
> fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
> movingSmoother->SetInput( movingNormalizer->GetOutput() );
>
>
> // fixedSmoother->Update();
> // movingSmoother->Update();
>
> registration->SetFixedImage( fixedSmoother->GetOutput() );
> registration->SetMovingImage( movingSmoother->GetOutput() );
>
>
> fixedNormalizer->Update();
> //movingNormalizer->Update();
>
> registration->SetFixedImageRegion(
> fixedNormalizer->GetOutput()->GetBufferedRegion() );
>
>
> fixedImageReader->Update();
> movingImageReader->Update();
>
>
> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
>
> const FixedImageType::SpacingType&
> fixedSpacing = fixedImage->GetSpacing();
> const FixedImageType::PointType&
> fixedOrigin = fixedImage->GetOrigin();
>
> FixedImageType::SizeType fixedSize =
> fixedImage->GetLargestPossibleRegion().GetSize();
>
> TransformType::InputPointType centerFixed;
>
> centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
> centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
>
>
>
>
> MovingImageType::Pointer movingImage = movingImageReader->GetOutput();
>
> const MovingImageType::SpacingType&
> movingSpacing = movingImage->GetSpacing();
> const MovingImageType::PointType&
> movingOrigin = movingImage->GetOrigin();
>
> MovingImageType::SizeType movingSize =
> movingImage->GetLargestPossibleRegion().GetSize();
>
> TransformType::InputPointType centerMoving;
>
> centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] / 2.0;
> centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] / 2.0;
>
>
> transform->SetCenter( centerFixed );
> transform->SetTranslation( centerMoving - centerFixed );
> transform->SetAngle( 0.0 );
>
> // Now we pass the current transform's parameters as the initial
> // parameters to be used when the registration process starts.
>
>
> registration->SetInitialTransformParameters( transform->GetParameters() );
>
>
> typedef OptimizerType::ScalesType OptimizerScalesType;
> OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
>// const double translationScale = 1.0 / 1000.0;
>const double translationScale = 1.0 /3390.0;
>
>// const double translationScale1 = 1.0 /110.0;
>// const double translationScale2 = 1.0 /128.0;
>
> optimizerScales[0] = 1.0;
> optimizerScales[1] = translationScale;
> optimizerScales[2] = translationScale;
> optimizerScales[3] = translationScale;
> optimizerScales[4] = translationScale;
>
> optimizer->SetScales( optimizerScales );
>
>
>
>
>
>
>
> optimizer->SetLearningRate( 10);
> optimizer->SetNumberOfIterations( 500 );
> optimizer->MaximizeOn();
>
>
> // Create the Command observer and register it with the optimizer.
> //
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
>
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cout << "ExceptionObject caught !" << std::endl;
> std::cout << err << std::endl;
> return -1;
> }
>
>
>
> OptimizerType::ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
> const double finalAngle = finalParameters[0];
> const double finalRotationCenterX = finalParameters[1];
> const double finalRotationCenterY = finalParameters[2];
> const double finalTranslationX = finalParameters[3];
> const double finalTranslationY = finalParameters[4];
>
>
>
> unsigned int numberOfIterations = optimizer->GetCurrentIteration();
>
> double bestValue = optimizer->GetValue();
>
>
> const double finalAngleInDegrees = finalAngle * 45.0 / atan(1.0);
>
> std::cout << "Result = " << std::endl;
> std::cout << " Angle (radians) = " << finalAngle << std::endl;
> std::cout << " Angle (degrees) = " << finalAngleInDegrees << std::endl;
> std::cout << " Center X = " << finalRotationCenterX << std::endl;
> std::cout << " Center Y = " << finalRotationCenterY << std::endl;
> std::cout << " Translation X = " << finalTranslationX << std::endl;
> std::cout << " Translation Y = " << finalTranslationY << std::endl;
> std::cout << " Iterations = " << numberOfIterations << std::endl;
> std::cout << " Metric value = " << bestValue << std::endl;
>
>
>
>
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
>
> TransformType::Pointer finalTransform = TransformType::New();
>
> finalTransform->SetParameters( finalParameters );
>
> ResampleFilterType::Pointer resample = ResampleFilterType::New();
>
> resample->SetTransform( finalTransform );
> resample->SetInput( movingImageReader->GetOutput() );
>
>
> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
> resample->SetOutputOrigin( fixedImage->GetOrigin() );
> resample->SetOutputSpacing( fixedImage->GetSpacing() );
> resample->SetDefaultPixelValue( 100 );
>
>
> typedef unsigned char OutputPixelType;
>
> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
>
> typedef itk::CastImageFilter<
> FixedImageType,
> OutputImageType > CastFilterType;
>
> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>
>
> WriterType::Pointer writer = WriterType::New();
> CastFilterType::Pointer caster = CastFilterType::New();
>
>
> writer->SetFileName( "result.png" );
> // writer->SetFileName( argv[3] );
>
>
> caster->SetInput( resample->GetOutput() );
> writer->SetInput( caster->GetOutput() );
> writer->Update();
>
>
>
> return 0;
>}
>
>
>
>------------------------------------------------------------------------
>
>_______________________________________________
>Insight-users mailing list
>Insight-users at itk.org
>http://www.itk.org/mailman/listinfo/insight-users
>
>
More information about the Insight-users
mailing list