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
