[Insight-users] Affine registration labels
Karthik Krishnan
Karthik.Krishnan at kitware.com
Wed Dec 14 11:23:09 EST 2005
You need to call SetInitialTransformParamters( .. ) before you start the
registration process.
The size mismatch it is complaining about it between the
m_InitialTransformParameters ivar in the ImageRegistartionMethod and the
number of parameters the transform has.
Its also incorrect to invoke registration without setting the initial
parameters cause you need to give the optimizers a point in parametric
space to start from.
you could just have
transform->SetIdentity()
registrtration->SetInitialTransformParameters( transform->GetParameters() );
karthik
Roger Alvaredo wrote:
> Hello,
>
> I am a new user of ITK. My purpose is to register two label map images
> (in 3D) with Affine Transformation.
>
> I take an example from ITK and I changed the transformation and the
> dimension of the images. The code compiles but when I launch it I get
> this error:
>
> itk::ExceptionObject (0105F158)
> Location: "Unknown"
> File:
> d:\itk_source\insighttoolkit-2.0.1\code\algorithms\itkImageRegistrationMet
> hod.txx
> Line: 216
> Description: itk::ERROR: ImageRegistrationMethod(00278618): Size
> mismatch betwee
> n initial parameter and transform
>
> I attach the code.
>
> Can somebody help me? Let me know if there is something wrong in the
> code or other advices.
>
> Thanks
>
> ----------------------
> Roger Alvaredo
>
>------------------------------------------------------------------------
>
>/*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile:ImageRegistration9.cxx,v $
> Language: C++
> Date: $Date: 2005/08/31 14:14:15 $
> Version: $Revision: 1.13 $
>
> 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.
>
>=========================================================================*/
>#if defined(_MSC_VER)
>#pragma warning ( disable : 4786 )
>#endif
>
>
>
>
>// Software Guide : BeginLatex
>//
>// This example illustrates the use of the image registration framework in
>// Insight to align two label maps. Common structures are assumed to
>// use the same label. The registration metric simply counts the
>// number of corresponding pixels that have the same label.
>//
>//
>// Software Guide : EndLatex
>
>
>// Software Guide : BeginCodeSnippet
>#include "itkImageRegistrationMethod.h"
>#include "itkAffineTransform.h"
>#include "itkMatchCardinalityImageToImageMetric.h"
>#include "itkNearestNeighborInterpolateImageFunction.h"
>#include "itkAmoebaOptimizer.h"
>// Software Guide : EndCodeSnippet
>
>#include "itkImage.h"
>
>#include "itkImageFileReader.h"
>#include "itkImageFileWriter.h"
>
>#include "itkResampleImageFilter.h"
>#include "itkCastImageFilter.h"
>#include "itkFileOutputWindow.h"
>
>//
>// The following piece of code implements an 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::AmoebaOptimizer 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->GetCachedValue() << " ";
>// std::cout << optimizer->GetCachedCurrentPosition() << std::endl;
> }
>};
>
>
>int main( int argc, char *argv[] )
>{
> if( argc < 3 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0];
> std::cerr << " fixedImageFile movingImageFile ";
> std::cerr << " outputImagefile" << std::endl;
> return 1;
> }
>
> itk::FileOutputWindow::Pointer fow = itk::FileOutputWindow::New();
> fow->SetInstance( fow );
>
> // The types of each one of the components in the registration methods should
> // be instantiated. First, we select the image dimension and the type for
> // representing image pixels.
> //
> const unsigned int Dimension = 3;
> typedef float PixelType;
>
>
> // The types of the input images are instantiated by the following lines.
> //
> typedef itk::Image< PixelType, Dimension > FixedImageType;
> typedef itk::Image< PixelType, Dimension > MovingImageType;
>
> // Software Guide : BeginLatex
> // The transform that will map one image space into the other is defined
> // below.
> // Software Guide : EndLatex
> // Software Guide : BeginCodeSnippet
> typedef itk::AffineTransform< double, Dimension > TransformType;
> // Software Guide : EndCodeSnippet
>
> // Software Guide : BeginLatex
> // An optimizer is required to explore the parameter space of the transform
> // in search of optimal values of the metric. The metric selected
> // does not require analytical derivatives of its cost function.
> // Software Guide : EndLatex
> // Software Guide : BeginCodeSnippet
> typedef itk::AmoebaOptimizer OptimizerType;
> // Software Guide : EndCodeSnippet
>
> // Software Guide : BeginLatex
> // The metric will compare how well the two images match each
> // other. Metric types are usually parametrized by the image types
> // as can be seen in the following type declaration. The metric
> // selected here is suitable for comparing two label maps where the
> // labels are consistent between the two maps. This metric
> // measures the percentage of pixels that exactly match or
> // mismatch.
> // Software Guide : EndLatex
> // Software Guide : BeginCodeSnippet
> typedef itk::MatchCardinalityImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;
> // Software Guide : EndCodeSnippet
>
> // Finally, the type of the interpolator is declared. The
> // interpolator will evaluate the moving image at non-grid
> // positions.
> // Software Guide : BeginLatex
> // Since we are registering label maps, we use a
> // NearestNeighborInterpolateImageFunction to ensure subpixel
> // values are not interpolated (to labels that do not exist).
> // Software Guide : EndLatex
> // Software Guide : BeginCodeSnippet
> typedef itk:: NearestNeighborInterpolateImageFunction<
> MovingImageType,
> double > InterpolatorType;
> // Software Guide : EndCodeSnippet
>
> // The registration method type is instantiated using the types of the
> // fixed and moving images. This class is responsible for interconnecting
> // all the components we have described so far.
> typedef itk::ImageRegistrationMethod<
> FixedImageType,
> MovingImageType > RegistrationType;
>
> // Each one of the registration components is created using its
> // \code{New()} method and is assigned to its respective
> // \doxygen{SmartPointer}.
> //
> // Software Guide : BeginCodeSnippet
> MetricType::Pointer metric = MetricType::New();
> TransformType::Pointer transform = TransformType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
> // Software Guide : EndCodeSnippet
>
> // Software Guide : BeginLatex
> // We are using a MatchCardinalityImageToImageMetric to compare two
> // label maps. This metric simple counts the percentage of
> // corresponding pixels that have the same label. This metric does
> // not provide analytical derivatives, so we will use an
> // AmoebaOptimizer to drive the registration. The AmoebaOptimizer
> // can only minimize a cost function, so we set the metric to count
> // the percentages of mismatches.
> // Software Guide : BeginLatex
> // Software Guide : BeginCodeSnippet
> metric->MeasureMatchesOff();
> // Software Guide : EndCodeSnippet
>
>
> // Each component is now connected to the instance of the registration method.
> // \index{itk::RegistrationMethod!SetMetric()}
> // \index{itk::RegistrationMethod!SetOptimizer()}
> // \index{itk::RegistrationMethod!SetTransform()}
> // \index{itk::RegistrationMethod!SetFixedImage()}
> // \index{itk::RegistrationMethod!SetMovingImage()}
> // \index{itk::RegistrationMethod!SetInterpolator()}
> //
>
>
> 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] );
>
>
> // In this example, the fixed and moving images are read from files. This
> // requires the \doxygen{ImageRegistrationMethod} to acquire its inputs to
> // the output of the readers.
> //
> registration->SetFixedImage( fixedImageReader->GetOutput() );
> registration->SetMovingImage( movingImageReader->GetOutput() );
>
> // The registration can be restricted to consider only a particular region
> // of the fixed image as input to the metric computation. This region is
> // defined by the \code{SetFixedImageRegion()} method. You could use this
> // feature to reduce the computational time of the registration or to avoid
> // unwanted objects present in the image affecting the registration outcome.
> // In this example we use the full available content of the image. This
> // region is identified by the \code{BufferedRegion} of the fixed image.
> // Note that for this region to be valid the reader must first invoke its
> // \code{Update()} method.
> //
> // \index{itk::ImageRegistrationMethod!SetFixedImageRegion()}
> // \index{itk::Image!GetBufferedRegion()}
> //
> fixedImageReader->Update();
> movingImageReader->Update();
>
> registration->SetFixedImageRegion(
> fixedImageReader->GetOutput()->GetBufferedRegion() );
>
>
> // The parameters of the transform are initialized by passing them in an
> // array. This can be used to setup an initial known correction of the
> // misalignment. In this particular case, a translation transform is
> // being used for the registration. The array of parameters for this
> // transform is simply composed of the translation values along each
> // dimension. Setting the values of the parameters to zero
> // initializes the transform as an \emph{identity} transform. Note that the
> // array constructor requires the number of elements as an argument.
> //
> // \index{itk::TranslationTransform!GetNumberOfParameters()}
> // \index{itk::RegistrationMethod!SetInitialTransformParameters()}
> //
>
>/**********************************/
>
>
> // Display initial transform paramters
> std::cout << "Initial Transform Parameters: " << std::endl;
> std::cout << transform->GetParameters() << std::endl;
>
>
>
>
>/************************************/
> double translationScale = 1.0 / 1.000;
>
> std::cout << "begin optimizer scales" << std::endl;
> typedef OptimizerType::ScalesType OptimizerScalesType;
> OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
>
> optimizerScales[0] = 1.0;
> optimizerScales[1] = 1.0;
> optimizerScales[2] = 1.0;
> optimizerScales[3] = 1.0;
> optimizerScales[4] = 1.0;
> optimizerScales[5] = 1.0;
> optimizerScales[6] = 1.0;
> optimizerScales[7] = 1.0;
> optimizerScales[8] = 1.0;
> optimizerScales[9] = translationScale;
> optimizerScales[10] = translationScale;
> optimizerScales[11] = translationScale;
>
> optimizer->SetScales( optimizerScales );
>
> std::cout << "end optimizer scales" << std::endl;
>
> OptimizerType::ParametersType simplexDelta( transform->GetNumberOfParameters() );
> std::cout << transform->GetNumberOfParameters()<< std::endl;
> simplexDelta.Fill( 1.0 );
> std::cout << simplexDelta << std::endl;
>
> optimizer->AutomaticInitialSimplexOff();
> optimizer->SetInitialSimplexDelta( simplexDelta );
>
> // Software Guide : EndCodeSnippet
>
>
> // Software Guide : BeginLatex
> // We also adjust the tolerances on the optimizer to define
> // convergence. Here, we used a tolerance on the parameters of
> // 0.25 (which will be a quarter of image unit, in this case
> // pixels). We also set the tolerance on the cost function value to
> // define convergence. The metric we are using returns the
> // percentage of pixels that mismatch. So we set the function
> // convergence to be 0.1%
> // Software Guide : EndLatex
> // Software Guide : BeginCodeSnippet
> optimizer->SetParametersConvergenceTolerance( 0.25 ); // quarter pixel
> optimizer->SetFunctionConvergenceTolerance(0.001); // 0.1%
> // Software Guide : EndCodeSnippet
>
> // Software Guide : BeginLatex
> // In the case where the optimizer never succeeds in reaching the desired
> // precision tolerance, it is prudent to establish a limit on the number of
> // iterations to be performed. This maximum number is defined with the
> // method \code{SetMaximumNumberOfIterations()}.
> //
> // \index{itk::Amoeba\-Optimizer!SetMaximumNumberOfIterations()}
> //
> // Software Guide : EndLatex
> // Software Guide : BeginCodeSnippet
> optimizer->SetMaximumNumberOfIterations( 200 );
> // Software Guide : EndCodeSnippet
>
>
>
> //
> // Create the Command observer and register it with the optimizer.
> //
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
> registration->SetMetric( metric );
> registration->SetOptimizer( optimizer );
> registration->SetTransform( transform );
> registration->SetInterpolator( interpolator );
>
> // The registration process is triggered by an invocation of the
> // \code{StartRegistration()} method. If something goes wrong during the
> // initialization or execution of the registration an exception will be
> // thrown. We should therefore place the \code{StartRegistration()} method
> // in a \code{try/catch} block as illustrated in the following lines.
> //
> try
> {
> // print out the initial metric value. need to initialize the
>
> // registration method to force all the connections to be established.
> std::cout << "Before initialize registration" << std::endl;
> registration->Initialize();
> // std::cout << "Initial Metric value = " << metric->GetValue( transform->GetParameters() ) << std::endl;
>
> // run the registration
> std::cout << "Registration... " << std::endl;
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cout << "ExceptionObject caught !" << std::endl;
> std::cout << err << std::endl;
> return -1;
> }
>
> // In a real application, you may attempt to recover from the error in the
> // catch block. Here we are simply printing out a message and then
> // terminating the execution of the program.
> //
>
> //
> // The result of the registration process is an array of parameters that
> // defines the spatial transformation in an unique way. This final result is
> // obtained using the \code{GetLastTransformParameters()} method.
> //
> // \index{itk::RegistrationMethod!GetLastTransformParameters()}
> //
> OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
>
>
> const unsigned int numberOfIterations = optimizer->GetOptimizer()->get_num_evaluations();
>
>
> const double bestValue = metric->GetValue(finalParameters);
>
>
>
> const double versorX = finalParameters[0];
> const double versorY = finalParameters[1];
> const double versorZ = finalParameters[2];
> const double finalTranslationX = finalParameters[3];
> const double finalTranslationY = finalParameters[4];
> const double finalTranslationZ = finalParameters[5];
>
>
> // Print out results
> //
> std::cout << std::endl << std::endl;
> std::cout << "Result = " << std::endl;
> std::cout << " versor X = " << versorX << std::endl;
> std::cout << " versor Y = " << versorY << std::endl;
> std::cout << " versor Z = " << versorZ << std::endl;
> std::cout << " Translation X = " << finalTranslationX << std::endl;
> std::cout << " Translation Y = " << finalTranslationY << std::endl;
> std::cout << " Translation Z = " << finalTranslationZ << std::endl;
> std::cout << " Iterations = " << numberOfIterations << std::endl;
> std::cout << " Metric value = " << bestValue << std::endl;
>
>
> // It is common, as the last step of a registration task, to use the
> // resulting transform to map the moving image into the fixed image space.
> // This is easily done with the \doxygen{ResampleImageFilter}. Please
> // refer to Section~\ref{sec:ResampleImageFilter} for details on the use
> // of this filter. First, a ResampleImageFilter type is instantiated
> // using the image types. It is convenient to use the fixed image type as
> // the output type since it is likely that the transformed moving image
> // will be compared with the fixed image.
> //
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
>
> // A transform of the same type used in the registration process should be
> // created and initialized with the parameters resulting from the
> // registration process.
> //
> // \index{itk::ImageRegistrationMethod!Resampling image}
> //
> TransformType::Pointer finalTransform = TransformType::New();
> finalTransform->SetParameters( finalParameters );
>
> // Then a resampling filter is created and the corresponding transform and
> // moving image connected as inputs.
> //
> ResampleFilterType::Pointer resample = ResampleFilterType::New();
> resample->SetTransform( finalTransform );
> resample->SetInput( movingImageReader->GetOutput() );
>
> // As described in Section \ref{sec:ResampleImageFilter}, the
> // ResampleImageFilter requires additional parameters to be
> // specified, in particular, the spacing, origin and size of the output
> // image. The default pixel value is also set to the standard label
> // for "unknown" or background. Finally, we need to set the
> // interpolator to be the same type of interpolator as the
> // registration method used (nearest neighbor).
> //
> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
> resample->SetOutputOrigin( fixedImage->GetOrigin() );
> resample->SetOutputSpacing( fixedImage->GetSpacing() );
> //resample->SetDefaultPixelValue( 0 );
> resample->SetInterpolator( interpolator );
>
> // The output of the filter is passed to a writer that will store the
> // image in a file. An \doxygen{CastImageFilter} is used to convert the
> // pixel type of the resampled image to the final type used by the
> // writer. The cast and writer filters are instantiated below.
> //
> typedef float OutputPixelType;
> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
> typedef itk::CastImageFilter<
> FixedImageType,
> OutputImageType > CastFilterType;
> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>
> // The filters are created by invoking their \code{New()}
> // method.
> //
> WriterType::Pointer writer = WriterType::New();
> CastFilterType::Pointer caster = CastFilterType::New();
>
>
> writer->SetFileName( argv[3] );
>
>
> // The \code{Update()} method of the writer is invoked in order to trigger
> // the execution of the pipeline.
> //
> caster->SetInput( resample->GetOutput() );
> writer->SetInput( caster->GetOutput() );
> std::cout << std::endl <<"Writing image... :" << argv[3] <<std::endl;
> writer->Update();
>
>
> //
> // The fixed image and the transformed moving image can easily be compared
> // using the \code{SquaredDifferenceImageFilter}. This pixel-wise
> // filter computes the squared value of the difference between homologous
> // pixels of its input images.
> //
>
>
> return 0;
>}
>
>// SoftwareGuide : BeginLatex
>// The example was run on two binary images. The first binary image was generated by running the
>// confidence connected image filter (section \ref{sec:ConfidenceConnected}) on
>// the MRI slice of the brain. The second was generated similarly after
>// shifting the slice by 13 pixels horizontally and 17 pixels
>// vertically. The Amoeba optimizer converged after 34 iterations
>// and produced the following results:
>//
>// \begin{verbatim}
>// Translation X = 12.5
>// Translation Y = 16.77
>// \end{verbatim}
>// These results are a close match to the true misalignment.
>// SoftwareGuide : EndLatex
>
>
>------------------------------------------------------------------------
>
>_______________________________________________
>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