[Insight-users] calculating eigen values and vector of matrix

[Joshua Cates]

Hi Leila,

You can find some additional examples of using vnl_eigensystem in
itkVectorGradientMagnitudeImageFilter.txx/.h.  Note that for 3x3 real
symmetric matrices I chose to code my own solver since the vnl_eigensystem
routines are not thread safe and pretty generic (i.e. slow).

Josh.


On Thu, 7 Apr 2005, Miller, James V (Research) wrote:

> Leila, 
>  Here is an example of using the vnl routines.  The code a eigenvalue
> decomposition of a matrix M and extracts the eigenvector corresponding
> to the smallest eigenvalue.  I use the vnl_real_eigensystem to compute
> the eigensystem for the real matrix M.  The results (eigenvalues and
> eigenvectors) are always "returned" as complex numbers matrices of
> complex numbers (eigenvalues stored in a vnl_diag_matrix).
>  To start this process with an itk::Matrix, you should be able to just
> the result of GetVnlMatrix() to the eigensolver.
>  
> #include <vnl/vnl_matrix.h>
> #include <vnl/vnl_math.h>
> #include <vnl/vnl_vector.h>
> #include <vnl/algo/vnl_real_eigensystem.h>
> 
> {
>   vnl_matrix<double> M(3, 3);
>   vnl_vector<double> params;
>   double lambda;
> 
>   M = .......;  // some matrix M
>  
>   // construct the eigenvector problem
>   vnl_real_eigensystem eig(M);
> 
>   // extract eigenvector corresponding to the smallest positive real eigenvalue
>   vnl_diag_matrix<vcl_complex<double> >::iterator it;
>   vnl_diag_matrix<vcl_complex<double> >::iterator small_it = eig.D.end();
>   double small = vcl_numeric_limits<double>::max();
>   const double eigenvalue_zero_tol = 1.0e-11;
>   
>   for (it = eig.D.begin(); it != eig.D.end(); ++it)
>     {
>     if ((*it).real() > eigenvalue_zero_tol
>         && vcl_fabs((*it).imag()) < eigenvalue_zero_tol
>         && (*it).real() < small)
>       {
>       small_it = it;
>       }
>     }
> 
>   params.set_size( eig.V.rows() ); // destroys previous values
>   for (i=0; i < eig.V.rows(); i++)
>     {
>     params(i) = eig.V(i, small_it - eig.D.begin()).real();
>     }
>   lambda = (*small_it).real();
> }
> 
> 	-----Original Message----- 
> 	From: Leila baghdadi [mailto:baghdadi at sickkids.ca] 
> 	Sent: Wed 4/6/2005 4:44 PM 
> 	To: Miller, James V (Research) 
> 	Cc: 
> 	Subject: RE: [Insight-users] calculating eigen values and vector of matrix
> 	
> 	
> 
> 	Hello Jim,
> 	
> 	I see what you mean but I am having difficulty using the vnl methods
> 	from my class,
> 	
> 	can you please give me a simple example,
> 	
> 	
> 	Thanks
> 	
> 	Leila
> 	
> 	On Wed, 2005-04-06 at 16:03, Miller, James V (Research) wrote:
> 	> Leila,
> 	> 
> 	> You'll have to call GetVnlMatrix() on your itk::Matrix and pass that to a vnl routine. There are several vnl eigensystem methods to choose from depending on the structure of your matrix (real, symmetric, ...)
> 	> 
> 	> Jim
> 	>
> 	>       -----Original Message-----
> 	>       From: insight-users-bounces at itk.org on behalf of Leila baghdadi
> 	>       Sent: Wed 4/6/2005 3:52 PM
> 	>       To: Insight-users at itk.org
> 	>       Cc:
> 	>       Subject: [Insight-users] calculating eigen values and vector of matrix
> 	>      
> 	>      
> 	>
> 	>       Hi everyone,
> 	>      
> 	>       I know that the vnl utilities should have methods for the above
> 	>      
> 	>       but can you just use itkMatrix,
> 	>      
> 	>       I do not see any public methods for calculating eigen values or vector
> 	>      
> 	>      
> 	>       any hints
> 	>      
> 	>       thanks
> 	>      
> 	>       Leila
> 	>      
> 	>       _______________________________________________
> 	>       Insight-users mailing list
> 	>       Insight-users at itk.org
> 	>       http://www.itk.org/mailman/listinfo/insight-users
> 	>      
> 	>
> 	
> 	
> 
> 

-- 
______________________________
 Josh Cates			
 Scientific Computing and Imaging Institute
 University of Utah
 http://www.sci.utah.edu/~cates