[Insight-users] PointSet to PointSet registration problem with LevenbergMarquardtOptimizer and VersorRigid3DTransform

[Neuhaus Jochen]

Hi List,

I need to register two sets of points with unknown point
correspondences. I use the code of example 2 in the software guide in
section 8.17, pages 492ff
(Examples/Patented/IterativeClosestPoint2.cxx), I just exchanged the
Euler3DTransform with a VersorRigid3DTransform.

When using artificial point sets that have just a uniform translation
between them, the optimizer finds the correct transformation.

When using real world point sets (from a tracking device) that are also
rotated and do not fit perfectly, the optimizer fails to find a good
transformation.

For testing, I established point correspondences manually and used the
itk::LandmarkTransformInitializer to initialize a
VersorRigid3DTransform.
The transform parameters were [0.422982, 0.403292, -0.581769, -27.3493,
155.878, 1790.09] (first the 3 versors, then the 3 translation
parameters). This resulted in a FRE (fiducial registration error) of
about 1.5mm.

Feeding the same point sets (with correct point correspondences) into
the optimizer resulted in these transform parameters:
[2.08919e-006, 0.00995692, 0.116826, -19.5648, 176.277, 1790.94]

As you can see, the versors were hardly changed. The
EuclideanDistancePointMetric had the following values: 
[56.2236, 51.0406, 28.7475, 109.227, 109.721, 78.9683] (this lead to a
FRE of about 50mm!)


Why does the optimizer stop when the cost function is that bad? As
proven with the itk::LandmarkTransformInitializer, a much better
solution exists.

I suspect that the parameter scaling could be the cause. Does anyone
have a hint of how to set the scales for a VersorRigid3DTransform?


My registration code is as following:

/* lots of type definitions */
  typedef itk::PointSet<mitk::ScalarType, 3> PointSetType;

  typedef itk::EuclideanDistancePointMetric< PointSetType, PointSetType>
MetricType;
  typedef itk::VersorRigid3DTransform< double > TransformType;
  typedef TransformType ParametersType;
  typedef itk::PointSetToPointSetRegistrationMethod< PointSetType,
PointSetType > RegistrationType;

  MetricType::Pointer metric = MetricType::New();  

  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();

  itk::LevenbergMarquardtOptimizer::Pointer optimizer =
itk::LevenbergMarquardtOptimizer::New();
  optimizer->SetUseCostFunctionGradient(false);

  RegistrationType::Pointer registration = RegistrationType::New();
  // Scale the translation components of the Transform in the Optimizer
  itk::LevenbergMarquardtOptimizer::ScalesType
scales(transform->GetNumberOfParameters());
  const double translationScale = 5000;
  const double rotationScale = 1.0; // dynamic range of rotations
  scales[0] = 1.0 / rotationScale;
  scales[1] = 1.0 / rotationScale;
  scales[2] = 1.0 / rotationScale;
  scales[3] = 1.0 / translationScale;
  scales[4] = 1.0 / translationScale;
  scales[5] = 1.0 / translationScale;
  //scales.Fill(0.01);
  unsigned long numberOfIterations = 80000;
  double gradientTolerance = 1e-10; // convergence criterion
  double valueTolerance = 1e-10; // convergence criterion
  double epsilonFunction = 1e-10; // convergence criterion
  optimizer->SetScales( scales );
  optimizer->SetNumberOfIterations( numberOfIterations );
  optimizer->SetValueTolerance( valueTolerance );
  optimizer->SetGradientTolerance( gradientTolerance );
  optimizer->SetEpsilonFunction( epsilonFunction );

  registration->SetInitialTransformParameters(
transform->GetParameters() );
    registration->SetMetric( metric );
  registration->SetOptimizer( optimizer );
  registration->SetTransform( transform );
  registration->SetFixedPointSet( targetPointSet );
  registration->SetMovingPointSet( sourcePointSet );

  registration->Update();

  std::cout << "ICP successful: Solution = " <<
transform->GetParameters() << std::endl;


Sorry for the long email. 
Thank you for your help!
  Jochen