eigenvalue/vectors for a group of square matrices (in 3D-form?)

Hi, Now I am running into calculating the eigenvalues of Hessian of an image. For every pixel there is the 2x2 matrix that i will need to use eig() on, which involves for-looping on each pixel.

Is there a more elegent way to perform a "group eig" on a 3D matrix of, say 2x2xN without using for loop?

Thanks for any idea and suggestion.

- yi

Answer 1

[Martin]

Yi,

Do you want to calculate the eigenvalues for the 2D Hessian or the 3D Hessian? For the 2D Hessian, you can use the following code to calculate the eigenvalues for all pixels at once.

eigVal1 = (Lxx + Lyy + sqrt((Lxx - Lyy).^2 + 4*Lxy.^2))/2;
eigVal2 = (Lxx + Lyy - sqrt((Lxx - Lyy).^2 + 4*Lxy.^2))/2;

This requires that you have already calculated the second order derivatives of the image and stored them in matrices Lxx, Lxy, Lyy.

[20 Sep 2011 - Edited to add link to 3D Hessian code]

For the 3D Hessian, it seems that there are numerical precision issues with the equivalent analytical formulae for the eigenvalues. I am currently using Dirk-Jan Kroon's code, which looks like an implementation of the QL algorithm ( http://www.mathworks.com/matlabcentral/fileexchange/24409 ). It's written as a Matlab C mex file and is pretty fast. I've just started using it, but it seems to give sensible enough results for my image processing application. You will need to have already calculated matrices containing the relevant 3D partial derivatives at each pixel. Usage is:

[eigVal1, eigVal2, eigVal3, eigVec1, eigVec2, eigVec3] = ...
  eig3volume(Lxx, Lxy, Lxz, Lyy, Lyz, Lzz);

Martin

[Somayeh Ebrahimkhani]

Hi

I have the same problem in the computation of eigenvectors and eigenvalues of 3-by-3 matrix. I tried to download the function (i.e., eig3volume), but the file is an empty.

Can you please write the code of eig3volume here?

Your help is appreciated....

Answer 2

[tomas]

Hi, I asked for something similar 2 weeks ago. Try this. It works fine for my purposes, may be it will work for yours as well.

Q1=rand(2,2,100);

Q2=reshape(Q1,size(Q1,1),size(Q1,3)*size(Q1,2)); Q3=mat2cell(Q2,size(Q1,1),ones(1,size(Q1,3))*size(Q1,2));

[eVect, eVal]=cellfun(@eig,Q3,'un',0);

[Yi]

Thanks for the idea. It does make the code more compact but unfortunately I didn't notice any significant improvement in speed.

For my specific problem (eigenvalue and eigenvector for 3D Hessian matrices) though, there is the analytical solution given here

http://en.wikipedia.org/wiki/Eigenvalue_algorithm#Eigenvalues_of_a_Symmetric_3x3_Matrix

and here:

http://www.groupsrv.com/science/about273201.html

which gives eigenvectors.

Direct calculation shortens the time to about 1/20 compared to eig() in for-loop.

- yi