
%% Apply Gaussian Filter using Deriche Recursive Filter Approximation
%
% Given a sigma (standard deviation) parameter, constructs an approximation
% to the requested gaussian filter (ie. simple, derivative, second
% derivative).
% 
% Chooses an approximate impulse response ha[n] such that:
%    ha[k] = sum_{i=1}^n a[i]*exp(-l[i]*k/sigma)  ; k>=0
% which is minimizes (nonlinear LM optimization) the squared error between 
% the approximation and the actual gaussian filter:
%    min sum_{i=0}^{10*sigma} ( h[i] - ha[i] )^2
%
% Uses fourth order (n=4) filter ha[k]
%
%References:
%1] Deriche, Rachid. "Recursively implementating the Gaussian and 
%   its derivatives." (1993): 24.


%% Input Arguments

% img - Input grayscale or color image.
%       If color, applies gaussian filter on each channel seperately.
%       NxMx(3) matrix.

% sigma - Standard deviation gaussian filter.
%         Scalar (>0).

% type - Optional type of gaussian filter.
%        String {default: 'simple','first','second'}


%% Output Arguments

% I - Filtered image.
%     NxMx(3) matrix.


%% Approximate Gaussian Filter

function [Ix, Iy] = deriche_gaussian_filter(img, sigma, type)

%general constants
DEBUG = false;

%parse input parameters
[img, sigma, type, isColor] = parse_inputs(img, sigma, type);


%--------------------------------------------------
%		 <step 1>
%--------------------------------------------------
%description: generate IIR filter coefficients

%use pre-computed parameter values from [1] for 4th order filter
if strcmp(type,'simple')
    %a0,a1,c0,c1,w0,w1,b0,b1
    param = [1.68, 3.735, -0.6803, -0.2598, 0.6318, 1.997, 1.783, 1.723];
    symmetry = 'symmetric';
elseif strcmp(type,'first')
    param = [-0.6472, -4.531, 0.6494, 0.9557, 0.6719, 2.072, 1.527, 1.516];
    symmetry = 'antisymmetric';
else % second derivative
    param = [-1.331, 3.661, 0.3225, -1.738, 0.748, 2.166, 1.24, 1.314];
    symmetry = 'symmetric';
end
causality = 'noncausal';


%compute coefficients of IIR filter
[cNum, cDem, ncNum, ncDem, ~] = deriche_compute_coefs(param, sigma, ...
                                                      'causality',causality, ...
                                                      'symmetry',symmetry);

                                                  
%--------------------------------------------------
%		 <step 2>
%--------------------------------------------------
%description: apply IIR filter 

%apply padding at borders to prevent image-border effects
padSize = 4*sigma;
img = padarray(img, [padSize, padSize], 'replicate', 'both');

%apply filter     
if strcmp(type,'first')
    
    if isColor
                        
        Iy = apply_color_filter(img, cNum, cDem, 1);
        Ix = apply_color_filter(img, cNum, cDem, 2);
        
        if strcmp(causality,'noncausal')
            
            Iy = Iy + apply_noncausal_color_filter(img, ncNum, ncDem, 1);
            Ix = Ix + apply_noncausal_color_filter(img, ncNum, ncDem, 2);           
            
        end
        
    else
        
        Iy = apply_filter(img, cNum, cDem, 1);
        Ix = apply_filter(img, cNum, cDem, 2);
        
        if strcmp(causality,'noncausal')
           
            Iy = Iy + apply_noncausal_filter(img, ncNum, ncDem, 1);
            Ix = Ix + apply_noncausal_filter(img, ncNum, ncDem, 2);           
            
        end
        
    end 
    
