[Insight-users] registration example 2 and 9

Mr Gaffe lagaffe74130 at yahoo.fr
Wed Nov 10 10:58:45 EST 2004


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 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 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 
> 
> 
> *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

Le nouveau Yahoo! Messenger est arrivé ! Découvrez toutes les nouveautés pour dialoguer instantanément avec vos amis.Téléchargez GRATUITEMENT ici !
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20041110/116fdb66/attachment.html


More information about the Insight-users mailing list