[Insight-users] FEMRegistrationFilter: multi-resolution techniques

[Luis Ibanez]

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



Hi Jens,


You have three easy options for creating the
subsampled volumes.


A) Use the attached program
    which as also been CVS commited under

     Insight/Examples/Filtering/SubsampleVolume.cxx


B) Use the MultiResolutionPyramidImageFilter

http://www.itk.org/Insight/Doxygen/html/classitk_1_1MultiResolutionPyramidImageFilter.html

    Select to have a single level and set the scale factors
    to the level you want.


C) Use the RecursiveMultiResolutionPyramidImageFilter

http://www.itk.org/Insight/Doxygen/html/classitk_1_1RecursiveMultiResolutionPyramidImageFilter.html

    Select to have a single level and set the scale factors
    to the level you want.



In all the cases, the easy way to go is to save the
resulting subsampled volumes into new files and then
play with the parameters of the FEM registration method
using these reduced files.  In these context you use
the FEM framework at single resolution until you find
good parameters for solving these first level.


About your subsampling schedule:

I would suggest you to not subsample Z at all.
Your original image only have 64 slices and you
may need the information in that direction.

You may get better results by using

  =========================
   level	resolution
  -------------------------
    0	128 x  96 x 64
    1	256 x 192 x 64
    2	512 x 384 x 64
  =========================



Please let us know if you have further questions.


Thanks



    Luis




-----------------------
Jens Fisseler wrote:

> Hi everybody!
> 
> I'm quite new to ITK and I'm currently using the FEMRegistrationFilter 
> to register two MRI liver volume datasets. The results this far are
> quite promising, but I've a question regarding the use of the
> multi-resolution technique.
> 
> My datasets have the resolution 512x384x64 and I'm using a three level
> pyramid:
> 
> ===================
>  level	resolution
> -------------------
>  0	128x96x16
>  1	256x192x32
>  2	512x384x64
> ===================
> 
> As Luis Ibanez wrote in
> http://www.itk.org/pipermail/insight-users/2003-November/005704.html :
> "If a level has not converged, there are few chances that the next level
> will compensate.", I'm trying to optimise the parameters for each level
> before going up to the next. But this is where my problems start. I
> don't quite know how to get access to the deformed volume at a lower
> level in order to assess the registration. Let me explain: If I want to
> assess the deformation at level 0, I would make the registration filter
> stop after this level by setting the second line in my parameter file to
> "1". When I write the deformed image to a file it has the expected
> dimensions but is all black, probably because a small portion of the
> original volume gets deformed, which is black in this area because the
> liver has been segmented.
> 
> Has anybody an idea how I can get access to downsampled deformed volume?
> Is my approach a good idea anyway?
> 
> This is a part of the parameter file I'm using:
> 
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> %
> % Parameters for the single- or multi-resolution techniques
> %
> 
> 3        % Number of levels in the multi-res pyramid (1 = single-res)
> 1        % Highest level to use in the pyramid
> 4 4 4    % Scaling at lowest level of pyramid
> 4 4 4    % Number of pixels per element
> 1.e5 1.e4 1.e4 % Elasticity (E)
> 1.e4 1.e4 1.e4 % Density x capacity (RhoC)
> 1 1 1    % Image energy scaling (gamma) - sets gradient step size
> 2 2 2    % NumberOfIntegrationPoints
> 4 4 4    % WidthOfMetricRegion
> 10 1 1   % MaximumIterations
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> %
> % Parameters for the registration
> %
> 
> 0 0.99 % Similarity metric (0=mean sq, 1 = ncc, 2=pattern int, 3=MI,
> 5=demons)
> 1.0    % Alpha
> 0      % DescentDirection (1 = max, 0 = min)
> 0      % DoLineSearch (0=never, 1=always, 2=if needed)
> 1.e1   % TimeStep
> 0.5    % Landmark variance
> 0      % Employ regridding / enforce diffeomorphism ( >= 1 -> true)
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> %
> % Information about the image inputs
> %
> 
> 512               % Nx (image x dimension)
> 384               % Ny (image y dimension)
> 64                % Nz (image z dimension - not used if 2D)
> ./liver_1.mhd     % ReferenceFileName
> ./liver_2.mhd     % TargetFileName
> 
> 
> Regards,
> 
> 	Jens
> 


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

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: SubsampleVolume.cxx,v $
  Language:  C++
  Date:      $Date: 2004/01/14 16:12:15 $
  Version:   $Revision: 1.1 $

  Copyright (c) Insight Software 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.

=========================================================================*/

