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

[Neuhaus Jochen]

Hi Bing,
Thanks for the hint. The gmmreg project offers deformable registration algorithms, but I'm looking specifically for rigid registration of two pointsets. 

I still did not manage to optimize the rotation part of the versor transform correctly :-( I could post my example pointsets, if that would help.

Thanks
 Jochen

> -----Ursprüngliche Nachricht-----
> Von: Bing Jian [mailto:bing.jian at gmail.com]
> Gesendet: Donnerstag, 3. September 2009 23:57
> An: Neuhaus Jochen
> Cc: insight-users at itk.org
> Betreff: Re: [Insight-users] PointSet to PointSet registration problem
> with LevenbergMarquardtOptimizer and VersorRigid3DTransform
> 
> 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
> >