Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.
Close
From: somi <seesomi@gm...>  20100901 20:29:38
Attachments:
Message as HTML
VNL_Root_Test.cxx

# 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) 
From: Peter Vanroose <peter_vanroose@ya...>  20100902 13:41:09
Attachments:
Message as HTML

Somi, I have the impression that this has to do with the precision and size of your coefficients, and more precisely of the magnitude difference between them. Your constant coefficient is around 4e18. If I scale the solutions by a factor 10, by dividing coefficient 3 by 100 and coefficient 5 by 10000, I obtain the correct result (to be interpreted 10x larger): ======================== VNL roots are 4554.6 +i 4.62223e32 4457.61 +i 1.44445e34 4554.6 +i 4.62223e32 4457.61 +i 1.44445e34 ======================== So try reducing too extreme values first (e.g. by dividing coefficient n by 10^n) before calling the solver. I presume that the original algorithm (by Kriegman and Ponce) has this same instability. Not sure whether this scaling solution should be introduced into the VNL interface layer, or whether this should better be left to the user who calls the method. Suggestions welcome, of course.  Peter.   Den ons 20100901 skrev somi <seesomi@...>: Från: somi <seesomi@...> Ämne: [Vxlusers] Bug in VNL vnl_rnpoly_solve ? Till: vxlusers@... Datum: onsdag 1 september 2010 22:29 Hi,I am getting incorrect solution when I try to solve a polynomial using the method "vnl_rnpoly_solve",that is provided in ITK. a) Polynomial coefficients (buggy case ) 1.0000e+00 0.0000e+00 4.0615e+09 4.0000e06 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) Polynomial coefficients (test case) 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. I have included a test program. Thanks,Somi Infogad bilaga följer  This SF.net Dev2Dev email is sponsored by: Show off your parallel programming skills. Enter the Intel(R) Threading Challenge 2010. http://p.sf.net/sfu/intelthreadsfd Infogad bilaga följer _______________________________________________ Vxlusers mailing list Vxlusers@... https://lists.sourceforge.net/lists/listinfo/vxlusers 
From: somi <seesomi@gm...>  20100902 13:50:48
Attachments:
Message as HTML

Hi Peter, Thanks for the reply. I checked the solution in atleast two programs (octave >roots and Numerical Recipies>zroots) and they give the correct solution for these coefficients. You are right, this has something to do with numerical instabilities during scaling. Unfortunately I can't change the coefficients. Let me explain the context in which I am using these coefficients. I have to compute the rigid transformation between two 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.0000e06 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) TEST 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 coefficients from step (d) Is there any way around it ? I could scale down points by 10, and again scale back the results. But I want to get a generic solution as I use this program in a library and I would have no way of telling automatically when to scale and when not to scale. Thanks, Somi On Thu, Sep 2, 2010 at 9:41 AM, Peter Vanroose <peter_vanroose@...>wrote: > Somi, > > I have the impression that this has to do with the precision and size of > your coefficients, > and more precisely of the magnitude difference between them. > Your constant coefficient is around 4e18. > If I scale the solutions by a factor 10, by dividing coefficient 3 by 100 > and coefficient 5 by 10000, I obtain the correct result (to be interpreted > 10x larger): > ======================== > VNL roots are > 4554.6 +i 4.62223e32 > 4457.61 +i 1.44445e34 > 4554.6 +i 4.62223e32 > 4457.61 +i 1.44445e34 > ======================== > > So try reducing too extreme values first (e.g. by dividing coefficient n by > 10^n) before calling the solver. > > I presume that the original algorithm (by Kriegman and Ponce) has this same > instability. > > Not sure whether this scaling solution should be introduced into the VNL > interface layer, or whether this should better be left to the user who calls > the method. Suggestions welcome, of course. > >  Peter. > > >  > >  Den *ons 20100901 skrev somi <seesomi@...>*: > > > Från: somi <seesomi@...> > Ämne: [Vxlusers] Bug in VNL vnl_rnpoly_solve ? > Till: vxlusers@... > Datum: onsdag 1 september 2010 22:29 > > > Hi, > I am getting incorrect solution when I try to solve a polynomial using the > method "vnl_rnpoly_solve", > that is provided in ITK. > > a) Polynomial coefficients (buggy case ) > > 1.0000e+00 0.0000e+00 4.0615e+09 4.0000e06 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) Polynomial coefficients (test case) > > > 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. > > > > I have included a test program. > > Thanks, > Somi > > Infogad bilaga följer > > >  > This SF.net Dev2Dev email is sponsored by: > > Show off your parallel programming skills. > Enter the Intel(R) Threading Challenge 2010. > http://p.sf.net/sfu/intelthreadsfd > > Infogad bilaga följer > > _______________________________________________ > Vxlusers mailing list > Vxlusers@...<http://mc/compose?to=Vxlusers@...>; > https://lists.sourceforge.net/lists/listinfo/vxlusers > > > 
From: Peter Vanroose <peter_vanroose@ya...>  20100902 15:11:11
Attachments:
Message as HTML

