[Insight-users] Re: Problems in changing the optimisation

[Luis Ibanez]

This is a multi-part message in MIME format.
--------------000109070402040304020901
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit



Hi Shekhar,

Comparing the performance of several optimizers
is certainly an interesting exercise. However, you
shouldn't expect all optimizers to be appropriate
for a particular optimization problem.

The choice of the optimizer and its parameters
depends on the characteristics of the image Metric
in the parametric space of the Transform.

I would suggest you to start with a simple
Translation transform instead of the AffineTransform.
When you use an AffineTransform, you have to deal
with the optimizationScales, that are probably the
most critical parameters in the entire optimization
problem.

Also, do not stat this exercise in 3D. It will be
unnecessarily long, you will need something faster
where you can run many tests in an hour.  Please
use a 2D images, 128x128 will be fine for the size.

Now... before you start playing with the optimizers,
please generate plot of the metric, so you know how
the Cost function looks like. The attached program,
illustrates how to do this. Take the output of the
program and plot the metric with something like
GNUplot or even Excel.  The aspect of this plot will
give you good indications on what optimizer could do
a good job and what its parameters may be.

Then, every time you test an optimizer, make sure
that you connect a Command/Observer to it. This
will allow you to get a trace of the transform
parameters and the metric found by the optimizer.
Take that trace and plot it on top of the Metric
graph that you got from the attached program.

The path followed by the optimizer in this plot
will indicate you if the Metric is too noisy,
or if the convergence criteria that you are using
is too stringent, or if the step length is too
large or too small.

Note that the optimizers may stop for several
different reasons. Check what lead them to stop.
Then you will know if you are being too restrictive
with any of the convergence criteria: number of
iterations, gradient magnitude, minimum step length...


When you compare the OnePlusOne optimizer, beware
that this optimizer only change is position when
it found a better spot. In the meantime it is
evaluating many location around its current position.
you may want to put a cout in the Metric itself
in order to have a better perception of the full
activity that is going on inside the optimizer.



Please let us know if you have further questions.



    Thanks


      Luis



------------------------------------------
Dwivedi, Shekhar (Research, Non-ge) wrote:

> hi luis,
> 
> I am using MultiResImageRegistration...  
> 
> On changing the optimizer(Regular step gradient) to any of the following :
> 
>  Amoeba,OneplusOneEvolutionary,LBFGS etc
> 
> the results are not proper (as compared to the results from Regular step gradient).
> for LBFGS , I am getting the error as too many samples mapped ouside the image, for others the result is incorrect.
> 
> the only change I have made is in calling the parameters of the optimizer,it is also not clear which 
> parameters are required to be called for a particular optimizer. As in Regular Step Gradient only two parameters
> are set.
> 
> For amoeba I have set :
> 
> Maximum no.of iterations : 100
> Parameter Convergence Tolerance : .0001
> Function Convergence Tolerance : .01	
> 
> if(registration->GetCurrentLevel() == 0)
> {
> 	optimizer->Set....
> }
> else
> {
> 	optimizer->Set..( optimizer->Get..() )
> }
> 
> For LBFGS : 
> 
> Maximum no.of function evaluations : 100
> Gradient Convergence Tolerance : .1
> Line Search Accuracy = 0.1
> Default Step Length  = 4.0
> 
> and have called them as in amoeba above.
> 
> Similarly I have called the parameters for few other optimizers also but m not getting the correct result.
> I am using Affine Transform , but I think if the same transform is working for the regular step gradient then
> it should hold true for other also (at least few).
> 
> I am also not clear whether changes to this particular file is only required.
> 
> The image which I am using is of size 128*128*128
> with pixel spacing of 1mm.
> 
> Please let me know the approach required to select the best optimizers.
> 
> Thank you.
> 
> Shekhar.
> 


--------------000109070402040304020901
Content-Type: text/plain;
 name="meanSquaredDifferences.cxx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="meanSquaredDifferences.cxx"


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkTranslationTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"


int main( int argc, char * argv[] )
{
  if( argc < 3 )
    {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << "  fixedImage  movingImage" << std::endl;
    return 1;
    }

  const     unsigned int   Dimension = 2;
  typedef   unsigned char  PixelType;

  typedef itk::Image< PixelType, Dimension >   ImageType;
  typedef itk::Image< PixelType, Dimension >   ImageType;


  typedef itk::ImageFileReader< ImageType >  ReaderType;

  ReaderType::Pointer fixedReader  = ReaderType::New();
  ReaderType::Pointer movingReader = ReaderType::New();

  fixedReader->SetFileName(  argv[ 1 ] );
  movingReader->SetFileName( argv[ 2 ] );

  try 
    {
    fixedReader->Update();
    movingReader->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception catched !" << std::endl;
    std::cerr << excep << std::endl;
    }


  typedef itk::MeanSquaresImageToImageMetric< ImageType, ImageType >  MetricType;

  MetricType::Pointer metric = MetricType::New();



  typedef itk::TranslationTransform< double, Dimension >  TransformType;

  TransformType::Pointer transform = TransformType::New();



//  typedef itk::LinearInterpolateImageFunction< 
  typedef itk::NearestNeighborInterpolateImageFunction< 
                                 ImageType, double >  InterpolatorType;

  InterpolatorType::Pointer interpolator = InterpolatorType::New();


  metric->SetInterpolator( interpolator );
  metric->SetTransform( transform );


  transform->SetIdentity();

  ImageType::ConstPointer fixedImage  = fixedReader->GetOutput();
  ImageType::ConstPointer movingImage = movingReader->GetOutput();

  metric->SetFixedImage(  fixedImage  );
  metric->SetMovingImage( movingImage );

  metric->SetFixedImageRegion(  fixedImage->GetBufferedRegion()  );

  try 
    {
    metric->Initialize();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception catched !" << std::endl;
    std::cerr << excep << std::endl;
    return -1;
    }


  MetricType::TransformParametersType displacement( Dimension );

  int rangex = 50;
  int rangey = 50;

  for( int dx = -rangex; dx <= rangex; dx++ )
    {
    for( int dy = -rangey; dy <= rangey; dy++ )
      {
      displacement[0] = dx;
      displacement[1] = dy;
      const double value = metric->GetValue( displacement );

      std::cout << dx << "   "  << dy << "   " << value << std::endl;
      }
    }

  std::cout << std::endl;


  return 0;
}


--------------000109070402040304020901--