[Insight-users] Mutual Info-based Registration on Ultrasound image

[Luis Ibanez]

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


Hi Thomas,

Let me start by saying that we all have a lot of respect
for those who dare to work with Ultra Sound images.  :-)

Ultrasound is more challenging for Mutual Information due
to the presence of speckle/noise. Since MI uses a reduced
number of samples (typically 50 to 100), there are good
chances that the Metric values will have a large variability
because samples will fall sometimes in bright speckle spots
and sometimes in dark speckle spot.

I will suggest you to start by using the Mattes implementation
of Mutual Information. This Metric produces a less noisy output
than the Viola-Wells implementation.  This will not yet solve
your problem though...


Are you trying to register UltraSound to UltraSound ?
or UltraSound to MRI ?

This will be very important in order to chose the metric and
anticipate its behavior.




In order to get some insight on how to tune the parameters
for your optimization, probably the best thing to do is to
study the characteristics of the Metric when computed over
your images.  Please find attached a small program  that allows
you to explore the values of the metric using a translation
transform. You want to plot these values in order to get a
feeling on the level of noise of the Metric, the presence
of local minima, and the overall reproducibility of the Metric
function itself. (you will have to convert the dimension
to 3D, since the example was configured for 2D)

Please run this program first using the *same* input image
as fixed and moving images, and let the translation be
evaluated for a range that is close to 1/2 of the image extent
(physical extent in millimeters).

Plot the values using your favorite plotting program. (anything
from GNUPlot to Excel).

For examples on how these plot will look like, please take
a  look at the course in Image Registration:

     http://www.cs.rpi.edu/courses/spring04/imagereg/

     In particular, look at lectures 8 and 9.


Once you get this plot, I'll suggest you to experiment with
the edge preserving smoothing filters described in the Software
Guide. Run these filter as pre-processing of the input images,
and then repeat the metric plots. Since the filters will reduce
the noise, and hopefully some of the US speckle, it is to expect
that the Metric plot will also be smoothed... and therefore
easier to optimize.

You could also give it a try to more traditional smoothing such
as Gaussians....

---

It is not productive to fight with the registration + optimization
until you find a way of generating relatively smooth Metric plots.
Note that this is just an exercise on learning how to tune the
parameters.  Once you figure out the parameters, you will not need
to plot the Metric anymore.

Make sure that origin and spacing are correctly set in your images
before you start computing all these metric plots.



   Please let us know if you have further questions.


     Thanks


        Luis


-------------------------
Thomas - Kuiran Chen wrote:

> Hi Luis, 
> 
> I've been trying to use ITK's implementation of viola's mutual information based registration (itk::MutualInformationImageToImageMetric) for ultrasound images.
> 
> The parameter settings of the matric in the example (ImageRegistration2) work well on MRI image, but when I tried with ultrasound images (I even tried two identical images), the result was interesting but somehow unexpected.  The bestValue from the optimizer was around 0.5 (mutual info entropy), and the program suggested a translation of (0.11, 0.09) along X and Y directions, and the number of iterations used was 200. 
> 
> The parameters need to be set are (as for the MRI image - BrainT1SliceBorder20.pgn, in the example):
> 
> 1. for itk::MutualInformationImageToImageMetric:
> - Fixed and moving image standard deviation, the example suggests 0.4;
> - Number of Spatial Sample Bins, the example uses 50;
> 
> 2. for itk::DiscreteGaussianImageFilter:
> - Fixed and moving image smoother various, the example sets 2.0
> 
> 3. for itk::GradientDescentOptimizer:
> - learning rate: set to 20.0 in the example
> - number of iterations: 200
> 
> Can you provide any hints/guidelines of how to play with these parameters for different image modality?  In particular, those work for ultrasound images?  Also is there any text I can find to explain these parameters in more details?
> 
> Thanks a lot,
> Thomas
> 
> 
> 
================================================================

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


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkTranslationTransform.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::MattesMutualInformationImageToImageMetric< ImageType, ImageType >  MetricType;

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



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

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



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

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


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

  metric->SetNumberOfHistogramBins( 20 );
  metric->SetNumberOfSpatialSamples( 10000 );

  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;
}


--------------030304000207020101040806--