[Insight-users] Mattes MI & 1+1 ES optimizer = overlap exception

[Christos Panagiotou]

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