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