[Insight-users] Help needed for 2D multi-modality Registration.
Luis Ibanez
luis.ibanez at kitware.com
Wed Mar 25 20:12:07 EDT 2009
Hi Reeju,
The problem seem to be that you missed to extend the array that
carry information related to the transform.
In particular:
ParametersType initialParameters(
transform->GetNumberOfParameters() );
initialParameters[0] = 0.0; // Initial offset in mm along X
initialParameters[1] = 0.0; // Initial offset in mm along Y
initialParameters[2] = 0.0; // vertical scaling
registration->SetInitialTransformParameters( initialParameters );
was for the three components of the translation transform,
but now that you are using an Affine transform in 3D, you
need 12 parameters. The first nine are the components of the
rotation+scaling+shearing matrix, while the last three are
the translation components. Something like:
ParametersType initialParameters(
transform->GetNumberOfParameters() );
initialParameters[0] = 1.0;
initialParameters[1] = 0.0;
initialParameters[2] = 0.0;
initialParameters[3] = 0.0;
initialParameters[4] = 1.0;
initialParameters[5] = 0.0;
initialParameters[6] = 0.0;
initialParameters[7] = 0.0;
initialParameters[8] = 1.0;
initialParameters[9] = 0.0;
initialParameters[10] = 0.0;
initialParameters[11] = 0.0;
registration->SetInitialTransformParameters( initialParameters );
You will also need to set up the OptimizerScales, in order
to balance the dynamic range of the matrix components (usually
in the range close to -1.0 : 1.0) with the one of the translation
componets (usually in the range -100:100).
This will be code like:
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales(
transform->GetNumberOfParameters() );
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
... // all with 1.0
optimizerScales[8] = 1.0;
optimizerScales[ 9] = translationScale;
optimizerScales[10] = translationScale;
optimizerScales[11] = translationScale;
optimizer->SetScales( optimizerScales );
This should be inserted *before* you call StartRegistration();
For an explanation of this scales, please look at the
ITK Software Guide
http://www.itk.org/ItkSoftwareGuide.pdf
in Section 8.6.1, in pdf-page 383.
Finally, the interpretation of the output parameters must also
be updated. You still have the code from the 3 parameters of
the translation transform:
ParametersType finalParameters =
registration->GetLastTransformParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];
double TranslationAlongZ = finalParameters[2];
while you should have now something like:
ParametersType finalParameters =
registration->GetLastTransformParameters();
double M00 = finalParameters[0];
double M01 = finalParameters[1];
double M02 = finalParameters[2];
...
double M22 = finalParameters[8];
double TranslationAlongX = finalParameters[ 9];
double TranslationAlongY = finalParameters[10];
double TranslationAlongZ = finalParameters[11];
For an explanation of the AffineTransform Parameters,
please look at the documentation in the ITK Doxygen page
for this Transform:
http://www.itk.org/Insight/Doxygen/html/classitk_1_1AffineTransform.html
Regards,
Luis
---------------------
Reeju Pokharel wrote:
> Hello Luis,
>
> I was wondering if you could help me with a 2D multi-modality
> registration. I have four images obtained from SEM and I need to
> register the three moving images with one fixed image. I have three
> transform parameters. The images are taken from different angles,
> therefore they behave as if they were taken from different imaging
> modalities. Since I'm new to itk, I was trying to learn the basics of
> registration by running the example on image registration that uses
> mutual information metric and gradientdescent optimizer. However, the
> example uses tranlation transform and since I have shearing in my image
> I need to use affine transform. Therefore, I changed the code and tried
> to use affine transform instead of translation transform. Now the code
> doesnt work and I do not know what I'm doing wrong here. Could you
> please give me some comments on this. Here is the code:
>
> #include "itkImageRegistrationMethod.h"
> #include "itkTranslationTransform.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"
> #include "itkCheckerBoardImageFilter.h"
> #include "itkTIFFImageIO.h"
> #include "itkCenteredTransformInitializer.h"
> #include "itkAffineTransform.h"
>
>
> #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( ! itk::IterationEvent().CheckEvent( &event ) )
> {
> return;
> }
> std::cout << optimizer->GetCurrentIteration() << " ";
> std::cout << optimizer->GetValue() << " ";
> std::cout << optimizer->GetCurrentPosition() << std::endl;
>
> // Print the angle for the trace plot
> vnl_matrix<double> p(2, 2);
> p[0][0] = (double) optimizer->GetCurrentPosition()[0];
> p[0][1] = (double) optimizer->GetCurrentPosition()[1];
> p[1][0] = (double) optimizer->GetCurrentPosition()[2];
> p[1][1] = (double) optimizer->GetCurrentPosition()[3];
> vnl_svd<double> svd(p);
> vnl_matrix<double> r(2, 2);
> r = svd.U() * vnl_transpose(svd.V());
> double angle = asin(r[1][0]);
> std::cout << " AffineAngle: " << angle * 45.0 / atan(1.0) <<
> 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::cerr << "[checkerBoardBefore] [checkerBoardAfter]" << std::endl;
> return EXIT_FAILURE;
> }
>
>
> // Software Guide : BeginCodeSnippet
> const unsigned int Dimension = 2;
> typedef unsigned short PixelType;
>
> typedef itk::Image< PixelType, Dimension > FixedImageType;
> typedef itk::Image< PixelType, Dimension > MovingImageType;
> // Software Guide : EndCodeSnippet
>
>
> // Software Guide : BeginCodeSnippet
> typedef float InternalPixelType;
> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
> // Software Guide : EndCodeSnippet
>
> // Software Guide : BeginCodeSnippet
> typedef itk::AffineTransform<
> double,
> Dimension > TransformType;
> // Software Guide : EndCodeSnippet
>
>
> //typedef itk::TranslationTransform< double, Dimension > 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 );
>
>
>
> // Software Guide : BeginCodeSnippet
> MetricType::Pointer metric = MetricType::New();
> registration->SetMetric( metric );
> // Software Guide : EndCodeSnippet
>
>
>
>
> // Software Guide : BeginCodeSnippet
> metric->SetFixedImageStandardDeviation( 0.4 );
> metric->SetMovingImageStandardDeviation( 0.4 );
> // Software Guide : EndCodeSnippet
>
>
>
> 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] );
>
>
>
> // Software Guide : BeginCodeSnippet
> typedef itk::CenteredTransformInitializer<
> TransformType,
> FixedImageType,
> MovingImageType >
> TransformInitializerType;
> TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
> initializer->SetTransform( transform );
> initializer->SetFixedImage( fixedImageReader->GetOutput() );
> initializer->SetMovingImage( movingImageReader->GetOutput() );
> initializer->MomentsOn();
> initializer->InitializeTransform();
> // Software Guide : EndCodeSnippet
>
>
>
> // Software Guide : BeginCodeSnippet
> typedef itk::NormalizeImageFilter<
> FixedImageType,
> InternalImageType
> > FixedNormalizeFilterType;
>
> typedef itk::NormalizeImageFilter<
> MovingImageType,
> InternalImageType
> > MovingNormalizeFilterType;
>
> FixedNormalizeFilterType::Pointer fixedNormalizer =
> FixedNormalizeFilterType::New();
>
> MovingNormalizeFilterType::Pointer movingNormalizer =
>
> MovingNormalizeFilterType::New();
>
>
> // Software Guide : BeginCodeSnippet
> typedef itk::DiscreteGaussianImageFilter<
> InternalImageType,
> InternalImageType
> > GaussianFilterType;
>
> GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
> GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
>
> fixedSmoother->SetVariance( 2.0 );
> movingSmoother->SetVariance( 2.0 );
> // Software Guide : EndCodeSnippet
>
>
>
> // Software Guide : BeginCodeSnippet
> fixedNormalizer->SetInput( fixedImageReader->GetOutput() );
> movingNormalizer->SetInput( movingImageReader->GetOutput() );
>
> fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
> movingSmoother->SetInput( movingNormalizer->GetOutput() );
>
> registration->SetFixedImage( fixedSmoother->GetOutput() );
> registration->SetMovingImage( movingSmoother->GetOutput() );
> // Software Guide : EndCodeSnippet
>
>
> fixedNormalizer->Update();
> FixedImageType::RegionType fixedImageRegion =
> fixedNormalizer->GetOutput()->GetBufferedRegion();
> registration->SetFixedImageRegion( fixedImageRegion );
>
> typedef RegistrationType::ParametersType ParametersType;
> ParametersType initialParameters( transform->GetNumberOfParameters() );
>
> initialParameters[0] = 0.0; // Initial offset in mm along X
> initialParameters[1] = 0.0; // Initial offset in mm along Y
> initialParameters[2] = 0.0; // vertical scaling
>
> registration->SetInitialTransformParameters( initialParameters );
>
>
>
> // Software Guide : BeginCodeSnippet
> const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();
>
> const unsigned int numberOfSamples =
> static_cast< unsigned int >( numberOfPixels *
> 0.01 );
>
> metric->SetNumberOfSpatialSamples( numberOfSamples );
> // Software Guide : EndCodeSnippet
>
>
>
> // Software Guide : BeginCodeSnippet
> optimizer->SetLearningRate( 15.0 );
> optimizer->SetNumberOfIterations( 1 );
> optimizer->MaximizeOn();
> // Software Guide : EndCodeSnippet
>
>
>
>
> // 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 EXIT_FAILURE;
> }
>
> ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
> double TranslationAlongX = finalParameters[0];
> double TranslationAlongY = finalParameters[1];
> double VerticalScaling = finalParameters[2];
>
>
> unsigned int numberOfIterations = optimizer->GetCurrentIteration();
>
> double bestValue = optimizer->GetValue();
>
>
> // Print out results
> //
> std::cout << std::endl;
> std::cout << "Result = " << std::endl;
> std::cout << " Translation X = " << TranslationAlongX << std::endl;
> std::cout << " Translation Y = " << TranslationAlongY << std::endl;
> std::cout << " Vertical Scaling = " << VerticalScaling << std::endl;
> std::cout << " Iterations = " << numberOfIterations << std::endl;
> std::cout << " Metric value = " << bestValue << std::endl;
> std::cout << " Numb. Samples = " << numberOfSamples << 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() );
>
> FixedImageType::Pointer fixedImage = fixedImageReader->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();
>
>
> itk::TIFFImageIO::Pointer tiffio = itk::TIFFImageIO::New();
> tiffio->SetUseCompression(0);
> writer->SetImageIO( tiffio);
> writer->SetFileName( argv[3] );
>
>
> caster->SetInput( resample->GetOutput() );
> writer->SetInput( caster->GetOutput() );
> writer->Update();
>
>
> // Generate checkerboards before and after registration
> //
> typedef itk::CheckerBoardImageFilter< FixedImageType >
> CheckerBoardFilterType;
>
> CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
>
> checker->SetInput1( fixedImage );
> checker->SetInput2( resample->GetOutput() );
>
> caster->SetInput( checker->GetOutput() );
> writer->SetInput( caster->GetOutput() );
>
> // Before registration
> TransformType::Pointer identityTransform = TransformType::New();
> identityTransform->SetIdentity();
> resample->SetTransform( identityTransform );
>
> if( argc > 4 )
> {
> writer->SetFileName( argv[4] );
> writer->Update();
> }
>
>
> // After registration
> resample->SetTransform( finalTransform );
> if( argc > 5 )
> {
> writer->SetFileName( argv[5] );
> writer->Update();
> }
>
>
>
>
> return EXIT_SUCCESS;
> }
>
>
>
> *Error Message*:
>
>
> ExceptionObject caught !
>
> itk::ExceptionObject (013DF87C)
> Location: "void __thiscall
> itk::MutualInformationImageToImageMetric<class itk::I
> mage<float,2>,class itk::Image<float,2> >::SampleFixedImageDomain(class
> std::vec
> tor<class itk::MutualInformationImageToImageMetric<class
> itk::Image<float,2>,cla
> ss itk::Image<float,2> >::SpatialSample,class std::allocator<class
> itk::MutualIn
> formationImageToImageMetric<class itk::Image<float,2>,class
> itk::Image<float,2>
> >::SpatialSample> > &) const"
> File:
> c:\miia\insighttoolkit-3.4.0\code\algorithms\itkMutualInformationImageToIm
> ageMetric.txx
> Line: 193
> Description: itk::ERROR: MutualInformationImageToImageMetric(018B46A0):
> All the
> sampled point mapped to outside of the moving image
>
>
>
> Thank you for your help. I really appreciate it.
>
> Reeju Pokharel.
>
>
More information about the Insight-users
mailing list