[vtkusers] ICP Error: vtkMath::Jacobi: Error extracting eigenfunctions

[David Gobbi]

Do you think that JacobiN was only failing when there was a NAN in the
input, or was it failing under more normal circumstances too?

 - David


2012/2/6 Dženan Zukić <dzenanz at gmail.com>:
> I think it is useful to put this finding on the mailing list for future
> reference.
>
> On Mon, Feb 6, 2012 at 18:28, Maarten Beek <beekmaarten at yahoo.com> wrote:
>>
>> Hate to admit, but I found a bug in my code as well. Same bug as Dzenan.
>> Both vtkPoints objects were created in similar for-loops. So one of the
>> arrays was not totally filled. Gives very unpredictable behavior (sometimes
>> not even crashing).
>> Switching the point sets however always caused crashes and should have
>> been the clue I shouldn't have missed.
>>
>> Remain the facts that vtkMath::vtkJacobiN doesn't really handle this error
>> nicely and that the algorithm used is potentially unstable.
>>
>> Maarten
>>
>> ________________________________
>> From: Dženan Zukić <dzenanz at gmail.com>
>> To: David Gobbi <david.gobbi at gmail.com>
>> Cc: Maarten Beek <beekmaarten at yahoo.com>
>> Sent: Friday, February 3, 2012 10:55:58 AM
>>
>> Subject: Re: [vtkusers] ICP Error: vtkMath::Jacobi: Error extracting
>> eigenfunctions
>>
>> In trying to create a self-contained test program, I discovered a bug in
>> my own code. I was supplying some junk points to the vtkLandmarkTransform
>> (one points array not completely filled with values).
>>
>> Now the execution in release mode doesn't give any warnings either. Maybe
>> something similar happened to you Maarten?
>>
>> 2012/2/1 David Gobbi <david.gobbi at gmail.com>
>>
>> Would it be possible for you to send me a short, simple program that
>> includes vtkLandmarkTransform that triggers the error?  Or can you try
>> to modify vtkMath::JacobiN in such a way that the error goes away,
>> using the Numerical Recipes discussion on the following link as a
>> guideline? http://www.nr.com/forum/showthread.php?p=4913
>>
>> There isn't much that I can do on my side unless I have some way of
>> reproducing the error...
>>
>> 2012/2/1 Dženan Zukić <dzenanz at gmail.com>:
>> > Same: problem is not encountered.
>> >
>> > 2012/2/1 David Gobbi <david.gobbi at gmail.com>
>> >>
>> >> Can you modify my test program so that main() uses vtkMath::JacobiN
>> >> instead of using its own JacobiN, and then try it again?
>> >>
>> >> 2012/2/1 Dženan Zukić <dzenanz at gmail.com>:
>> >> > This program does not crash in either debug or release mode. The
>> >> > optimization options are the same (as CMake sets them).
>> >> >
>> >> > I have modified it to read input for the files I sent you, but it
>> >> > goes
>> >> > through them without errors. I guess that whatever is happening, is
>> >> > not
>> >> > triggered in this simplified example.
>> >> >
>> >> > 2012/1/31 David Gobbi <david.gobbi at gmail.com>
>> >> >>
>> >> >> I haven't been able to reproduce the problem, but at work all I have
>> >> >> is
>> >> >> a Mac and the optimizations it does are probably different from
>> >> >> MSVC.
>> >> >>
>> >> >> Can you try different matrices with the attached program and see if
>> >> >> you can get it to fail?
>> >> >>
>> >> >>  - David
>> >> >>
>> >> >>
>> >> >> 2012/1/31 Dženan Zukić <dzenanz at gmail.com>:
>> >> >> > Here is an output from non-release version, same input as before.
>> >> >> >
>> >> >> > 2012/1/31 Dženan Zukić <dzenanz at gmail.com>
>> >> >> >>
>> >> >> >> Here are 3 different runs of my program, identical input.
>> >> >> >>
>> >> >> >> 2012/1/31 David Gobbi <david.gobbi at gmail.com>
>> >> >> >>>
>> >> >> >>> Thanks for the info.  Can you post your compiler optimization
>> >> >> >>> flags
>> >> >> >>> to
>> >> >> >>> the list?  There is definitely a problem with JacobiN that needs
>> >> >> >>> to
>> >> >> >>> be
>> >> >> >>> fixed.  The fix should be very easy to do if we can get a matrix
>> >> >> >>> that
>> >> >> >>> reliably causes JacobiN to fail when certain optimization flags
>> >> >> >>> are
>> >> >> >>> used, so if either of you can print out such a matrix (with all
>> >> >> >>> numbers printed with a precision of 20 digits), it would be a
>> >> >> >>> big
>> >> >> >>> help.
>> >> >> >>>
>> >> >> >>>  - David
>> >> >> >>>
>> >> >> >>>
>> >> >> >>> 2012/1/31 Dženan Zukić <dzenanz at gmail.com>:
>> >> >> >>> > I also saw this warning recently, from the block that "never
>> >> >> >>> > gets
>> >> >> >>> > called".
>> >> >> >>> > And it only happened in release version of my program.
>> >> >> >>> >
>> >> >> >>> > On Tue, Jan 31, 2012 at 00:48, Maarten Beek
>> >> >> >>> > <beekmaarten at yahoo.com>
>> >> >> >>> > wrote:
>> >> >> >>> >>
>> >> >> >>> >> Hi David,
>> >> >> >>> >>
>> >> >> >>> >> Ai, seems these crashes are hard to get rid off....
>> >> >> >>> >>
>> >> >> >>> >> I currently use VTK within Amira modules.
>> >> >> >>> >>
>> >> >> >>> >> I changed the type of vtkPoints from the default VTK_FLOAT to
>> >> >> >>> >> VTK_DOUBLE
>> >> >> >>> >> when I convert an Amira object to a VTK object. This reduces
>> >> >> >>> >> the
>> >> >> >>> >> occurrences
>> >> >> >>> >> of the errors, but doesn't solve the issue.
>> >> >> >>> >> (Ps: why is the default in vtkPoints still float, while
>> >> >> >>> >> everywhere
>> >> >> >>> >> else
>> >> >> >>> >> floats were changed into doubles a long while ago?)
>> >> >> >>> >>
>> >> >> >>> >> I don't really know much about different optimizations, so I
>> >> >> >>> >> just
>> >> >> >>> >> copy-paste what I found in MSVS:
>> >> >> >>> >>
>> >> >> >>> >> I build my VTK with:
>> >> >> >>> >> /O2 /Ob2 /FD /EHsc /MD /Fo"vtkCommon.dir\Release\\"
>> >> >> >>> >> /Fd"G:\vtk-5.6.0\build_vs90_noQt\bin\Release/vtkCommon.pdb"
>> >> >> >>> >> /W4
>> >> >> >>> >> /nologo /c
>> >> >> >>> >> /TP /errorReport:prompt/Zm 1000 and a couple of '/I'-s and
>> >> >> >>> >> '/D'-s
>> >> >> >>> >>
>> >> >> >>> >> I build my Amira modules with:
>> >> >> >>> >> /O2 /Ob1 /Oi /Oy /GF /FD /EHsc /MD /GS- /Gy /fp:fast /openmp
>> >> >> >>> >> /Fo"./../../obj/arch-Win64VC9-Optimize/GraphPackage/"
>> >> >> >>> >> /Fd"./../.././bin/arch-Win64VC9-Optimize/GraphPackage.pdb"
>> >> >> >>> >> /W3
>> >> >> >>> >> /nologo
>> >> >> >>> >> /c
>> >> >> >>> >> /Zi /wd4068 /wd4305 /wd4018 /wd4244 /wd4800 /wd4275 /wd4251
>> >> >> >>> >> /errorReport:prompt and a couple of '/I'-s and '/D'-s
>> >> >> >>> >>
>> >> >> >>> >> I am working on getting vtkMath::JacobiN to print the matrix
>> >> >> >>> >> in
>> >> >> >>> >> the
>> >> >> >>> >> if-loop that 'gets never called' using vtkGenericWarningMacro
>> >> >> >>> >> to
>> >> >> >>> >> avoid
>> >> >> >>> >> having to include Amira stuff. Having a hard time getting
>> >> >> >>> >> that
>> >> >> >>> >> to
>> >> >> >>> >> work. (
>> >> >> >>> >> char* test = "hello";vtkGenericWarningMacro( test ); won't do
>> >> >> >>> >> the
>> >> >> >>> >> trick).
>> >> >> >>> >> Any suggestions?
>> >> >> >>> >>
>> >> >> >>> >> Maarten
>> >> >> >>> >>
>> >> >> >>> >> ________________________________
>> >> >> >>> >> From: David Gobbi <david.gobbi at gmail.com>
>> >> >> >>> >> To: Maarten Beek <beekmaarten at yahoo.com>
>> >> >> >>> >> Cc: VTK list <vtkusers at vtk.org>
>> >> >> >>> >> Sent: Monday, January 30, 2012 2:31:38 PM
>> >> >> >>> >> Subject: Re: [vtkusers] ICP Error: vtkMath::Jacobi: Error
>> >> >> >>> >> extracting
>> >> >> >>> >> eigenfunctions
>> >> >> >>> >>
>> >> >> >>> >> Hi Maarten,
>> >> >> >>> >>
>> >> >> >>> >> I think there are two issues here, one is that
>> >> >> >>> >> vtkLandmarkTransform
>> >> >> >>> >> should have a check as you have described, but another
>> >> >> >>> >> potentially
>> >> >> >>> >> serious issue is that, as far as I understand it, JacobiN
>> >> >> >>> >> should
>> >> >> >>> >> never
>> >> >> >>> >> fail.  I did a google search and found a message from 2002 by
>> >> >> >>> >> the
>> >> >> >>> >> Numerical Recipes authors stating that the NR Jacobi code can
>> >> >> >>> >> become
>> >> >> >>> >> numerically unstable due to certain compiler optimizations.
>> >> >> >>> >> I'm
>> >> >> >>> >> guessing that certain compiler optimizations could similarly
>> >> >> >>> >> make
>> >> >> >>> >> the
>> >> >> >>> >> VTK Jacobi code become unstable.
>> >> >> >>> >> http://www.nr.com/forum/showthread.php?p=4913
>> >> >> >>> >>
>> >> >> >>> >> Can you provide info about your compiler/flags etc?  Also, is
>> >> >> >>> >> it
>> >> >> >>> >> possible to get your code to print a matrix that causes
>> >> >> >>> >> JacobiN
>> >> >> >>> >> to
>> >> >> >>> >> fail?
>> >> >> >>> >>
>> >> >> >>> >> - David
>> >> >> >>> >>
>> >> >> >>> >>
>> >> >> >>> >> On Mon, Jan 30, 2012 at 9:44 AM, Maarten Beek
>> >> >> >>> >> <beekmaarten at yahoo.com>
>> >> >> >>> >> wrote:
>> >> >> >>> >> > Hi all,
>> >> >> >>> >> >
>> >> >> >>> >> > I get the following warning when using ICP:
>> >> >> >>> >> >
>> >> >> >>> >> > vtkMath::Jacobi: Error extracting eigenfunctions
>> >> >> >>> >> >
>> >> >> >>> >> > It comes from vtkMath::vtkJacobiN that calculates the
>> >> >> >>> >> > eigenvalues,
>> >> >> >>> >> > -vectors
>> >> >> >>> >> > of a matrix.
>> >> >> >>> >> > The function returns 0, when this error occurs (otherwise
>> >> >> >>> >> > 1),
>> >> >> >>> >> > but
>> >> >> >>> >> > vtkLandmarkTransform doesn't check for the return value and
>> >> >> >>> >> > thus
>> >> >> >>> >> > uses
>> >> >> >>> >> > undefined values for the eigenvalues, -vectors.
>> >> >> >>> >> > This means my point set becomes invalid.
>> >> >> >>> >> > I believe the way to avoid this error is to 'giggle' the
>> >> >> >>> >> > coordinates
>> >> >> >>> >> > a
>> >> >> >>> >> > little when it happens, but how can I intervene when
>> >> >> >>> >> > vtkLandmarkTransform
>> >> >> >>> >> > stubbornly uses whatever values it gets from
>> >> >> >>> >> > vtkMath::vtkJacobiN?
>> >> >> >>> >> > Is there a way to intercept the error message?
>> >> >> >>> >> >
>> >> >> >>> >> > I also noticed that if I apply ICP on A and B1 resulting in
>> >> >> >>> >> > the
>> >> >> >>> >> > error
>> >> >> >>> >> > message, applying ICP on A and B2 (loading the same data
>> >> >> >>> >> > file
>> >> >> >>> >> > that
>> >> >> >>> >> > gives
>> >> >> >>> >> > me
>> >> >> >>> >> > B1) works fine. I don't see vtkMath::vtkJacobiN using a
>> >> >> >>> >> > random
>> >> >> >>> >> > number or
>> >> >> >>> >> > doing something similar, which makes me believe the error
>> >> >> >>> >> > message
>> >> >> >>> >> > I
>> >> >> >>> >> > get
>> >> >> >>> >> > is
>> >> >> >>> >> > due to floating point accuracy? Even more a reason to think
>> >> >> >>> >> > that
>> >> >> >>> >> > 'giggling'
>> >> >> >>> >> > might work (if I were able to intervene before the invalid
>> >> >> >>> >> > transform
>> >> >> >>> >> > is
>> >> >> >>> >> > applied to my data...)
>> >> >> >>> >> >
>> >> >> >>> >> > Thanks - Maarten
>> >> >> >>
>> >> >> >>
>> >> >> >
>> >> >
>> >> >
>> >
>> >
>>
>>
>>
>>
>