function G = GaussianR(I,sigma,padding)
%GAUSSIANR performs recursive gaussian filter
%   GaussianR takes 1d vector o 2d matrix and apply the Recursive Gaussian
%   Filter
%
%   G = GaussianR(I,sigma,padding)
%
%         I: Input 1d or 2d vector
%     sigma: Sigma for desidered gaussian (must be >=0.5)
%   padding: Desidered padding on input Specify is use padding or not. 
%            Choose 0 = no padding, 3 default padding or your desidered
%   
%
%   *See*:
%   # PkLab, _Implementazione di Gaussian smoothing ricorsivo_
%   <http://www.pklab.net/index.php?id=34>
%   # Young I.T. and L.J. Van Vliet, _Recursive Implementation of the Gaussian
%   Filter._ Signal Processing, 1995. 44(2): p. 139-151
%   <http://www.ph.tn.tudelft.nl/~lucas/publications/1995/SP95TYLV/SP95TYLV.html>
%   # Young I.T., Gerbrands J.J and L.J. Van Vliet, _Image Processing
%   Fundamentals 2.3_: p 57
%
% *Copyright (c) 2010 PkLab.net <http://www.pklab.net>*
% $Revision: 1.0 $  $Date: 2010/01/7 $

% Determine parameter q
if (sigma >= 2.5)
  q = 0.98711 * sigma - 0.96330;
elseif (sigma >= 0.5)
  q = 3.97156 - 4.14554 * sqrt(1 - 0.26891 * sigma);
else
  error('Gaussian sigma must be greater or equal than 0.5');
end;

% Determine coefficients
b0 = 1.57825 + (2.44413 * q) + (1.4281 * q^2) + (0.422205 * q^3);  
b1 = (2.44413 * q)  + (2.85619 * q^2) + (1.26661 * q^3);
b2 = -(1.4281 * q^2) - (1.26661 * q^3);
b3 = 0.422205 * q^3;
B  = 1 - ( ( b1 + b2 + b3 )/ b0 );

% Downscaling coefficient for efficient computations
B1 = b1/b0;
B2 = b2/b0;
B3 = b3/b0;

% Add border to reduce border error propagation
if(padding>0)
  border = 3;
  if(isvector(I))
    G = padarray(I,[0 border],'replicate','both');
  else
    G = padarray(I,[border border],'replicate','both');
  end
else
  border=0;
end
[R C] = size(G);

% Because gaussian is separable apply same 
% algorithm for rows and cols  

%---------------------------
% Apply algorithm along cols 

% Standard forward difference equation on COLS
% w[n] = B*G[n] + (b1*w[n-1] + b2*w[n-2] + b3*w[n-3])/b0;

% premultiply forward image with B coefficient can 
% take advance from matlab matrix computations
W = B * G;

% compute forward image using scaled 
% coefficient. This avoid NumCols division
% W[n] = W[n] + B1*w[n-1] + B2*w[n-2] + B3*w[n-3];

for r = 1:R
  W(r,2) = W(r,2) + B1 * W(r,1);
  W(r,3) = W(r,3) + B1 * W(r,2) + B2 * W(r,1);
  for c = 4:C
    W(r,c) = W(r,c) + ...
      B1 * W(r,c-1) + B2 * W(r,c-2) + B3 * W(r,c-3);
  end
end

% Backward difference equation on COLS
G = B * W;
for r = 1:R
  % Backward equation is applied from last col down to first col
  G(r,C-1) = G(r,C-1) + B1 * G(r,C);
  G(r,C-2) = G(r,C-2) + B1 * G(r,C-1) + B2 * G(r,C);
  for c = C-3:-1:1
    G(r,c) = G(r,c) + ...
      B1 * G(r,c+1) + B2 * G(r,c+2) + B3 * G(r,c+3);
  end;
end;
%
% At this point G contains gaussian filtered image along cols
%

%---------------------------
% Apply algorithm along rows
if(~isvector(I))
  W = B * G;
  % Forward difference equation on ROWS
  for c=1:C
    W(2,c) = W(2,c) + B1 * W(1,c);
    W(3,c) = W(3,c) + B1 * W(2,c) + B2 * W(1,c);
    for r=4:R
      W(r,c) = W(r,c) + ...
        B1 * W(r-1,c) + B2 * W(r-2,c) + B3 * W(r-3,c);
    end
  end

  % Backward difference equation on ROWS
  G = B * W;
  for c = 1:C
    % Backward equation is applied from last row down to first row
    G(R-1,c) = G(R-1,c) + B1 * G(R,c);
    G(R-2,c) = G(R-2,c) + B1 * G(R-1,c) + B2 * G(R,c);
    for r = R-3:-1:1
      G(r,c) = G(r,c) + ...
        B1 * G(r+1,c) + B2 * G(r+2,c) + B3 * G(r+3,c);
    end;
  end;
end;

% remove padding
if (border>0)
  if(isvector(I))
    G = G (border + 1:C-border);
  else
    G = G (border + 1:R-border, border + 1:C-border);
  end
end