[Insight-users] Help With Affine Transform

Huss hussaintameem at gmail.com
Mon Aug 2 12:57:49 EDT 2004


Hello Everyone,

I am using the MultiResMIRegistration package and I am trying to
implement the Affine Transform and gradient descent optimizer. For
some reason I am not able to implement the scaling for which I intend
to use the Affine Transform for. I would really appreciate if someone
can check and let me know whats wrong and if I made the correct
changes. I am sending you the MIMRegistration.txt and
MIMRegistration.h files and I have indicated the changes i made to
implement Affine.

MIMRegistration.h:
-----------------------------------------------------------------------------------------------

#ifndef _MIMRegistrator_h
#define _MIMRegistrator_h

#include "itkObject.h"
#include "itkMultiResolutionImageRegistrationMethod.h"
#include "itkAffineTransform.h"

// #include "itkQuaternionRigidTransform.h" ---> Removed

#include "itkMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
// #include "itkQuaternionRigidTransformGradientDescentOptimizer.h" -->Removed
// ---> new Added
#include "itkGradientDescentOptimizer.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkRecursiveMultiResolutionPyramidImageFilter.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkNormalizedMutualInformationHistogramImageToImageMetric.h"


#include "CommandIterationUpdate.h"

#include "itkArray.h"

namespace itk
{

template <typename TFixedImage, typename TMovingImage>
class MIMRegistrator : public Object
{
public:

  /** Standard class typedefs. */
  typedef MIMRegistrator Self;
  typedef Object Superclass;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self>  ConstPointer;

  float maxStep; 
  float minStep; 

  /** Run-time type information (and related methods). */
  itkTypeMacro(MIMRegistrator, Object);

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Fixed Image Type. */
  typedef TFixedImage FixedImageType;

  /** Moving Image Type. */
  typedef TMovingImage MovingImageType;

  /** Image dimension enumeration. */
  itkStaticConstMacro (ImageDimension, unsigned int,
TFixedImage::ImageDimension);

  /** Transform Type. */
 //  typedef QuaternionRigidTransform< double >       TransformType;
-----> Changed
  typedef AffineTransform<double, ImageDimension>        TransformType;

  /** Optimizer Type. */
 //  typedef QuaternionRigidTransformGradientDescentOptimizer   
OptimizerType;  ---> Changed
  typedef  GradientDescentOptimizer                 OptimizerType;

  /** Metric Type. */
  typedef MutualInformationImageToImageMetric< 
                                    FixedImageType, 
                                    MovingImageType >    MetricType;
  unsigned int nBins;
  

  /** Interpolation Type. */
  typedef LinearInterpolateImageFunction< 
                                    MovingImageType,
                                    double          >    InterpolatorType;

  /** Fixed Image Pyramid Type. */
  typedef RecursiveMultiResolutionPyramidImageFilter<
                                    FixedImageType,
                                    FixedImageType  >    FixedImagePyramidType;

  /** Moving Image Pyramid Type. */
  typedef RecursiveMultiResolutionPyramidImageFilter<
                                    MovingImageType,
                                    MovingImageType  >   MovingImagePyramidType;

  /** Registration Method. */
  typedef MultiResolutionImageRegistrationMethod< 
                                    FixedImageType, 
                                    MovingImageType >    RegistrationType;

  /** Transform parameters type. */
  typedef typename RegistrationType::ParametersType     ParametersType;

  /** DoubleArray type. */
  typedef Array<double>  DoubleArray;

  /** UnsignedIntArray type. */
  typedef Array<unsigned int> UnsignedIntArray;

  /** ShrinkFactorsArray type. */
  typedef FixedArray<unsigned
int,itkGetStaticConstMacro(ImageDimension)> ShrinkFactorsArray;

  /** Affine transform type. */
  typedef AffineTransform<double,itkGetStaticConstMacro(ImageDimension)>
  AffineTransformType;
  typedef typename AffineTransformType::Pointer AffineTransformPointer;

  /** Set the fixed image. */
  itkSetObjectMacro( FixedImage, FixedImageType );

  /** Set the moving image. */
  itkSetObjectMacro( MovingImage, MovingImageType );

  /** Set the number of resolution levels. */
  itkSetClampMacro( NumberOfLevels, unsigned short, 1,
    NumericTraits<unsigned short>::max() );

  /** Set the translation parameter scales. */
  itkSetClampMacro( TranslationScale, double, 0.0,
    NumericTraits<double>::max() );

  /** Set the image parzen window widths. */
  itkSetClampMacro( MovingImageStandardDeviation, double, 0.0,
    NumericTraits<double>::max() );
  itkSetClampMacro( FixedImageStandardDeviation, double, 0.0,
    NumericTraits<double>::max() );

