Hessian_2D.m

function [lambda1,lambda2] = Hessian_2D(image)

% =========================================================================
% sub function called by iterative particle reconstruction
% calculate the Hessian maxtrix of a 2D image
% and its eigenvalues for filtering.
% by Tianqi Guo, Aug 2017
% =========================================================================

% display('Hessian filtering')

% addpath('Eig3Folder/');

% Gaussian smoothing parameters for eigenvalue matrix
filter_sd = 0.65; 

% first order derivatives of intensity volume
gx = socdiff(image,[],2);
gy = socdiff(image,[],1);

% second order derivatives of intensity volume
gxx = socdiff(gx,[],2);
gxy = socdiff(gx,[],1);

gyx = socdiff(gy,[],2);
gyy = socdiff(gy,[],1);

% eigenvalues of the Hessian matrix at each voxel
% sorted as absolute values lambda1 < lambda2

H = cat(4,gxx,gxy,gyx,gyy);

H = reshape(H,[size(image(:),1),2,2]);

H = permute(H,[3 2 1]);

% eigenvalue for Hessian matrix and its absolute value
e = eig2(H);
e_abs = abs(e);

[~, ind] = sort(e_abs);
[m,n] = size(e);
e_sort = real(e(bsxfun(@plus,ind,(0:n-1)*m)));

lambda1 = reshape(e_sort(1,:),size(image));
lambda2 = reshape(e_sort(2,:),size(image));

% Gaussian smoothing of eigenvalue matrices to suppress noise
lambda1 = imgaussfilt(lambda1,filter_sd);
lambda2 = imgaussfilt(lambda2,filter_sd);


end