[Insight-users] Re: On incorrect matrix inverse

[Suyash P. Awate]

Dear ITK Users,

I found the cause of the error. It deals with my build of ITK in the SUSE 10 environment that uses gcc 4.0 by default. When I build itk and the applications on SUSE 9.3 using gcc 3.4 (definitely not gcc 4.0) or Fedora Core 2 (again not using gcc 4.0), I did NOT get any vnl error.

I confirmed that precision is NOT an issue because the condition number for the matrix with truncated values is still very close to 1. I found this via matlab. So numerical "instabilities" in inverse-finding should not occur.

So it seems that there is something strange about the way the vnl code/libraries interact with gcc 4.0.

Thanks again for all your help.

Best regards,
Suyash.
  ----- Original Message ----- 
  From: Miller, James V (Research) 
  To: Suyash P. Awate ; jiafucang ; jjomier at cs.unc.edu ; insight-users at itk.org 
  Sent: Friday, November 11, 2005 11:24 AM
  Subject: RE: [Insight-users] Re: On incorrect matrix inverse


  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 [mailto:insight-users-bounces+millerjv=crd.ge.com at itk.org]On Behalf Of Suyash P. Awate
    Sent: Friday, November 11, 2005 1:06 PM
    To: jiafucang; jjomier at cs.unc.edu; insight-users at itk.org
    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: jiafucang 
      To: suyash at cs.utah.edu ; jjomier at cs.unc.edu 
      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 

      and added a bug 2450

      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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20051219/cd8a998e/attachment.htm