[Insight-users] pointset to image registration exception

[Dean Inglis]

my current pipeline is hopefully set up to register a set of 3D surface 
points
obtained from a motion tracking/probing system to a 3D CT scan.
I am using an inverted gradient magnitude of the CT scan as the moving
image and setting the values associated with the points to be a constant:
the mean of the inverted gradient magnitude image + a factor * the std dev.
I need the transformation to allow for rotation, translation and (assuming
isotropic) scaling so I am using the Similarity3DTransform class.
With the code below I am getting an exception error that I dont know how
to resolve:

itk::ExceptionObject (00C7FACC)
Location: "class itk::CovariantVector<double,3> *__thiscall 
itk::ImportImageCont
ainer<unsigned long,class itk::CovariantVector<double,3> 
 >::AllocateElements(uns
igned long) const"
File: 
d:\developer\vs10\x86\static\install\include\insighttoolkit-4.0\itkImportI
mageContainer.txx
Line: 191
Description: Failed to allocate memory for image.

ITK is built statically, visual studio 2010, nmake, release, cmake 2.8.4 on
Win 7.

Dean

#include "itkPointSetToImageRegistrationMethod.h"
#include "itkPointSet.h"
#include <iostream>
#include <fstream>
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkIntensityWindowingImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkMinimumMaximumImageCalculator.h"
#include "itkStatisticsImageFilter.h"
#include "itkSimilarity3DTransform.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkPointSetToImageRegistrationMethod.h"
#include "itkVersorRigid3DTransformOptimizer.h"
#include "itkNormalizedCorrelationPointSetToImageMetric.h"
//#include "itkCommandIterationUpdate.h"


//
// Observer to the optimizer
//
class CommandIterationUpdate : public itk::Command
{
public:
  typedef  CommandIterationUpdate   Self;
  typedef  itk::Command             Superclass;
  typedef  itk::SmartPointer<Self>  Pointer;
  itkNewMacro( Self );
protected:
  CommandIterationUpdate() {};
public:
  typedef   itk::RegularStepGradientDescentOptimizer  OptimizerType;
  typedef   const OptimizerType *                     OptimizerPointer;

  void Execute(itk::Object *caller, const itk::EventObject & event)
    {
    Execute( (const itk::Object *)caller, event);
    }

  void Execute(const itk::Object * object, const itk::EventObject & event)
    {
    OptimizerPointer optimizer =
      dynamic_cast< OptimizerPointer >( object );
    if( typeid( event ) != typeid( itk::IterationEvent ) )
      {
      return;
      }

    OptimizerType::DerivativeType gradient = optimizer->GetGradient();
    OptimizerType::ScalesType     scales   = optimizer->GetScales();

    double magnitude2 = 0.0;

    for(unsigned int i=0; i<gradient.size(); i++)
      {
      const double fc = gradient[i] / scales[i];
      magnitude2 += fc * fc;
      }

    const double gradientMagnitude = vcl_sqrt( magnitude2 );

    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << gradientMagnitude << "   ";
    std::cout << optimizer->GetCurrentPosition() << std::endl;
    }
};

