Add files via upload

Hessian_2D.m

@@ -0,0 +1,53 @@
+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

ParabolaRoots.m

@@ -0,0 +1,57 @@
+function roots = ParabolaRoots(varargin)
+% function roots = ParabolaRoots(P)
+%
+% Find roots of second order polynomials
+%
+% INPUT
+%   P: (n x 2) array, each row corresponds to coefficients of each
+%   polynomial, P(:,1)*x^2 + P(:,2)*x + P(:,3)
+% OUTPUT
+%   roots: (n x 2) array, each row correspond to the roots of P
+%
+% To adjust the parameter below which the the discriminant is considerered
+% as nil, use
+%   ParabolaRoots(P, tol)
+% Adjusting tol is useful to avoid the real roots become complex due to
+% numerical accuracy. The default TOL is 0
+%
+% See also: roots, CardanRoots, eig2
+%
+% Author: Bruno Luong <brunoluong@yahoo.com>
+% History:
+%     Original 27-May-2010
+
+% Adjustable parameter
+tol = 0;
+
+if nargin<3
+    P = varargin{1};
+    a = P(:,1);
+    b = P(:,2);
+    c = P(:,3);
+    if nargin>=2
+        tol = varargin{2};
+    end
+else
+    [a b c] = deal(varargin{1:3});
+    if nargin>=4
+        tol = varargin{2};
+    end
+end
+
+if ~isequal(a,1)
+    b = b./a;
+    c = c./a;
+end
+
+b = 0.5*b;
+
+delta = b.^2 - c;
+delta(abs(delta)<tol) = 0;
+sqrtdelta = sqrt(delta);
+
+roots = [sqrtdelta -sqrtdelta];
+roots = bsxfun(@minus, roots, b);
+
+end
+