[Insight-users] registration example 2 and 9
Luis Ibanez
luis.ibanez at kitware.com
Wed Nov 10 09:50:11 EST 2004
Hi Lagaffe,
The error is likely related to the fact that you commented
out the lines where the input images are connected to the
CenteredTransformInitializer
>
> /*
> initializer->SetFixedImage( fixedImageReader->GetOutput() );
> initializer->SetMovingImage( movingImageReader->GetOutput() );
> */
>
Therefore, the initialization of the Transform is not being
performed correctly and you start your optimization process
with corrupted numerical values such as:
> ---------------
> Initial Transform Parameters
> [1, 0, 0, 1, -1.00276e+008, -4.25387e+008, 4.90708e+007, 3.85172e+008]
>
Also, note that you *MUST* call Update() in the readers before you use
the TransformInitializer object. The reason is that the Initializer is
not a pipeline object and therefore it will not trigger an update in its
inputs. This object is similar to the image calculator in this regard.
If you notice that the optimizer is stopping too early without
converging, that indicates that you are using:
A) a value too large for the step length of the optimizer
B) inappropriate values for the parameter scaling
If you continue having problems while fine tunning those parameters
please post to the list the log output of the Command/Observer.
From the plot of the Metric values and the Transform parameters at
every iteration it is usually easier to figure out what parameters
must be changed in your optimization setup.
Please let us know if you have further questions,
Thanks,
Luis
---------------
Mr Gaffe wrote:
> Hi all,
> I try to combine registration example: ImageRegistration2 and 9 to
> perform an CenteredAffine+MutualInfo registration.
> The code bellow compile well, but when I try with:
> BrainT1SliceBorder20.png and BrainProtonDensitySliceR10X13Y17S12.png I
> have the following error:
>
> ---------------
> Initial Transform Parameters
> [1, 0, 0, 1, -1.00276e+008, -4.25387e+008, 4.90708e+007, 3.85172e+008]
> ExceptionObject caught !
> itk::ExceptionObject (00FCFD30)
> Location: "Unknown"
> File:
> D:\InsightToolkit-1.8.1\Code\Algorithms\itkMutualInformationImageToImageMetric.txx
> Line: 155
> Description: itk::ERROR: MutualInformationImageToImageMetric(012F4780):
> All the sampled point mapped to outside of the moving image
> -----------------
>
> ->The Typedef is set to float since Images are unsigned chart is it the
> problem ?
> If I change the initialisation of the CenteredTransformInitializer as:
> initializer->SetFixedImage( fixedImageReader->GetOutput() );
> initializer->SetMovingImage( movingImageReader->GetOutput() );
> There is no error, but the after 9 iterations it stop with a wrong image
> results....
>
> Is it just a problem of initial typeddef for my fixed and moved image ?
>
> thanks,
>
> lagaffe
>
> ---------------------------------------------------------------------------------------------------------------------
> /*
> * Code from ImageRegistration2 and ImageRegistration9
> */
> #include "itkImageRegistrationMethod.h"
> //#include "itkMeanSquaresImageToImageMetric.h"
> #include "itkMutualInformationImageToImageMetric.h"
> #include "itkNormalizeImageFilter.h"
> #include "itkDiscreteGaussianImageFilter.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkRegularStepGradientDescentOptimizer.h"
> #include "itkImage.h"
>
> #include "itkCenteredTransformInitializer.h"
> #include "itkCenteredAffineTransform.h"
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkSquaredDifferenceImageFilter.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::RegularStepGradientDescentOptimizer 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::endl;
> std::cerr << " outputImagefile [differenceOutputfile]
> [differenceBeforeRegistration] " << std::endl;
> std::cerr << " [stepLength] [maxNumberOfIterations] "<< std::endl;
> return 1;
> }
> const unsigned int Dimension = 2;
> typedef float PixelType;
> typedef float InternalPixelType;
> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
> typedef itk::Image< PixelType, Dimension > FixedImageType;
> typedef itk::Image< PixelType, Dimension > MovingImageType;
>
> typedef itk::CenteredAffineTransform<
> double,
> Dimension > TransformType;
> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
> /*
> typedef itk::MeanSquaresImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;
> */
> typedef itk::MutualInformationImageToImageMetric<
> InternalImageType,
> InternalImageType > MetricType;
> typedef itk:: LinearInterpolateImageFunction<
> MovingImageType,
> double > InterpolatorType;
> typedef itk::ImageRegistrationMethod<
> FixedImageType,
> MovingImageType > RegistrationType;
> MetricType::Pointer metric = MetricType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
> TransformType::Pointer transform = TransformType::New();
>
> registration->SetMetric( metric );
> registration->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator );
> registration->SetTransform( transform );
> 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] );
> 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() );
> registration->SetFixedImage( fixedSmoother->GetOutput() );
> registration->SetMovingImage( movingSmoother->GetOutput() );
>
> fixedNormalizer->Update();
> registration->SetFixedImageRegion(
> fixedNormalizer->GetOutput()->GetBufferedRegion() );
>
> typedef itk::CenteredTransformInitializer<
> TransformType,
> FixedImageType,
> MovingImageType >
> TransformInitializerType;
> TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
> initializer->SetTransform( transform );
> /*
> initializer->SetFixedImage( fixedImageReader->GetOutput() );
> initializer->SetMovingImage( movingImageReader->GetOutput() );
> */
>
> initializer->SetFixedImage( fixedSmoother->GetOutput() );
> initializer->SetMovingImage( movingSmoother->GetOutput() );
> initializer->MomentsOn();
> initializer->InitializeTransform();
> std::cout << "Initial Transform Parameters " << std::endl;
> std::cout << transform->GetParameters() << std::endl;
>
> // -----------------------
> registration->SetInitialTransformParameters(
> transform->GetParameters() );
>
> double translationScale = 1.0 / 1000.0;
> if( argc > 8 )
> {
> translationScale = atof( argv[8] );
> }
> 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] = translationScale;
> optimizerScales[5] = translationScale;
> optimizerScales[6] = translationScale;
> optimizerScales[7] = translationScale;
> optimizer->SetScales( optimizerScales );
> double steplength = 0.1;
> if( argc > 6 )
> {
> steplength = atof( argv[6] );
> }
>
> unsigned int maxNumberOfIterations = 300;
> if( argc > 7 )
> {
> maxNumberOfIterations = atoi( argv[7] );
> }
> optimizer->SetMaximumStepLength( steplength );
> optimizer->SetMinimumStepLength( 0.001 );
> optimizer->SetNumberOfIterations( maxNumberOfIterations );
> optimizer->MinimizeOn();
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
> OptimizerType::ParametersType finalParameters =
> registration->GetLastTransformParameters();
> const double finalRotationCenterX = finalParameters[4];
> const double finalRotationCenterY = finalParameters[5];
> const double finalTranslationX = finalParameters[6];
> const double finalTranslationY = finalParameters[7];
> const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
> const double bestValue = optimizer->GetValue();
> // Software Guide : EndCodeSnippet
>
> // Print out results
> //
> std::cout << "Result = " << 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() );
> 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();
>
> writer->SetFileName( argv[3] );
>
> caster->SetInput( resample->GetOutput() );
> writer->SetInput( caster->GetOutput() );
> writer->Update();
>
> typedef itk::SquaredDifferenceImageFilter<
> FixedImageType,
> FixedImageType,
> OutputImageType > DifferenceFilterType;
> DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
> WriterType::Pointer writer2 = WriterType::New();
> writer2->SetInput( difference->GetOutput() );
>
> if( argc > 4 )
> {
> difference->SetInput1( fixedImageReader->GetOutput() );
> difference->SetInput2( resample->GetOutput() );
> writer2->SetFileName( argv[4] );
> writer2->Update();
> }
>
> if( argc > 5 )
> {
> writer2->SetFileName( argv[5] );
> difference->SetInput1( fixedImageReader->GetOutput() );
> difference->SetInput2( movingImageReader->GetOutput() );
> writer2->Update();
> }
>
> return 0;
> }
>
>
> ------------------------------------------------------------------------
> Créez gratuitement votre Yahoo! Mail avec *100 Mo de stockage !*
> Créez votre Yahoo! Mail
> <http://fr.rd.yahoo.com/mail/taglines/*http://fr.rd.yahoo.com/evt=25917/*http://fr.rd.yahoo.com/mail/mail_taglines_100/default/*http://fr.benefits.yahoo.com/>
>
> *Le nouveau Yahoo! Messenger est arrivé !* Découvrez toutes les
> nouveautés pour dialoguer instantanément avec vos amis. Téléchargez
> GRATUITEMENT ici !
> <http://fr.rd.yahoo.com/mail/taglines/*http://fr.rd.yahoo.com/evt=26111/*http://fr.rd.yahoo.com/messenger/mail_taglines/default/*http://fr.messenger.yahoo.com>
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> 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