[Insight-users] Bug in VNL vnl_rnpoly_solve ?

[somi]

Hi,
I compute the rigid transformation between 3D point set v1 and v2
 v2 = s R v1 + T , based on Horn's derivation of the closed form least
squares formulation
(for horns method see :
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.65.971&rep=rep1&type=pdf)

For this I have to solve a polynomial as an intermediate step.
However the function " vnl_rnpoly_solve" included in the ITK (vnl library)
doesn't seem to be returning the correct
roots for a particular case.


a) BUGGY CASE: I used the following points:

+9  -253 -1187
-45 -222 -740
-98 -223 -750

For this the polynomial coefficients would be (while registering above
points to self) :

1.0000e+00   0.0000e+00  -4.0615e+09   4.0000e-06   4.1220e+18

The roots of the polynomial should be
   4.5546e+04
   4.4576e+04
  -4.5546e+04
  -4.4576e+04

But vnl_rnpoly_solve , doesn't return any root.


b) CORRECT CASE:  I then scaled the points by 0.1 to get the modified points

+0.9 -25.3 -118.7
-4.5 -22.2 -74.0
-9.8 -22.3 -75.0

For this the polynomial coefficients would be (while registering above
points to self) :

1.0000e+00   0.0000e+00  -4.0615e+07   0.0000e+00   4.1220e+14

For this the roots should be:
  -4554.6
  -4457.6
   4554.6
   4457.6

Which is what I get using VNL.




This is how I get the coefficients:
a) center the point sets at their respective centroids
b) compute the cross covariance terms of the two centered datasets
c) Form the matrix "N" as in horns method
d)  form the quartic coefficients of the characteristic equation of
     matrix N for solving for the rotation in terms of quaternions
e) Solve the polynomial with co-efficients from step (d)

I have included a test program. To make the code simple I only included the
polynomial solving part in the test code.

Thanks,
Somi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100831/471bb280/attachment.htm>
-------------- next part --------------
# This is the root ITK CMakeLists file.
CMAKE_MINIMUM_REQUIRED(VERSION 2.4)
IF(COMMAND CMAKE_POLICY)
  CMAKE_POLICY(SET CMP0003 NEW)
ENDIF(COMMAND CMAKE_POLICY)


# This project is designed to be built outside the Insight source tree.
PROJECT(VNL_Root_Test)

# Find ITK.
FIND_PACKAGE(ITK REQUIRED)
INCLUDE(${ITK_USE_FILE})

ADD_EXECUTABLE(VNL_Root_Test VNL_Root_Test.cxx )

TARGET_LINK_LIBRARIES(VNL_Root_Test ITKCommon)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: VNL_Root_Test.cxx
Type: text/x-c++src
Size: 3088 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100831/471bb280/attachment.cxx>