sdgd_filt2D(data_in, x, mu, sigma, sigma_width)

function [data_out, kernel] = sdgd_filt2D(data_in, x, mu, sigma, sigma_width)
% Function [data_out, kernel] = sdgd_filt2D(data_in, x, mu, sigma, sigma_width)
% Second-Derivative-in-the-Gradient-Direction edge detection filter with 2D kernel.
% This function convolves input image with set 2D Gaussian kernels then combines 
% results in accordance with formula given in reference paper.
%
%INPUTS:
%   DATA_IN     - input image
%   X           - size of Gaussian kernel.
%   MU          - mean of the Gaussian
%   SIGMA       - standard deviation of a Gaussian function
%   SIGMA_WIDTH - defines where to cut the Gaussian kernel tail (or width of the kernel in sigma).
%   Bigger number- more of Gaussian included in the kernel. Defaul value= 3 or 98% of Gaussian 
%
%OUTPUTS:
%   DATA_OUT- filtered data output
%   KERNEL - stucture with derivative kernels used for filtering:
%   {kernel.dGx, kernel.dGxx, kernel.dGxy}
%
%WARNING:   Derivatives for some pixels may be zero. This produce an error "divideByZero" and such 
%           pixels gets value =NaN. It is recomended to check for NaNs after filtering and replacing
%           it with value of convenience -> zero, global minimum,  etc...
%           To find all NaNs, Inf and -Inf and set them to zero you may use the following line:
%                       data_out(~isfinite(data_out)) = 0;
%
%EXAMPLES:
%X = imread('tire.tif'); imagesc(X); colormap gray; title('Original image')
%%This image can be found at \MATLAB704\toolbox\images\imdemos
%
%%Try one of the following:
%[dd, kernel] = sdgd_filt2D(X, 5);  dd = dd ./(max(dd(:)));
%figure; imagesc(dd> 0.1); colormap gray; colorbar
%[dd, kernel] = sdgd_filt2D(X, 11); dd = dd ./(max(dd(:)));
%figure; imagesc(dd> 0.1); colormap gray; colorbar
%
%%OR
%[dd, kernel] = sdgd_filt2D(X, 15, 0, 1, 2.5); 
%dd = dd ./(max(dd(:))); 
%figure; subplot(221); imagesc(dd> 0.1); %colormap gray
%title('Edge detected and thresholded image')
%subplot(222); surfl(kernel.dGx); colormap gray; shading interp
%title('kernel dGx'); axis tight
%subplot(223); surfl(kernel.dGxx); colormap gray; shading interp
%title('kernel dGxx'); axis tight
%subplot(224); surfl(kernel.dGxy); colormap gray; shading interp
%title('kernel dGxy'); axis tight
%
%REF: 
%1   Marr D., Hildreth E.C. in  Theory of edge Detection In Proceedings 
%    of the Royal Society, 207, pp. 187-217, 1980.
%2   Verbeek P.W., Vliet L.J., On the Location Error of Curved Edges
%    in Low-pass Filtered 2-D and 3-D Images, IEEE Trans. on Pattern 
%    Analysis and Machine Intelligence, v.16, n.7, pp.726-733, 1994
%
% Other m-files required: GaussDx2
% Subfunctions: none
% See also: sdgd_filt1D, plus_filt1D.m, plus_filt2D.m

%---------------------------------------------------|
% 	Sergei Koptenko, Resonant Medical, Montreal, Qc.|
% ph: 514.985.2442 ext265, www.resonantmedical.com  |
%      sergei.koptenko{at}resonantmedical.com       |
%Started--------Aug/07/2004-------------------------|
%Updated--------Nov/27/2005-------------------------|
%Updated--------Jun/07/2006-------------------------|

%% CHECK FOR MISSING VARIABLES
if nargin == 4, sigma_width = 3; 
    elseif nargin == 3 , sigma_width = 3; sigma = 1;
        elseif nargin == 2, sigma_width = 3; sigma = 1; mu = 0;
            elseif nargin <2, disp('function needs at least two variables!'); return;
end

warning off MATLAB:divideByZero % some image's pixels may have a_x=0 and a_y=0
data_in = double(data_in); % to avoid warning "CONV2 on values of class UINT8 is obsolete."
%% DERIVARIVES KERNELS
dGx = GaussDx2(x, mu, sigma, sigma_width); % First derivative dGx
dGy = dGx';                             % First derivative dGy
dGxx = conv2(dGx, dGx, 'same');         % Second Derivative dGxx
dGyy = dGxx';                           % Second Derivative dGyy
dGxy = conv2(dGx, dGy, 'same');           % Second Derivative dGxy 
%% STORE KERNELS
kernel.dGx = dGx;
kernel.dGxx = dGxx;
kernel.dGxy = dGxy;

%% DERIVATIVES
a_x = conv2(data_in, dGx, 'same');
a_y = conv2(data_in, dGy, 'same'); 
a_xx = conv2(data_in, dGxx, 'same');
a_yy = conv2(data_in, dGyy, 'same');
a_xy = conv2(data_in, dGxy, 'same');
%% OUTPUT OF THE FILTER
data_out = (a_xx .* a_x.^2 + 2* a_xy .* a_x .* a_y + a_yy .* a_y .^2 )./ (a_x .^2 + a_y .^2);
warning on MATLAB:divideByZero