[Insight-users] [Insight-developers] instantiating itkVector : How to solve a linear system with VNL

[Luis Ibanez]

Hi Richard,


Here are some options for solving linear systems with VNL:


         For the problem :     A x = b


         Where "A" and "b" are known,

         and you want to solve for "x".


You can use:


A) The expression:

       #include "vnl_qr.h"
       x = vnl_qr< double >( A ).solve( b );


or


B)   #include "vnl_svd.h"
      vnl_svd< double > svd(A);
      vnl_matrix<double> invA = svd.inverse();

      x = invA * b;


or

C)  #include "vnl_svd.h"
     x = vnl_svd< double >( A ).solve( b );


Depending on the characteristics of your matrix,
and its size you may get different performance
with the methods above.


You may want to try them all and profile them
to see what gives you better performance.


Please let us know what you find.

--

BTW: Just for the record, please note the
      vnl_inverse is only implemented for
      matrix sizes up to 4x4.



      Thanks


         Luis



---------------------
Richard Beare wrote:
> Hi,
> I've started playing with the vnl routines, which is convenient for
> the moment. I'll see about optimizing later if it seems worthwhile.
> 
> One question though - I'm now at the point of wanting to solve Ax=b ,
> with small, non sparse matrices. It seems that I should be able to use
> vnl_linear_system (in order to avoid inverting A explicitly), but I
> can't find any examples of how to use it. Any suggestions?
> 
> Thanks
> 
> On Sun, Dec 7, 2008 at 1:02 PM, Richard Beare <richard.beare at gmail.com> wrote:
> 
>>It is difficult to be certain about performance issues at this stage,
>>but I can make some guesses. The critical issue is having matrix maths
>>available.
>>
>>I'm actually coding up the SIFT algorithms - there was an IJ
>>publication about this a while back, but it doesn't appear
>>particularly good.
>>
>>The matrix maths part is applied on a per feature basis, rather than a
>>per voxel basis, so I don't expect it to be a big part of the
>>execution time.
>>
>>Is there a variable length itk vector that can be used in conjunction
>>with a variable size matrix while allowing matrix algebra?
>>
>>On Sun, Dec 7, 2008 at 9:55 AM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>>
>>>Richard,
>>>
>>>You may have to consider how visible this class is supposed to be
>>>and what performance implications you can anticipate from its
>>>construction (e.g. is the type of class that you will use to
>>>create thousands of instances in an inner loop ?).
>>>
>>>Using one of vcl vector classes that you can resize at run-time
>>>is definitely an option. It is worth to write a small use case
>>>example to see which one will work better.
>>>
>>>Note as well that ITK has the variable length vector, whose
>>>size can be defined at run time.
>>>
>>>
>>>  Regards,
>>>
>>>
>>>      Luis
>>>
>>>--------------------
>>>Richard Beare wrote:
>>>
>>>>Am I likely to be better off using underlying vcl classes directly?
>>>>Will that give me more flexibility?
>>>>
>>>>On Sat, Dec 6, 2008 at 4:26 AM, Luis Ibanez <luis.ibanez at kitware.com>
>>>>wrote:
>>>>
>>>>
>>>>>Hi Richard,
>>>>>
>>>>>Thanks for the clarification.
>>>>>
>>>>>
>>>>>This is good news and bad news  :-)
>>>>>
>>>>>
>>>>>
>>>>>The good news is that you know the number at compile time,
>>>>>since that is what the itkVector instantiation will require.
>>>>>
>>>>>The bad news is that the expression that you need to evaluate
>>>>>for computing the vector size is a non-trivial one.
>>>>>
>>>>>
>>>>>One potential way of generating that computation would be to
>>>>>have a Traits helper class, at the image of the NumericTraits
>>>>>class.
>>>>>
>>>>>In this helper class you could define explicit instantiations
>>>>>such that an internal enum (or static int) will be initialized
>>>>>to the value of 3^(N-1).
>>>>>
>>>>>This is a bit manual, but I'm assuming that you would not need
>>>>>more that dimensions in the range N=1,2,3,4,5 ... (?).
>>>>>So you will be done with the Traits in a couple of copy/pasting
>>>>>sessions, and in less than 50 lines of code.
>>>>>
>>>>>There are probably smarter ways of getting this to work,
>>>>>but more coffee may be required for figuring them out...
>>>>>
>>>>>
>>>>>  Regards,
>>>>>
>>>>>
>>>>>      Luis
>>>>>
>>>>>
>>>>>
>>>>>----------------------
>>>>>Richard Beare wrote:
>>>>>
>>>>>
>>>>>>Hi,
>>>>>>
>>>>>>In this case I'll know it at compile time - it is dependent on the
>>>>>>dimensionality of the image. Actually, to be specific, it is given by
>>>>>>3^(ImType::ImageDimension - 1).
>>>>>>
>>>>>>It would probably look neater in the code if I could do it in a way
>>>>>>that was not dependent on compile time knowledge, but I'm not sure
>>>>>>that is possible. I'll need to do something equivalent with matrix
>>>>>>types too.
>>>>>>
>>>>>>On Fri, Dec 5, 2008 at 12:50 AM, Luis Ibanez <luis.ibanez at kitware.com>
>>>>>>wrote:
>>>>>>
>>>>>>
>>>>>>
>>>>>>>Richard,
>>>>>>>
>>>>>>>
>>>>>>>Would you know the size of the neighborhood at compilation time ?
>>>>>>>
>>>>>>>or only at run-time ?
>>>>>>>
>>>>>>>
>>>>>>>Please let us know,
>>>>>>>
>>>>>>>
>>>>>>>Thanks
>>>>>>>
>>>>>>>
>>>>>>> Luis
>>>>>>>
>>>>>>>
>>>>>>>----------------------
>>>>>>>Richard Beare wrote:
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>>Hi,
>>>>>>>>I need to create a vector type with length equal to the size of a
>>>>>>>>neighborhood - i.e. 3^ImageDimension. I suspect that there needs to be
>>>>>>>>some sort of trick with macros, because calling pow is no good. Any
>>>>>>>>suggestions?
>>>>>>>>
>>>>>>>>Thanks
>>>>>>>>_______________________________________________
>>>>>>>>Insight-developers mailing list
>>>>>>>>Insight-developers at itk.org
>>>>>>>>http://www.itk.org/mailman/listinfo/insight-developers
>>>>>>>>
>>>>>>>
>