[Insight-users] Re: On incorrect matrix inverse

[Suyash P. Awate]

I forgot to mention that I always have the -O3 optimization flag ON, while compiling.

Regards,
Suyash.
  ----- Original Message ----- 
  From: Suyash P. Awate 
  To: Miller, James V (Research) ; jiafucang ; jjomier at cs.unc.edu ; insight-users at itk.org 
  Sent: Monday, December 19, 2005 8:50 PM
  Subject: Re: [Insight-users] Re: On incorrect matrix inverse


  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/20051220/bb16f498/attachment.html