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

[Leila baghdadi]

Hi Josh and Jim,

Thanks very much both of you for helping me out, that saves a lot of
hours,

my matrices are 4x4 and because I am implementing this as part of
surface registration in model-based segmentation (i.e, deformable
models), I am a bit concerned about the speed issue that Josh mentioned,
(i.e, there are a fair number of eigen values and vectors which need to
be calculated)

any comments


Leila

On Thu, 2005-04-07 at 12:53, Joshua Cates wrote:
> 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
> > 	>      
> > 	>
> > 	
> > 	
> > 
> >