<br>Hi Neuhaus,<br><br>In order to make a good selection of scaling parameters <br>for a Versor3DTransform, you should answer the following <br>two questions:<br><br><br>A) How large is the translation that you expect, <br>
will be needed for registering the points ?<br> <br> 1mm ? 10mm ? 2m ?<br><br><br>B) How large of a rotation do you anticipate to be<br> necessary for registering the two datasets ?<br><br> 1 degree ? 10 degrees ? 90 degrees ?<br>
<br> BTW: Note that you will have to convert this <br> to radians.<br><br><br>Once you answer these two question, what you<br>want to do is to setup the rotation scales to 1.0<br>and the translation scales to:<br><br>
ExpectedRotation / ExpectedTranslation<br> <br><br>Note also that, if you need a rotation larger than<br>(about) 30 degrees, it is unlikely that a registration<br>will be able to solve it. You will need to initialize<br>
the Versor with that rotation.<br><br><br>Finally, <br>you may find useful to look at the following<br>Insight Journal paper:<br><br><a href="http://www.insight-journal.org/browse/publication/304">http://www.insight-journal.org/browse/publication/304</a><br>
<a href="http://hdl.handle.net/1926/1497">http://hdl.handle.net/1926/1497</a><br>"A Novel Information-Theoretic Point-Set Measure Based on the Jensen-Havrda-Charvat-Tsallis Divergence"<br>by Tustison N., Awate S., Gee J.<br>
<br><br> <br> Regards,<br><br><br> Luis<br><br><br><br>--------------------------------------------------------------------------------------<br><div class="gmail_quote">On Wed, Sep 2, 2009 at 7:21 AM, Neuhaus Jochen <span dir="ltr"><<a href="mailto:j.neuhaus@dkfz-heidelberg.de">j.neuhaus@dkfz-heidelberg.de</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Hi List,<br>
<br>
I need to register two sets of points with unknown point<br>
correspondences. I use the code of example 2 in the software guide in<br>
section 8.17, pages 492ff<br>
(Examples/Patented/IterativeClosestPoint2.cxx), I just exchanged the<br>
Euler3DTransform with a VersorRigid3DTransform.<br>
<br>
When using artificial point sets that have just a uniform translation<br>
between them, the optimizer finds the correct transformation.<br>
<br>
When using real world point sets (from a tracking device) that are also<br>
rotated and do not fit perfectly, the optimizer fails to find a good<br>
transformation.<br>
<br>
For testing, I established point correspondences manually and used the<br>
itk::LandmarkTransformInitializer to initialize a<br>
VersorRigid3DTransform.<br>
The transform parameters were [0.422982, 0.403292, -0.581769, -27.3493,<br>
155.878, 1790.09] (first the 3 versors, then the 3 translation<br>
parameters). This resulted in a FRE (fiducial registration error) of<br>
about 1.5mm.<br>
<br>
Feeding the same point sets (with correct point correspondences) into<br>
the optimizer resulted in these transform parameters:<br>
[2.08919e-006, 0.00995692, 0.116826, -19.5648, 176.277, 1790.94]<br>
<br>
As you can see, the versors were hardly changed. The<br>
EuclideanDistancePointMetric had the following values:<br>
[56.2236, 51.0406, 28.7475, 109.227, 109.721, 78.9683] (this lead to a<br>
FRE of about 50mm!)<br>
<br>
<br>
Why does the optimizer stop when the cost function is that bad? As<br>
proven with the itk::LandmarkTransformInitializer, a much better<br>
solution exists.<br>
<br>
I suspect that the parameter scaling could be the cause. Does anyone<br>
have a hint of how to set the scales for a VersorRigid3DTransform?<br>
<br>
<br>
My registration code is as following:<br>
<br>
/* lots of type definitions */<br>
typedef itk::PointSet<mitk::ScalarType, 3> PointSetType;<br>
<br>
typedef itk::EuclideanDistancePointMetric< PointSetType, PointSetType><br>
MetricType;<br>
typedef itk::VersorRigid3DTransform< double > TransformType;<br>
typedef TransformType ParametersType;<br>
typedef itk::PointSetToPointSetRegistrationMethod< PointSetType,<br>
PointSetType > RegistrationType;<br>
<br>
MetricType::Pointer metric = MetricType::New();<br>
<br>
TransformType::Pointer transform = TransformType::New();<br>
transform->SetIdentity();<br>
<br>
itk::LevenbergMarquardtOptimizer::Pointer optimizer =<br>
itk::LevenbergMarquardtOptimizer::New();<br>
optimizer->SetUseCostFunctionGradient(false);<br>
<br>
RegistrationType::Pointer registration = RegistrationType::New();<br>
// Scale the translation components of the Transform in the Optimizer<br>
itk::LevenbergMarquardtOptimizer::ScalesType<br>
scales(transform->GetNumberOfParameters());<br>
const double translationScale = 5000;<br>
const double rotationScale = 1.0; // dynamic range of rotations<br>
scales[0] = 1.0 / rotationScale;<br>
scales[1] = 1.0 / rotationScale;<br>
scales[2] = 1.0 / rotationScale;<br>
scales[3] = 1.0 / translationScale;<br>
scales[4] = 1.0 / translationScale;<br>
scales[5] = 1.0 / translationScale;<br>
//scales.Fill(0.01);<br>
unsigned long numberOfIterations = 80000;<br>
double gradientTolerance = 1e-10; // convergence criterion<br>
double valueTolerance = 1e-10; // convergence criterion<br>
double epsilonFunction = 1e-10; // convergence criterion<br>
optimizer->SetScales( scales );<br>
optimizer->SetNumberOfIterations( numberOfIterations );<br>
optimizer->SetValueTolerance( valueTolerance );<br>
optimizer->SetGradientTolerance( gradientTolerance );<br>
optimizer->SetEpsilonFunction( epsilonFunction );<br>
<br>
registration->SetInitialTransformParameters(<br>
transform->GetParameters() );<br>
registration->SetMetric( metric );<br>
registration->SetOptimizer( optimizer );<br>
registration->SetTransform( transform );<br>
registration->SetFixedPointSet( targetPointSet );<br>
registration->SetMovingPointSet( sourcePointSet );<br>
<br>
registration->Update();<br>
<br>
std::cout << "ICP successful: Solution = " <<<br>
transform->GetParameters() << std::endl;<br>
<br>
<br>
Sorry for the long email.<br>
Thank you for your help!<br>
Jochen<br>
_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at: <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
</blockquote></div><br>