int main( int argc, char * argv[] )
{
  // define the object that will hold the 3D probe points
  const unsigned int Dimension = 3;
  typedef unsigned char FixedPointDataType;

  // associate a default value with each point
  // we will use an inverted gradient magnitude image to register the
  // points to, so it should be a low salar value
  // for example, if the gradient magnitude image (GMI) background is zero,
  // then the inverted GMI will have a high value
  FixedPointDataType defaultValue = 0;

  typedef itk::PointSet< FixedPointDataType, Dimension > FixedPointSetType;
  FixedPointSetType::Pointer fixedPointSet = FixedPointSetType::New();

  typedef FixedPointSetType::PointsContainer  FixedPointsContainer;
  FixedPointsContainer::Pointer fixedPointContainer  = 
FixedPointsContainer::New();

  typedef FixedPointSetType::PointType FixedPointType;
  FixedPointType fixedPoint;

  // Read the file containing coordinates of fixed points.
  std::ifstream   fixedFile;

  fixedFile.open( argv[1] );
  if( fixedFile.fail() )
  {
    std::cerr << "Error opening points file with name : " << std::endl;
    std::cerr << argv[1] << std::endl;
    return EXIT_FAILURE;
  }

  unsigned int pointId = 0;
  fixedFile >> fixedPoint;
  while( !fixedFile.eof() )
  {
    fixedPointContainer->InsertElement( pointId, fixedPoint );
    fixedFile >> fixedPoint;
    std::cout << "point " << pointId << ": " << fixedPoint[0] << ", " << 
fixedPoint[1] << ", " << fixedPoint[2] << std::endl;
    pointId++;
  }
  fixedPointSet->SetPoints( fixedPointContainer );
  std::cout << "number of fixed Points = " << 
fixedPointSet->GetNumberOfPoints() << std::endl;

  // read in the original dicom image
  std::string inputFileName = argv[2];
  // Setup image types
  typedef itk::Image< float,  3 >         FloatImageType;
  typedef itk::Image< unsigned short, 3 > UnsignedShortImageType;
  typedef itk::Image< unsigned char, 3 >  UnsignedCharImageType;

  // image reader
  typedef itk::ImageFileReader< UnsignedShortImageType > ImageReaderType;
  ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName( inputFileName );
  reader->Update();

  std::cout << "read image " << inputFileName << std::endl;

  // gradient magnitude
  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< 
UnsignedShortImageType, FloatImageType >  GradientFilterType;
  GradientFilterType::Pointer gradientFilter = GradientFilterType::New();
  gradientFilter->SetInput( reader->GetOutput() );
  gradientFilter->SetSigma( 1.0 );
  gradientFilter->Update();

  std::cout << "grad mag calculated " << std::endl;

  typedef itk::MinimumMaximumImageCalculator< FloatImageType > 
CalculatorType;
  CalculatorType::Pointer calculator = CalculatorType::New();
  calculator->SetImage( gradientFilter->GetOutput() );
  calculator->ComputeMaximum();
  float maximum = calculator->GetMaximum();

  std::cout << "max grad mag calculated " << maximum << std::endl;

  // rescale the intensities
  typedef itk::IntensityWindowingImageFilter< FloatImageType, FloatImageType 
 > IntensityFilterType;
  IntensityFilterType::Pointer intensityFilter = IntensityFilterType::New();
  intensityFilter->SetInput( gradientFilter->GetOutput() );
  intensityFilter->SetWindowMinimum( 0 );
  intensityFilter->SetWindowMaximum( maximum );
  intensityFilter->SetOutputMinimum( 255 );
  intensityFilter->SetOutputMaximum( 0 );
  intensityFilter->Update();

  std::cout << "intensity inversion and rescaling calculated " << std::endl;

  // cast to unsigned  char
  typedef itk::CastImageFilter< FloatImageType, UnsignedCharImageType > 
CastFilterType;
  CastFilterType::Pointer castFilter = CastFilterType::New();
  castFilter->SetInput( intensityFilter->GetOutput() );
  castFilter->Update();

  std::cout << "cast to uchar calculated " << std::endl;

  // write the image out
  typedef itk::ImageFileWriter<  UnsignedCharImageType > ImageWriterType;
  ImageWriterType::Pointer writer = ImageWriterType::New();
  writer->SetFileName( "testout2.mhd" );
  writer->SetInput( castFilter->GetOutput() );
  writer->Update();

  typedef itk::StatisticsImageFilter< UnsignedCharImageType > 
StatisticsImageFilterType;
  StatisticsImageFilterType::Pointer statisticsImageFilter
          = StatisticsImageFilterType::New ();
  statisticsImageFilter->SetInput( castFilter->GetOutput() );
  statisticsImageFilter->Update();

  std::cout << "Mean: " << statisticsImageFilter->GetMean() << std::endl;
  std::cout << "Std.: " << statisticsImageFilter->GetSigma() << std::endl;

  // set the values associated with the points to be a factor of the stdev + 
mean of the cast output
  float defaultScaleFactor = 0.0;
  defaultValue = static_cast< FixedPointDataType >( defaultScaleFactor * 
statisticsImageFilter->GetSigma() + statisticsImageFilter->GetMean() );
  std::cout << "default fixed point value: " << static_cast<int>( 
defaultValue )
            << " scale: " << defaultScaleFactor
            << " cast mean: " << statisticsImageFilter->GetMean()
            << " cast stddev: " << statisticsImageFilter->GetSigma()
            << std::endl;

  for(int i = 0 ; i < fixedPointSet->GetNumberOfPoints(); ++i )
  {
    fixedPointSet->SetPointData( i, defaultValue );
  }


  //-----------------------------------------------------------
  // Set up a Transform - does rotation (using versors), translation and 
isotropic scaling
  // for anisotropic scaling use itkScaleVersor3DTransform
  //-----------------------------------------------------------
  typedef double CoordinateRepresentationType;
  typedef itk::Similarity3DTransform< CoordinateRepresentationType   > 
TransformType;
  TransformType::Pointer transform = TransformType::New();

  // Start from an Identity transform (in a normal case, the user
  // can probably provide a better guess than the identity...
  transform->SetIdentity();



  //-----------------------------------------------------------
  // Set up Metric (also known as a cost function)
  // must inherit from itkPointSetToImageMetric
  // other options are:
  //  itkMeanReciprocalSquareDifferencePointSetToImageMetric
  //  itkMeanSquaresPointSetToImageMetric
  //  itkNormalizedCorrelationPointSetToImageMetric
  //-----------------------------------------------------------
  typedef itk::NormalizedCorrelationPointSetToImageMetric< 
FixedPointSetType, UnsignedCharImageType > MetricType;
  MetricType::Pointer metric = MetricType::New();
  typedef MetricType::TransformType TransformBaseType;
  typedef TransformBaseType::ParametersType ParametersType;

  //------------------------------------------------------------
  // Set up transform parameters
  //------------------------------------------------------------
  ParametersType parameters( transform->GetNumberOfParameters() );
  // initialize the offset/vector part
  for( unsigned int k = 0; k <Dimension; k++ )
  {
  parameters[k]= 0.0f;
  }

   //------------------------------------------------------------
  // Set up an Interpolator
  //------------------------------------------------------------
  typedef itk::LinearInterpolateImageFunction< UnsignedCharImageType, 
CoordinateRepresentationType > InterpolatorType;
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  //------------------------------------------------------------

  //------------------------------------------------------------
  // Optimizer Type
  //------------------------------------------------------------
  typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
  OptimizerType::Pointer optimizer = OptimizerType::New();

  unsigned long   numberOfIterations =   50;
  //double          translationScale   =  1.0; // not used
  double          maximumStepLenght  =  1.0;  // no step will be larger than 
this
  double          minimumStepLenght  =  0.01; // convergence criterion
  double          gradientTolerance  =  1e-6; // convergence criterion

  // Scale the translation components of the Transform in the Optimizer
  //OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
  //scales.Fill( 1.0 );

  //optimizer->SetScales( scales );
  optimizer->SetNumberOfIterations( numberOfIterations );
  optimizer->SetMinimumStepLength( minimumStepLenght );
  optimizer->SetMaximumStepLength( maximumStepLenght );
  optimizer->SetGradientMagnitudeTolerance( gradientTolerance );
  optimizer->MinimizeOn();

  typedef CommandIterationUpdate    IterationObserverType;
  IterationObserverType::Pointer        iterationObserver;
  //optimizer->AddObserver( itk::IterationEvent(), iterationObserver );

  //------------------------------------------------------------
  // registration method
  //------------------------------------------------------------
  typedef itk::PointSetToImageRegistrationMethod< FixedPointSetType, 
UnsignedCharImageType > RegistrationType;
  RegistrationType::Pointer registration = RegistrationType::New();

  //------------------------------------------------------
  // Connect all the components required for Registration
  //------------------------------------------------------
  registration->SetMetric(        metric        );
  registration->SetOptimizer(     optimizer     );
  registration->SetTransform(     transform     );
  registration->SetFixedPointSet( fixedPointSet );
  registration->SetMovingImage(   castFilter->GetOutput()   );
  registration->SetInterpolator(  interpolator  );

  try
    {
    registration->StartRegistration();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cout << e << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "output written " << std::endl;
  return EXIT_SUCCESS;
}