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

[Bing Jian]

Hi Jochen,
   If you need to register two point sets with unknown point
correspondences,  you may want to take a look at
http://code.google.com/p/gmmreg/
and the links therein.

Thanks,
Bing

On Wed, Sep 2, 2009 at 7:21 AM, Neuhaus
Jochen<j.neuhaus at dkfz-heidelberg.de> wrote:
> 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
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the ITK FAQ at: http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>