function [res, G_sp] = bilateral(M, R, sigma_col, sigma_sp)
%
% function res = bilateral(M, R, sigma_col, sigma_sp)
%
% Bilateral filtering 
%
% M: input image
% R: radius of the filter; 
%         the size of the filter is actually [2*R-1 x 2*R-1]
% 
% sigma_col: sigma of the filter in the color space
% sigma_sp: sigma of the filter in the spatial domain

% Direct implementation of the convolution
%      g(x,y) = 1/K(x,y) * sum_{n \in N(x,y)}
%      M(n)*(G(n;(x,y),sigma_sp)+G(M(n);M(x,y),sigma_col))
% where G(t;m,s) is the Gaussian having mean m and sigma s,
% and 1/K(x,y) is the normalizing factor, i.e., 
%        K(x,y) = sum_{n \in N(x,y)}(G(n;(x,y),sigma_sp)
%                                    +G(M(n);M(x,y),sigma_col))
% and N(x,y) is the neighborhood of (x,y)
%

res = zeros(size(M));

% zero padding
M1 = [zeros(R-1, size(M,2)+2*(R-1));...
      zeros(size(M,1), R-1) M zeros(size(M,1), R-1);...
      zeros(R-1, size(M,2)+2*(R-1))];

M1 = double(M1);
M = double(M);

t = ((-(R-1):R-1)'*ones(1,2*R-1)).^2;
sp_dist = sqrt(t+t');
G_sp = exp(-sp_dist.^2/sigma_sp^2);
for i=1:size(M,1), % rows
    i1 = i+R-1;
    for j=1:size(M,2), % columns
        j1 = j+R-1;
        col_dist = M1(i1-(R-1):i1+R-1, j1-(R-1):j1+R-1)...
            - M(i,j)*ones(2*R-1);
        G_col = exp(-col_dist.^2/sigma_col^2);
        t = M1(i1-(R-1):i1+R-1, j1-(R-1):j1+R-1).*G_col.*G_sp;
        norm_fact = sum(G_col(:).*G_sp(:));
        res(i,j) = sum(t(:))/norm_fact;
    end;
end;
