[Insight-developers] SVD compiled with optimziation

Karthik Krishnan Karthik.Krishnan at kitware.com
Mon Nov 28 20:30:25 EST 2005 

Hello,

I ran into the same problem today, using itk::ImageMomemtsCalculator. (I'm using gcc40 on debian)

Turns out, that uses vnl_svd which uses netlib/dsvdc.c and all these are related.

Funnily enough, turns out that is some numerical issue in netlib/dsvdc when optimization is turned on. 

Did you compile with the -O3 flag (usually the Release mode on CMake) ? For now, I changed the Utilities/vxl/netlib/Makefile to -O2 on the build line that compiles dsvdc.c
and everything was working great.

I haven't snooped around to see what exactly is causing the problem, but there seem to be others in the VXL community who've run into similar issues. 

The following snippet in dsvdc.c and the cvs logs from the vxl sources.

<snip>
----------------------------
revision 1.6
date: 2002/03/13 15:51:28;  author: fsm;  state: Exp;  lines: +12 -2
hack to avoid excess precision
----------------------------
/*
 * Calling this ensures that the operands are spilled to
 * memory and thus avoids excessive precision when compiling
 * for x86 with heavy optimization (gcc). It is better to do
 * this than to turn on -ffloat-store.
 */
static int fsm_ieee_doubles_differ(double *x, double *y)
 {
   return *x != *y;
 }
</snip>


In any case comparing two doubles for equality sounds questionable.

What's bothersome is that the inverse with -O3 builds is not marginally wrong, but quite wrong for a 3x3 rotation matrix ?

Has anybody run into this issue. Is it limited only to gcc ?

I've changed my nightly builds to reflect this. All the associated (23 tests or so fail now on the release builds)
http://www.itk.org/Testing/Sites/Sabbath.kitware/zRel22-Linux-gcc40/20051128-0100-Nightly/Test.html

Thanks
karthik






---------------------------------------------
Suyash, 
 
You might not have enough precision in the matrix values that you are reading from the command line.  Try generating the values with more
precision from whatever tool is generating them.  Try going out to about 12 decimal places.
 
I have seen this before where one program computes a matrix that is supposed to be invertible, writes the matrix in ascii, another program reads that matrix and it is no longer invertible.  
 
Jim
 
 

-----Original Message-----
From: insight-users-bounces+millerjv=crd.ge.com at itk.org <http://www.itk.org/mailman/listinfo/insight-users> [mailto:insight-users-bounces+millerjv=crd.ge.com at itk.org <http://www.itk.org/mailman/listinfo/insight-users>]On Behalf Of Suyash P. Awate
Sent: Friday, November 11, 2005 1:06 PM
To: jiafucang; jjomier at cs.unc.edu <http://www.itk.org/mailman/listinfo/insight-users>; insight-users at itk.org <http://www.itk.org/mailman/listinfo/insight-users>
Subject: [Insight-users] Re: On incorrect matrix inverse


 
Hi everybody,
 
Thanks for the generous help.
 
I was also trying to get the inverse concerning with the registration application, like Fucang (see email attached).
 
Here is my piece of code that gave me the wrong result:
 
------------
const unsigned int Dimension = 3;
typedef double PixelType;
typedef itk::Matrix <PixelType, Dimension, Dimension> MatrixType;
MatrixType matrix;
  matrix.Fill (0);
  matrix (0,0)= atof (argv[4]);
  matrix (0,1)= atof (argv[5]);
  matrix (0,2)= atof (argv[6]);
  matrix (1,0)= atof (argv[7]);
  matrix (1,1)= atof (argv[8]);
  matrix (1,2)= atof (argv[9]);
  matrix (2,0)= atof (argv[10]);
  matrix (2,1)= atof (argv[11]);
  matrix (2,2)= atof (argv[12]);
vnl_matrix_fixed <PixelType, Dimension, Dimension> inverse;
inverse= matrix.GetInverse();
-----------
matrix
 0.99925   0.0371382 -0.0109982
-0.0373772 0.99905   -0.0223916
 0.0101562 0.0227859  0.999689

warnings:
InsightToolkit-2.2.0/Utilities/vxl/core/vnl/algo/vnl_svd.txx: suspicious return value (2) from SVDC
InsightToolkit-2.2.0/Utilities/vxl/core/vnl/algo/vnl_svd.txx: M is 3x3
M = [ ...
 0.9992500000000  0.0371382000000 -0.0109982000000
-0.0373772000000  0.9990500000000 -0.0223916000000
 0.0101562000000  0.0227859000000  0.9996890000000  ]

matrix inverse using vnl in itk (incorrect)
-0.982557 -0.0644525 -0.174434
-0.0677511 0.997618   0.0130162
-0.17318  -0.0246073  0.984583

matrix inverse in matlab (correct)
 0.9992   -0.0374    0.0102
 0.0371    0.9991    0.0228
-0.0110   -0.0224    0.9997
------------
 
I had also tried using the AffineTransform but for some reason could not get that to work and so I wrote this simple matrix code. Perhaps that was due to the same problem that Fucang reported.
Currently, however, my problem is different. As we can see I am not using any Transform class at all. Moreover I am getting warnings from the vnl_svd class, so vnl is definitely involved here.
 
Julien: I am giving 9 numbers (truncated: 0.99925   0.0371382 -0.0109982 -0.0373772 0.99905 -0.0223916 0.0101562 0.0227859  0.999689) to the matrix in itk and also the same numbers to the matrix in matlab. I got these numbers from the output of a separate program.
Just FYI, I have compiled ITK on Suse 10.0 using GCC 4.0. I wonder if that has anything to do with this !
 
Best regards,
Suyash.
 

 
----- Original Message ----- 
From:  <mailto:jiafucang at asisz.com <http://www.itk.org/mailman/listinfo/insight-users>> jiafucang 
To:  <mailto:suyash at cs.utah.edu <http://www.itk.org/mailman/listinfo/insight-users>> suyash at cs.utah.edu <http://www.itk.org/mailman/listinfo/insight-users> ;  <mailto:jjomier at cs.unc.edu <http://www.itk.org/mailman/listinfo/insight-users>> jjomier at cs.unc.edu <http://www.itk.org/mailman/listinfo/insight-users> 
Sent: Friday, November 11, 2005 7:10 AM
Subject: On incorrect matrix inverse


Hello, Suyash, Jomier,
 
I have reported similar problem on
 <http://public.kitware.com/pipermail/insight-users/2005-October/015441.html> http://public.kitware.com/pipermail/insight-users/2005-October/015441.html 
 
and added a bug 2450
 
 <http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24 <http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24>> http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24 <http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24>
 
I think it is because ITK does not computer inverse when setparameters used, and not VNL problem.
 
HTH
 
Fucang

More information about the Insight-developers mailing list