## Nice implementation of bilateral filtering

Nice implementation, e.g. the progress bar.

% BFILTER2 Two dimensional bilateral filtering.

% This function implements 2-D bilateral filtering using

% the method outlined in:

%

% C. Tomasi and R. Manduchi. Bilateral Filtering for

% Gray and Color Images. In Proceedings of the IEEE

% International Conference on Computer Vision, 1998.

%

% B = bfilter2(A,W,SIGMA) performs 2-D bilateral filtering

% for the grayscale or color image A. A should be a double

% precision matrix of size NxMx1 or NxMx3 (i.e., grayscale

% or color images, respectively) with normalized values in

% the closed interval [0,1]. The half-size of the Gaussian

% bilateral filter window is defined by W. The standard

% deviations of the bilateral filter are given by SIGMA,

% where the spatial-domain standard deviation is given by

% SIGMA(1) and the intensity-domain standard deviation is

% given by SIGMA(2).

%

% Douglas R. Lanman, Brown University, September 2006.

% dlanman@brown.edu, http://mesh.brown.edu/dlanman

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Pre-process input and select appropriate filter.

function B = bfilter2(A,w,sigma)

% Verify that the input image exists and is valid.

if ~exist(‘A’,’var’) || isempty(A)

error(‘Input image A is undefined or invalid.’);

end

if ~isfloat(A) || ~sum([1,3] == size(A,3)) || …

min(A(:)) < 0 || max(A(:)) > 1

error([‘Input image A must be a double precision ‘,…

‘matrix of size NxMx1 or NxMx3 on the closed ‘,…

‘interval [0,1].’]);

end

% Verify bilateral filter window size.

if ~exist(‘w’,’var’) || isempty(w) || …

numel(w) ~= 1 || w < 1

w = 5;

end

w = ceil(w);

% Verify bilateral filter standard deviations.

if ~exist(‘sigma’,’var’) || isempty(sigma) || …

numel(sigma) ~= 2 || sigma(1) <= 0 || sigma(2) <= 0

sigma = [3 0.1];

end

% Apply either grayscale or color bilateral filtering.

if size(A,3) == 1

B = bfltGray(A,w,sigma(1),sigma(2));

else

B = bfltColor(A,w,sigma(1),sigma(2));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Implements bilateral filtering for grayscale images.

function B = bfltGray(A,w,sigma_d,sigma_r)

% Pre-compute Gaussian distance weights.

[X,Y] = meshgrid(-w:w,-w:w);

G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));

% Create waitbar.

h = waitbar(0,’Applying bilateral filter…’);

set(h,’Name’,’Bilateral Filter Progress’);

% Apply bilateral filter.

dim = size(A);

B = zeros(dim);

for i = 1:dim(1)

for j = 1:dim(2)

% Extract local region.

iMin = max(i-w,1);

iMax = min(i+w,dim(1));

jMin = max(j-w,1);

jMax = min(j+w,dim(2));

I = A(iMin:iMax,jMin:jMax);

% Compute Gaussian intensity weights.

H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));

% Calculate bilateral filter response.

F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);

B(i,j) = sum(F(:).*I(:))/sum(F(:));

end

waitbar(i/dim(1));

end

% Close waitbar.

close(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Implements bilateral filter for color images.

function B = bfltColor(A,w,sigma_d,sigma_r)

% Convert input sRGB image to CIELab color space.

if exist(‘applycform’,’file’)

A = applycform(A,makecform(‘srgb2lab’));

else

A = colorspace(‘Lab<-RGB’,A);

end

% Pre-compute Gaussian domain weights.

[X,Y] = meshgrid(-w:w,-w:w);

G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));

% Rescale range variance (using maximum luminance).

sigma_r = 100*sigma_r;

% Create waitbar.

h = waitbar(0,’Applying bilateral filter…’);

set(h,’Name’,’Bilateral Filter Progress’);

% Apply bilateral filter.

dim = size(A);

B = zeros(dim);

for i = 1:dim(1)

for j = 1:dim(2)

% Extract local region.

iMin = max(i-w,1);

iMax = min(i+w,dim(1));

jMin = max(j-w,1);

jMax = min(j+w,dim(2));

I = A(iMin:iMax,jMin:jMax,:);

% Compute Gaussian range weights.

dL = I(:,:,1)-A(i,j,1);

da = I(:,:,2)-A(i,j,2);

db = I(:,:,3)-A(i,j,3);

H = exp(-(dL.^2+da.^2+db.^2)/(2*sigma_r^2));

% Calculate bilateral filter response.

F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);

norm_F = sum(F(:));

B(i,j,1) = sum(sum(F.*I(:,:,1)))/norm_F;

B(i,j,2) = sum(sum(F.*I(:,:,2)))/norm_F;

B(i,j,3) = sum(sum(F.*I(:,:,3)))/norm_F;

end

waitbar(i/dim(1));

end

% Convert filtered image back to sRGB color space.

if exist(‘applycform’,’file’)

B = applycform(B,makecform(‘lab2srgb’));

else

B = colorspace(‘RGB<-Lab’,B);

end

% Close waitbar.

close(h);