  /** Set the number of spatial samples. */
  itkSetClampMacro( NumberOfSpatialSamples, unsigned short, 1,
    NumericTraits<unsigned short>::max() );

  /** Set the number of iterations per level. */
  itkSetMacro( NumberOfIterations, UnsignedIntArray );

  /** Set the learning rate per level. */
  itkSetMacro( LearningRates, DoubleArray );

  /** Set the initial transform parameters. */
  itkSetMacro( InitialParameters, ParametersType );

  /** Set the fixed and moving image shrink factors. */
  itkSetMacro( FixedImageShrinkFactors, ShrinkFactorsArray );
  itkSetMacro( MovingImageShrinkFactors, ShrinkFactorsArray );

  /** Method to execute the registration. */
  virtual void Execute();

  /** Get number of parameters. */
  unsigned long GetNumberOfParameters()
    { return m_Transform->GetNumberOfParameters(); }

  /** Get computed transform parameters. */
  const ParametersType& GetTransformParameters();

  /** Get computed affine transform. */
  AffineTransformPointer GetAffineTransform();

  /** Initialize registration at the start of new level. */
  void StartNewLevel();

protected:
  MIMRegistrator();
  ~MIMRegistrator();

private:
  MIMRegistrator( const Self& ); //purposely not implemented
  void operator=( const Self& ); //purposely not implemented

  typename FixedImageType::Pointer            m_FixedImage;
  typename MovingImageType::Pointer           m_MovingImage;
  typename TransformType::Pointer             m_Transform;
  typename OptimizerType::Pointer             m_Optimizer;
  typename MetricType::Pointer                m_Metric;
  typename InterpolatorType::Pointer          m_Interpolator;
  typename FixedImagePyramidType::Pointer     m_FixedImagePyramid;
  typename MovingImagePyramidType::Pointer    m_MovingImagePyramid;
  typename RegistrationType::Pointer          m_Registration;
  CommandIterationUpdate::Pointer             m_Observer;

  unsigned short                              m_NumberOfLevels;
  double                                      m_TranslationScale;
  double                                      m_MovingImageStandardDeviation;
  double                                      m_FixedImageStandardDeviation;
  unsigned short                              m_NumberOfSpatialSamples;

  UnsignedIntArray                            m_NumberOfIterations;
  DoubleArray                                 m_LearningRates;

  ShrinkFactorsArray                          m_MovingImageShrinkFactors;
  ShrinkFactorsArray                          m_FixedImageShrinkFactors;

  ParametersType                              m_InitialParameters;
  AffineTransformPointer                      m_AffineTransform;

  unsigned long                               m_Tag;


};

} // namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "MIMRegistrator.txx"
#endif

#endif

---------------------------------------
MIMRegistration.txx


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

  //-----> added
  minStep = 0.5; 
  maxStep = 25;  

  // 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_Transform->SetIdentity();

  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);
  //  m_NumberOfIterations.Fill( 10 ); -------> changed 10 to 1000 below
  m_NumberOfIterations.Fill( 1000 );
  m_LearningRates = DoubleArray(1);
  m_LearningRates.Fill( 1e-4 );

  // m_InitialParameters = ParametersType(
m_Transform->GetNumberOfParameters() );----->changed
  // m_InitialParameters.Fill( 0.0 ); ----->changed
  // m_InitialParameters[3] = 1.0;  ---->changed

  //----> Replaced Above lines with the one below 
  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 = 4; j < 7; j++ )---> Changed as below
  for ( int j = 9; j < 12; j++ )
    {
    scales[j] = m_TranslationScale;
    }

  m_Optimizer->SetScales( scales );
  m_Optimizer->MaximizeOn();

  // Setup the metric
  m_Metric->SetMovingImageStandardDeviation( m_MovingImageStandardDeviation );
  m_Metric->SetFixedImageStandardDeviation( m_FixedImageStandardDeviation );
  m_Metric->SetNumberOfSpatialSamples( m_NumberOfSpatialSamples );

  // 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_InitialParameters ); ------> changed
  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->GetRotationMatrix() );
------> Changed
  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->SetNumberOfIterations( m_NumberOfIterations[level] );
    }

  if ( m_LearningRates.Size() >= level + 1 )
    {
    m_Optimizer->SetLearningRate( m_LearningRates[level] );
    }

  std::cout << " No. Iterations: " 
            << m_Optimizer->GetNumberOfIterations()
            << " Learning rate: "
            << m_Optimizer->GetLearningRate()
            << std::endl;

}


} // namespace itk


#endif


More information about the Insight-users mailing list