else
    
    if isColor
                        
        I1 = apply_color_filter(img, cNum, cDem, 1);
        
        if strcmp(causality,'noncausal')
            
            I1 = I1 + apply_noncausal_color_filter(img, ncNum, ncDem, 1);
            
        end
                
        I2 = apply_color_filter(I1, cNum, cDem, 2);
        
        if strcmp(causality,'noncausal')
            
            I2 = I2 + apply_noncausal_color_filter(I1, ncNum, ncDem, 2);           
            
        end
        
    else
        
        I1 = apply_filter(img, cNum, cDem, 1);
        
        if strcmp(causality,'noncausal')
            
            I1 = I1 + apply_noncausal_filter(img, ncNum, ncDem, 1);
            
        end
                
        I2 = apply_filter(I1, cNum, cDem, 2);
        
        if strcmp(causality,'noncausal')
            
            I2 = I2 + apply_noncausal_filter(I1, ncNum, ncDem, 2);           
            
        end
        
    end 
    
    Ix = I2;
    Iy = Ix;
   
end

%crop padding from image
Ix = Ix(padSize+1:end-padSize, padSize+1:end-padSize, :);
Iy = Iy(padSize+1:end-padSize, padSize+1:end-padSize, :);


%--------------------------------------------------
%		 <step 3>
%--------------------------------------------------
%description: show debug information

%display impulse response
if DEBUG
    N = 100;
    figure(1);
    if strcmp(causality,'noncausal')
        [cH,~]  = impz(cNum,[1,cDem],N);
        [ncH,~] = impz([0,ncNum],[1,ncDem],N);
        h = [flipud(ncH); zeros(N-1,1)] + [zeros(N-1,1); cH];
        stem(-N+1:N-1,h);
    else
        [h,t] = impz(cNum,[1,cDem]);
        stem(t, h);
    end
    title('Impulse Response');
end

    
end %end of deriche gaussian filter


%% Helper Functions

function outImg = apply_color_filter(inImg, n, d, dim)
%applies the filter:
%
% y[n] = n0*x[n] + n1*x[n-1] + ... + nN*x[n-N] 
%                - d1*y[n-1] - ... - dN*y[n-N]

outImg = inImg;
for i = 1:3,
    outImg(:,:,i) = filter(n, [1,d], inImg(:,:,i), [], dim);
end

end %end of apply color filter

function outImg = apply_noncausal_color_filter(inImg, n, d, dim)

if dim == 1
    outImg = flipud( apply_color_filter(flipud(inImg), [0,n], d, dim) );  
else
    outImg = fliplr( apply_color_filter(fliplr(inImg), [0,n], d, dim) );
end

end% end of apply non-causal filter

function outImg = apply_filter(inImg, n, d, dim)
%applies the filter:
%
% y[n] = n0*x[n] + n1*x[n-1] + ... + nN*x[n-N] 
%                - d1*y[n-1] - ... - dN*y[n-N]

outImg = filter(n, [1,d], inImg, [], dim);

end %end of apply filter

function outImg = apply_noncausal_filter(inImg, n, d, dim)

if dim == 1
    outImg = flipud( apply_filter(flipud(inImg), [0,n], d, dim) );  
else
    outImg = fliplr( apply_filter(fliplr(inImg), [0,n], d, dim) );
end

end% end of apply non-causal filter


%% Input Parser

function [img, sigma, type, isColor] = parse_inputs(img, sigma, type)

% parse input arguments
p = inputParser;
p.addRequired('img', @(x) ...
				validateattributes(x,{'numeric'},{'real','nonnan'})); 
p.addRequired('sigma',  @(x) ...
				validateattributes(x,{'numeric'},{'real','scalar','positive'}));   
p.addRequired('type', @(x) ischar(x) && (strcmp(x,'simple') || strcmp(x,'first') || strcmp(x,'second')));			       
p.parse(img, sigma, type);
inputs = p.Results;
img    = inputs.img;
sigma  = inputs.sigma;
type   = inputs.type;

nComps = size(img,3);

if nComps == 3
    isColor = true;
elseif nComps == 1
    isColor = false;
else
    error('dericheGuassFilter:invalidImg','Invalid image: should be grayscale or color image');
end

img = double(img);

end %end of input parser








