[Insight-users] pointset to image registration tuning

[Dean Inglis]

I am registering a set of (19) points obtained from a touch probe on
a subject's face.  I have a 3D CT scan of the subject's head that I have
processed into a cost image by the following sequence:
1) ImageFileReader - read in an unsigned short 3D CT image (background has 
value 0)
2) ShrinkImageFilter - shrink all dimensions by a factor of 2
3) ImageRegionConstIteratorWithIndex - find the first pixel index that has a 
value of 0
4) ConfidenceConnectedImageFilter - threshold to unsigned char using the 
seed from 3)
5) SliceBySliceImageFilter / BinaryFillholeImageFilter - fill all holes 
slice by slice with a value of 0
- now there is an image of the background (255) with the entire region 
within the surface of the
subject set to 0, sinuses and ear canals etc removed
6) GradientMagnitudeRecursiveGaussianImageFilter -  generates an image of 
the subject's surface
7) MinimumMaximumImageCalculator  - find the maximum gradient magnitude 
corresponding to strong edges
8) IntensityWindowingImageFilter - rescale the gradient magnitude image to 
0 - 255
9) CastImageFilter - cast to unsigned char

I then have a pointset to image registration pipeline defined in the code 
below.

I have tried setting the transform center to the fixed points center
and an initial translation to the difference between the fixed points and 
center of gravity
of the image.  I have passed in the points both in the original perturbed 
positions (non-registered)
and in a manually registered (using paraview transformation filter) 
configuration.
The pipeline should give a near zero output in the latter case, but does 
not! I seem
to have hit a wall on tuning the optimizer and transform, as no matter what 
I have tried,
(vary step lengths, number of iterations etc) I never get a satisfactory 
(not even close!) registration.
Any insight would be greatly appreciated.

Dean

///////////////////////////////////

  const unsigned int Dimension = 3;
  typedef unsigned char FixedPointDataType;
  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;
  double cx = 0.;
  double cy = 0.;
  double cz = 0.;
  while( !fixedFile.eof() )
  {
    fixedFile >> fixedPoint;
    fixedPointContainer->InsertElement( pointId, fixedPoint );
    std::cout << "point " << pointId << ": " << fixedPoint[0] << ", " << 
fixedPoint[1] << ", " << fixedPoint[2] << std::endl;
    cx += fixedPoint[0];
    cy += fixedPoint[1];
    cz += fixedPoint[2];
    pointId++;
  }
  fixedPointSet->SetPoints( fixedPointContainer );
  int numFixedPts = fixedPointSet->GetNumberOfPoints();
  cx /= numFixedPts;
  cy /= numFixedPts;
  cz /= numFixedPts;
  std::cout << "number of fixed Points = " << numFixedPts << std::endl;

  std::string inputFileName = argv[2];
  // Setup image types
  typedef itk::Image< float,  3 >         FloatImageType;
  typedef itk::Image< unsigned char, 3 >  UnsignedCharImageType;

  // image reader
  // read in the cost image
  //
  typedef itk::ImageFileReader< UnsignedCharImageType > ImageReaderType;
  ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName( inputFileName );
  reader->ReleaseDataFlagOff();
  reader->Update();

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

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

  typedef itk::Vector<double,3>  VectorType;
  typedef itk::ImageMomentsCalculator<UnsignedCharImageType> 
MomentCalculatorType;

    /* Compute the moments */
    MomentCalculatorType::Pointer moments = MomentCalculatorType::New();
    moments->SetImage( reader->GetOutput() );
    moments->Compute();
    moments->Print( std::cout );

    VectorType ccg = moments->GetCenterOfGravity();

    std::cout << "image cofg: " << ccg[0] << ", " << ccg[1] << ", " << 
ccg[2] << std::endl;

  //-----------------------------------------------------------
  // Set up a Transform - does rotation (using versors), adn translation, 
scaling not required
  //-----------------------------------------------------------
  typedef double CoordinateRepresentationType;
  typedef itk::VersorRigid3DTransform< CoordinateRepresentationType   > 
TransformType;
  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();

  TransformType::InputPointType centerFixed;
  centerFixed[0] = cx;
  centerFixed[1] = cy;
  centerFixed[2] = cz;
  transform->SetCenter(centerFixed);

  TransformType::OutputVectorType translation;
  translation[0] = ccg[0] - cx;
  translation[1] = ccg[1] - cy;
  translation[2] = ccg[2] - cz;
  transform->SetTranslation( translation );

  //-----------------------------------------------------------
  // Set up Metric (also known as a cost function)
  // must inherit from itkPointSetToImageMetric
  // other options are:
  //  itkMeanReciprocalSquareDifferencePointSetToImageMetric
  //  itkMeanSquaresPointSetToImageMetric
  //  itkNormalizedCorrelationPointSetToImageMetric
  //-----------------------------------------------------------

  typedef itk::MeanReciprocalSquareDifferencePointSetToImageMetric< 
FixedPointSetType, UnsignedCharImageType > MetricType;
  MetricType::Pointer metric = MetricType::New();
  metric->SetLambda( 10 );

  typedef MetricType::TransformType TransformBaseType;
  typedef TransformBaseType::ParametersType ParametersType;

  //------------------------------------------------------------
  // Set up transform parameters - not used
  //------------------------------------------------------------
  /*
    ParametersType parameters( transform->GetNumberOfParameters() );

  //  Versor type
    typedef    TransformType::VersorType      VersorType;
    VersorType versor;

    parameters[0] = versor.GetX();
    parameters[1] = versor.GetY();
    parameters[2] = versor.GetZ();
    parameters[3] = 0.0;             // Translation
    parameters[4] = 0.0;
    parameters[5] = 0.0;
    parameters[6] = 1.0;             // Scale

    transform->SetParameters( parameters );
  */
   //------------------------------------------------------------
  // Set up an Interpolator
  //------------------------------------------------------------
  typedef itk:: NearestNeighborInterpolateImageFunction<
                                    UnsignedCharImageType, 
CoordinateRepresentationType     >    InterpolatorType;
  InterpolatorType::Pointer interpolator = InterpolatorType::New();

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

  unsigned long   numberOfIterations =   500;
  double          maximumStepLength  =  1.0;  // no step will be larger than 
this
  double          minimumStepLength  =  0.001; // 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 );

  // set up the translation scales: translations are much larger than 
rotations
  scales[3] = 1./100.;
  scales[4] = 1./100.;
  scales[5] = 1./100.;

  optimizer->SetScales( scales );
  optimizer->SetNumberOfIterations( numberOfIterations );
  optimizer->SetMinimumStepLength( minimumStepLength );
  optimizer->SetMaximumStepLength( maximumStepLength );
  optimizer->SetGradientMagnitudeTolerance( gradientTolerance );
  optimizer->MinimizeOff();
  optimizer->MaximizeOn();

  typedef IterationCallback< OptimizerType >   IterationCallbackType;
  IterationCallbackType::Pointer callback = IterationCallbackType::New();
  callback->SetOptimizer( optimizer );

  //------------------------------------------------------------
  // 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(  reader->GetOutput()  );
  registration->SetInterpolator(  interpolator  );
  registration->SetInitialTransformParameters( transform->GetParameters());

  std::cout << "starting registration .. " << std::endl;

  try
    {
    registration->StartRegistration();