//  Software Guide : BeginLatex
//
//  This example illustrates how to perform subsampling of a volume using ITK
//  classes.  In order to avoid aliasing artifacts, the volume must be
//  processed by a low-pass filter before resampling.  Here we use the
//  \doxygen{RecursiveGaussianImageFilter} as low-pass filter.
//
//  Software Guide : EndLatex 


// Software Guide : BeginCodeSnippet
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRecursiveGaussianImageFilter.h"
#include "itkCastImageFilter.h"


int main( int argc, char * argv[] )
{
  if( argc < 6 )
    {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << "  inputImageFile  outputImageFile factorX factorY factorZ" << std::endl; 
    return 1;
    }

  const     unsigned int    Dimension = 3;

  typedef   unsigned char   InputPixelType;

  typedef   float           InternalPixelType;
  typedef   unsigned char   OutputPixelType;

  typedef itk::Image< InputPixelType,    Dimension >   InputImageType;
  typedef itk::Image< InternalPixelType, Dimension >   InternalImageType;
  typedef itk::Image< OutputPixelType,   Dimension >   OutputImageType;

  typedef itk::ImageFileReader< InputImageType  >  ReaderType;
  typedef itk::ImageFileWriter< OutputImageType >  WriterType;

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();

  reader->SetFileName( argv[1] );
  writer->SetFileName( argv[2] );

  const double factorX = atof( argv[3] );
  const double factorY = atof( argv[4] );
  const double factorZ = atof( argv[5] );

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


  InputImageType::ConstPointer inputImage = reader->GetOutput();


  typedef itk::CastImageFilter< InputImageType,
                                InternalImageType >   CastFilterType;

  CastFilterType::Pointer caster = CastFilterType::New();
  caster->SetInput( inputImage );

  typedef itk::RecursiveGaussianImageFilter< 
                                  InternalImageType,
                                  InternalImageType > GaussianFilterType;

  GaussianFilterType::Pointer smootherX = GaussianFilterType::New();
  GaussianFilterType::Pointer smootherY = GaussianFilterType::New();
  GaussianFilterType::Pointer smootherZ = GaussianFilterType::New();

  smootherX->SetInput( caster->GetOutput() );
  smootherY->SetInput( smootherX->GetOutput() );
  smootherZ->SetInput( smootherY->GetOutput() );

  const InputImageType::SpacingType& inputSpacing = inputImage->GetSpacing();

  const double sigmaX = inputSpacing[0] * factorX;
  const double sigmaY = inputSpacing[1] * factorY;
  const double sigmaZ = inputSpacing[2] * factorZ;

  smootherX->SetSigma( sigmaX );
  smootherY->SetSigma( sigmaY );
  smootherZ->SetSigma( sigmaZ );

  smootherX->SetDirection( 0 );
  smootherY->SetDirection( 1 );
  smootherZ->SetDirection( 2 );

  smootherX->SetNormalizeAcrossScale( false );
  smootherY->SetNormalizeAcrossScale( false );
  smootherZ->SetNormalizeAcrossScale( false );

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

  std::cout << "Image Smoothed" << std::endl;

  InternalImageType::ConstPointer smoothedImage = smootherZ->GetOutput();

  typedef itk::ResampleImageFilter<
                  InternalImageType, OutputImageType >  ResampleFilterType;

  ResampleFilterType::Pointer resampler = ResampleFilterType::New();

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

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

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

  resampler->SetInterpolator( interpolator );

  resampler->SetDefaultPixelValue( 0 ); // value for regions without source


  OutputImageType::SpacingType spacing;

  spacing[0] = inputSpacing[0] * factorX;
  spacing[1] = inputSpacing[1] * factorY;
  spacing[2] = inputSpacing[2] * factorZ;

  resampler->SetOutputSpacing( spacing );

  // Use the same origin
  resampler->SetOutputOrigin( inputImage->GetOrigin() );


  InputImageType::SizeType   inputSize = inputImage->GetLargestPossibleRegion().GetSize();
  typedef InputImageType::SizeType::SizeValueType SizeValueType;
  InputImageType::SizeType   size;

  size[0] = static_cast< SizeValueType >( inputSize[0] / factorX );
  size[1] = static_cast< SizeValueType >( inputSize[1] / factorY );
  size[2] = static_cast< SizeValueType >( inputSize[2] / factorZ );

  resampler->SetSize( size );

  resampler->SetInput( smoothedImage );

  writer->SetInput( resampler->GetOutput() );

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

  transform->SetIdentity();

  resampler->SetTransform( transform );

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

  std::cout << "Resampling Done !" << std::endl;

// Software Guide : EndCodeSnippet

  return 0;
}


--------------010308090605040005090309--