[Insight-users] Mutual Information rotation & translation and the
metric value (with the code)
Lydia Ng
lydiang at gmail.com
Fri Jul 15 16:26:04 EDT 2005
Hi Seniha,
This is most likely a parameter scaling problem. That is, a unit
change in rotation has a dramatically bigger effect that a unit change
in translation.
In ImageRegistration5.cxx there is a section of code (line 330
onwards) that shows how to set the "scales" for the optimizer. You
could try adding something similar in your example.
- Lydia
On 7/14/05, Seniha Esen Yuksel <eseny99 at yahoo.com> wrote:
>
> Hi,
> I am trying to write a mutual information registration code (Viola Wells)
> that implements both centered rotation and translation. Therefore, I tried
> to combine
> the examples ImageRegistration2.cxx (mutual information using translation
> transform) and ImageRegistration5.cxx (centeredrigid2D transform using
> MeanSquaresImageToImageMetric). However, when I run my code to register the
> given images
> BrainProtonDensitySliceBorder20.png and
> BrainProtonDensitySliceR10X13Y17.png, the angle
> values are extremely high (-3007.19) and I couldn't find where the problem
> is.
>
> Another question is, if I understood correctly, in these examples, we are
> taking
> the final registration parameters as the best values. However, the last
> metric value
> is not necessarily the best metric value, ie. we are not looking for the
> best metric
> that maximizes the mutual information but only taking the last one when we
> do "finalTransform->SetParameters( finalParameters )". Please correct me if
> I misunderstood this part or I would appreciate if you could explain why it
> was not done as such or how I can do it.
>
>
> Thank you very much,
> Best regards,
> Esen
>
>
>
> Seniha Esen Yuksel
> Research Assistant
> Computer Vision and Image Processing Lab,
> Lutz Hall Rm 414,
> University of Louisville,
> Louisville, KY
> work ph: 502- 8521528
> home ph: 502- 8549856
> email: esen at cvip.uofl.edu // esenyuksel at ieee.org
> web: http://www.cvip.louisville.edu/~esen/
>
> __________________________________________________
> Do You Yahoo!?
> Tired of spam? Yahoo! Mail has the best spam protection around
> http://mail.yahoo.com
> /*=========================================================================
>
> 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");
>
>
>
> 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() );
>
>
>
>
> optimizer->SetLearningRate( 20);
> 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] );
>
> //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