[Insight-users] pointsettoimageregistration - setup & transform results

[Dean Inglis]

I am registering a set of 3D points taken using a touch probe
over a subject's face to a 3D CT scan of the subject.  My pipeline
compiles and now runs without crashing but the results are
not acceptable: the output points are not positioned remotely
correctly.   I have a few issues that I would could use some
guidance on

1) does the pipeline posted below look reasonable or are there
better components that should be considered (e.g., different metric,
interpolator etc)

2) what parameters need to be tuned and how should they be
initialized?

3) since the points need to be rotated and translated, I am
using  a VersorRigid3DTransform.  Do I need to initialize it
with the center of the fixed points?  When I transform the input
points to get the final output points, is it sufficient to just do
    FixedPointType inputPt  = fixedPointSet->GetPoint( i );
    FixedPointType outputPt = transform->TransformPoint( inputPt );
or do I need an inverse transform and do I also need to take into
account center?

Dean

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

  // set the values associated with the points to be a factor of the stdev + 
mean of the cast output
  float defaultScaleFactor = -1.0;
  defaultValue = static_cast< FixedPointDataType >( defaultScaleFactor * 
statisticsImageFilter->GetSigma() + statisticsImageFilter->GetMean() );

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

  typedef itk::ShrinkImageFilter <UnsignedCharImageType, 
UnsignedCharImageType> ShrinkImageFilterType;
  ShrinkImageFilterType::Pointer shrinkFilter = 
ShrinkImageFilterType::New();
  shrinkFilter->SetInput(castFilter->GetOutput());
  shrinkFilter->SetShrinkFactor(0, 2); // shrink the dimension by a factor 
of 2 (i.e. 100 gets changed to 50)
  shrinkFilter->SetShrinkFactor(1, 2);
  shrinkFilter->SetShrinkFactor(2, 2);
  shrinkFilter->Update();

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

  typedef itk::MeanSquaresPointSetToImageMetric< FixedPointSetType, 
UnsignedCharImageType > MetricType;
  MetricType::Pointer metric = MetricType::New();
  typedef MetricType::TransformType TransformBaseType;
  typedef TransformBaseType::ParametersType ParametersType;

    ParametersType parameters( transform->GetNumberOfParameters() );

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

    parameters[0] = versor.GetX();   // Rotation axis * sin(t/2)
    parameters[1] = versor.GetY();
    parameters[2] = versor.GetZ();
    parameters[3] = 0.0;             // Translation
    parameters[4] = 0.0;
    parameters[5] = 0.0;

    transform->SetParameters( parameters );

  typedef itk:: NearestNeighborInterpolateImageFunction<
                                    UnsignedCharImageType, 
CoordinateRepresentationType     >    InterpolatorType;
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  //interpolator->DebugOn();

  typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
  OptimizerType::Pointer optimizer = OptimizerType::New();

  unsigned long   numberOfIterations =   50;
  double          maximumStepLength  =  1.0;  // no step will be larger than 
this
  double          minimumStepLength  =  0.01; // convergence criterion
  double          gradientTolerance  =  1e-6; // convergence criterion

  OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
  scales.Fill( 1.0 );

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


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

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

  registration->SetMetric(        metric        );
  registration->SetOptimizer(     optimizer     );
  registration->SetTransform(     transform     );
  registration->SetFixedPointSet( fixedPointSet );
  registration->SetMovingImage(   shrinkFilter->GetOutput()   );
  registration->SetInterpolator(  interpolator  );
  registration->SetInitialTransformParameters( transform->GetParameters());