[Insight-users] Re: Mattes MI & 1+1 ES optimizer = overlap exception
Luis Ibanez
luis.ibanez at kitware.com
Mon, 03 May 2004 11:20:19 -0400
Hi Christos,
Thanks for your detailed message,
A) Number of samples: Since the metric needs to estimate
the Joint histogram of the two image, it probably makes
sense to use the 1-20% of the larger image. If you
apply the ratio to the smaller image, it is likely that
the number of sampled points will not be enough for being
representative of the intensity distribution in the larger
image.
B) About the exception you are getting:
" Too many samples map outside moving
image buffer: 50703 / 360448 "
It means that the Affine transform that maps the moving
image is not generating enough overlap with the fixed
image.
I couldn't find in your code, where you are initializing
the AffineTransform. I would have expected something like:
m_AffineTransform->SetIdentity();
(maybe I missed it...)
Even better than that, I would strongly suggest you to use
the CenteredAffineTransform and the CenteredTransformInitializer.
This initializer will take care of matching the image centers,
by taking into account their respective origin, pixel spacing
and number of pixels along every dimension.
The fact that you are not even getting the first GetValue()
result, lead us to suspect that you simply have a poor initialized
transform. Note that Identity is not neecesarily a good
initialization. It depends on the origin and pixel spacing of
both of your images. A fundamental exercise that you should
perform before even starting your registration is to draw in
a piece of paper the physical extent of both images in a common
coordinate system. This will give you a clear idea of what
your initialization should be for the transform. This is illustrated
multiple times in the "resampling" section of the SoftwareGuide.
A similar effect can be obtained by simply taking a resample
image filter, feeding it with an itk::IdentityTransform and
mapping the Moving image into the Fixed image space. That
will give you a clear view of the current spatial relationship
between the two images.
C) It is not worth to start the registration until you find a
reasonable initialization for the transform...
Regards,
Luis
---------------------------
Christos Panagiotou wrote:
> Dear Luis
>
> i updated the one + one es optimizer from the cvs... thanks for letting
> me know about the oneplusone optimizer example (ImageRegistration11.cxx)
>
> i try to test the optimizer and the mattes metric on a registration of
> one 256x256x55 mri volume against a 18x128x31 pet volume...
>
> Mattes histogram metric
> 1.
> well i choose the translation scale to be the diagonial of the biggest
> volume in mm. So my spacing of the MRI volume is 1.0x1.0x1.5 ->
> 1.0 / sqrt( 256^2 + 256^2 + (256*1.5)^2) * 10
>
> -- you have pointed out that i need to multiply by the *magic factor* 10.
>
> 2.
> according to the itk guide the number of spatial samples in the metric
> should be 1 - 20 % of the pixels (pixels of both images? or the largest?)
> lets say 256*256*55 = 3604480 where 10% is 360448 spatial samples
>
> 3.
> number of bins = 50 (default )
>
> 4.
> 1+1 ES optimizer
>
> initialization of m_Optimizer->SetNormalVariateGenerator(m_Generator);
> where m_Generator->Initialize(12345);
>
> Epsilon : i ve tried 0.000001 - 1.0
> Initial search radius: i ve tried 10 - 30
> GrowthFactor: 1.05 - 1.30
> ShrinkFactor: 0.80 - 0.98
>
>
>
> I have modified the MultiResMIRegistration application so i can pass the
> optimizer values through a .txt file ( i ve modified actually
> SimpleAppInputParser.txx, .h mainly, default was for gradient descent
> optimizer as described in the itk guide) so i can pass the 4 parameters
> of the 1+1 optimizer (epsilon, initial radius, shrink , growth). I tell
> you this so i can remind you that
> the MultiResMiRegistration CENTERS the volumes before the registration.
> So some overlap definately exists.
>
> Let me post some results:
>
> ./MultiResMIRegistration input_mriAdult_pet_1plus1ES.txt
> Parsing input ...
> Fixed image filename: mr_adult.mha
> Big Endian: 0
> Image Size: [256, 256, 55]
> Image Spacing: [1, 1, 1.5]Moving image filename: pet1.mha
> Big Endian: 0
> Image Size: [128, 128, 31]
> Image Spacing: [1, 1, 1]Permute order: [0, 1, 2]Flip axes: [0, 0,
> 0]Number of levels: 1
> Fixed image shrink factors: [1, 1, 1]Moving image shrink factors: [1, 1,
> 1]Number of iterations: [4000]Initial Search Radius: 15
> Epsilon: 1e-06
> Growth Factor: 1.3
> Shrink Factor: 0.7
> Translation scale: 371.32
> PGM directory: pettopd
>
> Preprocess the images ...
> Register the images ...
> --- Starting level 0
> No. Iterations: 4000
> Caught an exception:
>
> itk::ExceptionObject (0x8cb3800)
> Location: "Unknown"
> File:
> /usr/local/include/InsightToolkit/Algorithms/itkMattesMutualInformationImageToImageMetric.txx
>
> Line: 546
> Description: itk::ERROR:
> MattesMutualInformationImageToImageMetric(0x8cb0410): Too many samples
> map outside moving image buffer: 50703 / 360448
>
>
> even if i dramaticaly change the parameters i get the same ration exactly
>
> Initial Search Radius: 40
> Epsilon: 0.001
> Growth Factor: 1.7
> Shrink Factor: 0.3
> Translation scale: 371.32
>
> It seems somehow that something is not initialized well as the system
> does not produce and ouput from the execute method... (no
> m_Optimizer->GetValue() is returned)
> the algorithm throws the exception immediately...
>
> well i dont know whats wrong....
>
> i know that this post is already long enough... i ll just paste the
> modified <MIMRegistrator.txx>
>
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: MIMRegistrator.txx,v $
> Language: C++
> Date: $Date: 2002/10/01 01:36:06 $
> Version: $Revision: 1.3 $
>
> Copyright (c) 2002 Insight 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.
>
> =========================================================================*/
> #ifndef _MIMRegistrator_txx
> #define _MIMRegistrator_txx
>
> #include "MIMRegistrator.h"
>
> #include "itkCommand.h"
>
>
>
> namespace itk
> {
>
>
>
>
> template <typename TFixedImage, typename TMovingImage>
> MIMRegistrator<TFixedImage,TMovingImage>
> ::MIMRegistrator()
> {
> // Images need to be set from the outside
> m_FixedImage = NULL;
> m_MovingImage = NULL;
>
>
> // Set up internal registrator with default components
> m_Transform = TransformType::New();
> m_Optimizer = OptimizerType::New();
> m_Metric = MetricType::New();
> m_Interpolator = InterpolatorType::New();
> m_FixedImagePyramid = FixedImagePyramidType::New();
> m_MovingImagePyramid = MovingImagePyramidType::New();
> m_Registration = RegistrationType::New();
> CommandIterationUpdate::Pointer m_Observer;
> m_Observer = CommandIterationUpdate::New();
>
> m_Generator = OptimizerNormalGeneratorType::New();
>
> m_Generator->Initialize(12345);
>
> m_Registration->SetTransform( m_Transform );
> m_Registration->SetOptimizer( m_Optimizer );
> m_Registration->SetMetric( m_Metric );
> m_Registration->SetInterpolator( m_Interpolator );
> m_Registration->SetFixedImagePyramid( m_FixedImagePyramid );
> m_Registration->SetMovingImagePyramid( m_MovingImagePyramid );
>
>
>
> m_Optimizer->AddObserver(IterationEvent(), m_Observer );
>
> m_AffineTransform = AffineTransformType::New();
>
> // Setup an registration observer
>
> typedef SimpleMemberCommand<Self> CommandType;
> typename CommandType::Pointer command = CommandType::New();
> command->SetCallbackFunction( this, &Self::StartNewLevel );
> m_Tag = m_Registration->AddObserver( IterationEvent(), command );
>
> // Default parameters
> m_NumberOfLevels = 1;
> m_TranslationScale = 1.0;
> m_MovingImageStandardDeviation = 0.4;
> m_FixedImageStandardDeviation = 0.4;
> m_NumberOfSpatialSamples = 50;
>
> m_FixedImageShrinkFactors.Fill( 1 );
> m_MovingImageShrinkFactors.Fill( 1 );
>
> m_NumberOfIterations = UnsignedIntArray(1);
> //Set elements of the array to the specified value
> m_NumberOfIterations.Fill( 1000 );
>
> m_InitialParameters = ParametersType(m_Transform->GetParameters());
>
>
> }
>
> template <typename TFixedImage, typename TMovingImage>
> MIMRegistrator<TFixedImage,TMovingImage>
> ::~MIMRegistrator()
> {
> m_Registration->RemoveObserver( m_Tag );
>
> }
>
>
> template <typename TFixedImage, typename TMovingImage>
> void
> MIMRegistrator<TFixedImage,TMovingImage>
> ::Execute()
> {
>
> // Setup the optimizer
> typename OptimizerType::ScalesType scales(
> m_Transform->GetNumberOfParameters() );
> scales.Fill(1.0);
>
> for ( int j = 9; j < 12; j++ )
> {
> scales[j] = m_TranslationScale;
> }
>
> m_Optimizer->SetScales( scales );
>
> m_Optimizer->MaximizeOn();
>
> m_Optimizer->Initialize(15); // Initial search radius
> m_Optimizer->SetEpsilon(m_Epsilon);
> m_Optimizer->SetGrowthFactor(m_GrowthFactor);
> m_Optimizer->SetShrinkFactor(m_ShrinkFactor);
>
> m_Optimizer->SetNormalVariateGenerator(m_Generator);
>
>
>
> m_Metric->SetNumberOfHistogramBins( 50 );
> m_Metric->SetNumberOfSpatialSamples( 360448 );
>
>
> // Setup the image pyramids
> m_FixedImagePyramid->SetNumberOfLevels( m_NumberOfLevels );
> m_FixedImagePyramid->SetStartingShrinkFactors(
> m_FixedImageShrinkFactors.GetDataPointer() );
>
> m_MovingImagePyramid->SetNumberOfLevels( m_NumberOfLevels );
> m_MovingImagePyramid->SetStartingShrinkFactors(
> m_MovingImageShrinkFactors.GetDataPointer() );
>
> // Setup the registrator
> m_Registration->SetFixedImage( m_FixedImage );
> m_Registration->SetMovingImage( m_MovingImage );
> m_Registration->SetNumberOfLevels( m_NumberOfLevels );
>
> m_Registration->SetInitialTransformParameters(
> m_Transform->GetParameters() );
>
> m_Registration->SetFixedImageRegion( m_FixedImage->GetBufferedRegion() );
>
> try
> {
> m_Registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cout << "Caught an exception: " << std::endl;
> std::cout << err << std::endl;
> throw err;
> }
>
> }
>
>
> template <typename TFixedImage, typename TMovingImage>
> const
> typename MIMRegistrator<TFixedImage,TMovingImage>
> ::ParametersType &
> MIMRegistrator<TFixedImage,TMovingImage>
> ::GetTransformParameters()
> {
> return m_Registration->GetLastTransformParameters();
> }
>
>
> template <typename TFixedImage, typename TMovingImage>
> typename MIMRegistrator<TFixedImage,TMovingImage>
> ::AffineTransformPointer
> MIMRegistrator<TFixedImage,TMovingImage>
> ::GetAffineTransform()
> {
> m_Transform->SetParameters(
> m_Registration->GetLastTransformParameters() );
>
>
> m_AffineTransform->SetMatrix( m_Transform->GetMatrix() );
> m_AffineTransform->SetOffset( m_Transform->GetOffset() );
>
> return m_AffineTransform;
> }
>
>
>
> template <typename TFixedImage, typename TMovingImage>
> void
> MIMRegistrator<TFixedImage,TMovingImage>
> ::StartNewLevel()
> {
> std::cout << "--- Starting level " << m_Registration->GetCurrentLevel()
> << std::endl;
>
> unsigned int level = m_Registration->GetCurrentLevel();
> if ( m_NumberOfIterations.Size() >= level + 1 )
> {
> m_Optimizer->SetMaximumIteration( m_NumberOfIterations[level] );
> }
>
> std::cout << " No. Iterations: "
> << m_NumberOfIterations[level]
> << std::endl;
>
> }
>
>
>
> } // namespace itk
>
>
> #endif
>
>
> i apologise once again for the length of the post, but i would like to
> see this working someday :)
> thanks for your patience
>
> christos
>
>
>