[Insight-users] registration example 2 and 9
Luis Ibanez
luis.ibanez at kitware.com
Thu Nov 18 09:08:21 EST 2004
Hi Lagaffe,
1) Yes, it is recommended that you normalize your images
before passing them to the MutualInformation metric.
This normalization can be done with the NormalizeImageFilter
http://www.itk.org/Insight/Doxygen/html/classitk_1_1NormalizeImageFilter.html
as it is illustrated in the ITK Software Guide
http://www.itk.org/ItkSoftwareGuide.pdf
in section 8.4.1, pdf-page 255.
2) Did you called update in the smoothing filters *BEFORE*
you execute the Transform initializer ?
As we mention in the previous email, the Transform
initializer *is not* a pipeline object. You must pass
to it the output of filters that have been updated.
3) It is fundamental that you solve the initialization
problem first. There is nothing that the registration
method can do if you are still feeding it with initial
parameters such as the ones you posted in your previous
email:
> > ---------------
> > Initial Transform Parameters
> > [1, 0, 0, 1, -1.00276e+008, -4.25387e+008, 4.90708e+007,
> 3.85172e+008]
One good way to verify your initialization is to skip the full
registration process (by now) and use the *initial transform*
for resampling the moving image, then save it in a file. If that
image doesn't look close to the fixed image, then you are not
ready yet for starting your registration process.
Please let us know if you have further questions,
Thanks
Luis
-----------------------------------------------
-----------------
Mr Gaffe wrote:
> Thanks Luis,
> Sorry, but If you look again at the code you will see that I dit the
> initialization of CenteredTransformInitializer with the normalize and
> smooth images.
> initializer->SetFixedImage( fixedSmoother->GetOutput() );
> initializer->SetMovingImage( movingSmoother->GetOutput() );
> But if I replace these two lines with the orginals images:
> initializer->SetFixedImage( fixedImageReader->GetOutput() );
> initializer->SetMovingImage( movingImageReader->GetOutput() );
> I do not have an exception but it does not converge correctly.
>
> And I have to normalize and smooth my images to correctly use the
> Mutual information has you done in example 2....
>
> Just for information these are the value if I initialize without
> normalize and smooth the image. with the initial fixed and moving image.
> Initial Transform Parameters
> [1, 0, 0, 1, 111.204, 131.591, -12.7911, -12.6914]
> 0 0.360832 [0.993629, -0.00122832, 0.00831941, 1.01269, 111.204,
> 131.591, -12.7794, -12.5935]
> 1 0.220862 [0.991423, -0.00269869, 0.00381388, 1.00791, 111.204,
> 131.591, -12.7428, -12.6268]
> 2 0.458038 [0.991053, -0.00423641, 0.00555051, 1.0107, 111.204,
> 131.591, -12.7455, -12.6022]
> 3 0.43253 [0.98915, -0.00600989, 0.00460683, 1.00976, 111.204,
> 131.591, -12.7564, -12.6076]
> 4 0.391366 [0.988866, -0.0053906, 0.00362396, 1.0088, 111.204,
> 131.591, -12.7547, -12.6199]
> 5 0.255813 [0.988056, -0.00601446, 0.00374833, 1.0088, 111.204,
> 131.591, -12.7609, -12.6203]
> 6 0.23651 [0.987756, -0.00630565, 0.00309831, 1.00838, 111.204,
> 131.591, -12.7629, -12.6261]
> 7 0.312396 [0.987953, -0.00558894, 0.00279363, 1.008, 111.204,
> 131.591, -12.7579, -12.6298]
> 8 0.157922 [0.987289, -0.00548913, 0.00273997, 1.00797, 111.204,
> 131.591, -12.7529, -12.6261]
> 9 0.259104 [0.986384, -0.0058869, 0.00200088, 1.00744, 111.204,
> 131.591, -12.7518, -12.6287]
> 10 0.37124 [0.98628, -0.00594178, 0.00209053, 1.00749, 111.204,
> 131.591, -12.753, -12.6278]
> Result =
> Center X = 111.204
> Center Y = 131.591
> Translation X = -12.753
> Translation Y = -12.6278
> Iterations = 12
> Metric value = 0.520743
>
>
>
> thanks again for your help
> Lagaffe
>
> */Luis Ibanez <luis.ibanez at kitware.com>/* wrote:
>
>
>
> 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 ob! ject 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->SetFixedIma! ge( 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 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 * obje! ct, 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->SetNumberOfIt! erations( 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
> >
> >
> > *Le nouveau Yahoo! Messenger est arrivé !* Découvrez toutes les
> > nouveautés pour dialoguer instantanément avec vos amis. Téléchargez
> > GRATUITEMENT ici !
> >
> >
> >
> >
> >
> ------------------------------------------------------------------------*>
>
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users at itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
>
>
>
>
> **
>
> ** **
>
> ------------------------------------------------------------------------
> **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