> Is there any way around it ? I could scale down points by 10, > and again scale back the results. But I want to get a generic > solution as I use this program in a library and I would have > no way of telling automatically when to scale and when not to scale. Yes, of course, I see the problem. A "generic" way would be to compare the first and last nonzero coefficient, and scale by a factor which makes them approximately equal (i.e., scale by the nthe root of that factor), then scale back the result. I was thinking of adding a method "solve_robust()" to the vnl_rnpoly_solve class which does something like this, but I'm not feeling comfortable enough to do this right away for the general mdimensional case. But for your particular (generic) use, it should not be too difficult: something like: double factor = coeffs[4]/coeffs[0]; // take 4th root: factor = sqrt(factor); factor = sqrt(factor); double mfactor = 1.0; for (int i=0; i<=4; ++i) coeffs[i]/=mfactor, mfactor*=factor; ... vnl_rnpoly_solve solver(l); vcl_vector<vnl_vector<double>*> re = solver.real(); vcl_vector<vnl_vector<double>*> im = solver.imag(); ... (*(*rp)) *= factor; // globally scale the vnl_vector back I've put a working implementation of exactly this in the last SVN version of core/vnl/algo/tests/test_rnpoly_roots.cxx  Peter. 
From: somi <seesomi@gm...>  20100902 15:39:38
Attachments:
Message as HTML

Thanks Peter. That seems to work for me ( I tried bunch of 4th order polynomials, and I got the correct solutions) On Thu, Sep 2, 2010 at 11:11 AM, Peter Vanroose <peter_vanroose@...>wrote: > > Is there any way around it ? I could scale down points by 10, > > and again scale back the results. But I want to get a generic > > solution as I use this program in a library and I would have > > no way of telling automatically when to scale and when not to scale. > > Yes, of course, I see the problem. > > A "generic" way would be to compare the first and last nonzero coefficient, > and scale by a factor which makes them approximately equal (i.e., scale by > the nthe root of that factor), then scale back the result. > > I was thinking of adding a method "solve_robust()" to the vnl_rnpoly_solve > class which does something like this, but I'm not feeling comfortable enough > to do this right away for the general mdimensional case. But for your > particular (generic) use, it should not be too difficult: something like: > > double factor = coeffs[4]/coeffs[0]; > // take 4th root: > factor = sqrt(factor); factor = sqrt(factor); > double mfactor = 1.0; > for (int i=0; i<=4; ++i) coeffs[i]/=mfactor, mfactor*=factor; > ... > vnl_rnpoly_solve solver(l); > vcl_vector<vnl_vector<double>*> re = solver.real(); > vcl_vector<vnl_vector<double>*> im = solver.imag(); > ... > (*(*rp)) *= factor; // globally scale the vnl_vector back > > I've put a working implementation of exactly this in the last SVN version > of > core/vnl/algo/tests/test_rnpoly_roots.cxx > > >  Peter. > > > > > > > > > > > > 