Igor Solovey

This function is an implementation of the ">>fast<< bilateral filter" - don't think of fast as some kind of a boast, but as a name for an approximation/acceleration technique. If you read the paper I referenced (Paris/Durand), you will see that a principal claim is that you can subsample the resulting 4D volume with a relatively small loss of accuracy.

My code presumes you are doing that (without any checks, which is a legitimate criticism), and hence if you don't subsample, you'll be allocating two N*M*K*256-sized arrays, which is bad news.

It is impossible to fairly compare your implementation with mine since you do not downsample your data, which is an approximation step. So my code is painfully slow without subsampling, and will generate an approximation when subsampling is used. Moreover, I think your Gaussian filter gets cut off at 1 stdev, while I extend mine to 3.5 (probably an overkill). You can, however, compare the result of running bilateral3 with samSxy=samSz=samR=1 and, for example, samSxy=samSz=samR=3 or 4.

The claim of "fast bilateral filter" technique is that the error (as compared to a true implementation, i.e. yours) will be small, but the speedup will be considerable.

Long story short: go ahead and try " tic,Db= bilateral3(D,3,3,5,2,2);toc ", " tic,Db2= bilateral3(D,3,3,5,4,4);toc " and compare the results/runtimes. It should be relatively easy to modify gkernel function in bilateral3.m to match your implementation of a Gaussian filter if you wish to compare directly to the output of your Bilateral3D.

Thanks Taemin, you're right. I wasn't aware that you can just supply a 1D vector as one of the arguments of convn.

%normalized gaussian kernel
a=[.0003 .1065 .7866 .1065 .0003];
n=length(a);
%make a test matrix
A = randn(64,64,64,64);
tic
b=a;
B3=convn(b,A);
b=a';
B3=convn(b,B3);
b=reshape(a,[1 1 n]);
B3=convn(b,B3);
b=reshape(a,[1 1 1 n]);
B3=convn(b,B3);
toc

Thanks for clearing that up Igor. Five stars for the effort!

And here my version of a 3D bilateral filter:

function IM_out=Bilateral3D(IM,VoxRes,opt)
% Preprocess a grayscale 3D image using bilateral filter.
%
% INPUT ARGUMENTS:
% - IM : M-by-N-by-D grayscale image. Note that during filtering IM
% is cast into a double format.
% - VoxRes : VoxRes=[dx dy dz], where dx, dy and dz are voxel dimensions
% along the x, y and z directions, respectively. The default
% setting is VoxRes=[1 1 1].
% - opt : opt=[Wd Wr R n ref pad]
% - Wd : standard deviation of the spatial-domain filter. Wd
% must be a positive integer. For instance, Wd=n
% actually corresponds to the standard deviation of
% n*dx units. Wd=3 is default.
% - Wr : standard deviation of the range (i.e. intensity)
% filter. Assuming IM is represented in an 8 bit
% format, the default setting for Wr is 5.
% - R : net filter radius. R must be a positive integer
% greater than or equal to 1. For instance, specifying
% R=n is equivalent to using filter window that has
% the radius of n*dx units. R=6 is default.
% - pad : padding option:
% pad=0: pad with zeros (default)
% pad=1: symmetric
% pad=2: replicate
% - n : number of filter iterations. n=1 is default.
% - ref : This parameter is used only when n>1. Set ref=1 to
% compute range filter responses (RFRs) with respect
% to intensities of the original input image (default).
% Set ref=2 to compute RFRs with respect to the
% intensities of the most recent iterate.
%
% OUTPUT:
% - IM_out : output image (in double format).
%
% REFERENCES:
% [1] Tomasi C., Manduchi, R. (1998) 'Bilateral filtering for gray and
% color images'
%
% AUTHOR : Anton Semechko (a.semechko@gmail.com)
% DATE : Nov.2011
%

% Verify format of the input agruments
if nargin<2 || isempty(VoxRes), VoxRes=[1 1 1]; end
if nargin<3 || isempty(opt), opt=[3 5 6 1 1 0]; end
opt=CheckInputArgs(IM,VoxRes,opt);

% Filter options
global Wr R pad h ix iy iz idx IM_ref
Wd=opt(1); Wr=opt(2); R=opt(3); pad=opt(4); n=opt(5); ref=opt(6);
if pad==1
pad='symmetric';
elseif pad==2
pad='replicate';
end

% Construct the spatial-domain filter
VoxRes=VoxRes/VoxRes(1); % normalize VoxRes wrt to x-dir
hx=0:VoxRes(1):R; hx=[fliplr(hx),hx(2:end)];
hy=0:VoxRes(2):R; hy=[fliplr(hy),hy(2:end)];
hz=0:VoxRes(3):R; hz=[fliplr(hz),hz(2:end)];
[hx,hy,hz]=meshgrid(hx,hy,hz);
D=hx.^2+hy.^2+hz.^2;
h=exp(-0.5*D/(Wd.^2));
h(D>R.^2)=0;

% Get the non-zero indices of h
idx=find(D<=R.^2);
[iy,ix,iz]=ind2sub(size(D),idx);
clear hx hy hz D

% Center the indices
c=floor(size(h)/2)+1;
iy=iy-c(1); ix=ix-c(2); iz=iz-c(3);

if ref==1
IM_ref=padarray(double(IM),[R R R],pad);
else
IM_ref=[];
end

% Apply the filter
for i=1:n, IM=CompFiltProd(IM); end
IM_out=IM;

%==========================================================================
function IM_out=CompFiltProd(IM)

global Wr R pad h ix iy iz idx IM_ref

if ~strcmpi(class(IM),'double'), IM=double(IM); end

% Pad the image
[M,N,D]=size(IM);
IM=padarray(IM,[R R R],pad);

% Progress bar
msg=sprintf('Filtering image %3.0f%%...',0);
h_bar = waitbar(0,msg);
set(h_bar,'Name','3D BILATERAL FILTER','Color','w');
drawnow

% Apply the filter
IM_out=zeros(size(IM));
Fnet=zeros(size(IM));
Nitr=numel(ix);
for i=1:Nitr

im=circshift(IM,[iy(i) ix(i) iz(i)]);

if ~isempty(IM_ref)
im_ref=circshift(IM_ref,[iy(i) ix(i) iz(i)]);
else
im_ref=im;
end

f=h(idx(i))*exp(-0.5*((im_ref-IM).^2)/(Wr^2));
Fnet=Fnet+f;
IM_out=IM_out+f.*im;

if ishandle(h_bar)
msg=sprintf('Filter in progress: %3.0f%% complete ...',i/Nitr*100);
waitbar(i/Nitr,h_bar,msg)
end

end
clear IM im f im_ref
IM_out=IM_out./Fnet;

% Remove the padding
IM_out=IM_out((R+1):(M+R),(R+1):(N+R),(R+1):(D+R));
if ishandle(h_bar), delete(h_bar); end

%==========================================================================
function opt=CheckInputArgs(IM,vr,opt)
% Verify the consistency of the input arguments

siz=size(IM);
if numel(siz)~=3 || siz(1)<siz(3) || siz(1)<siz(3) || ~isnumeric(IM)
error('Incorrect entry for 1st input argument. IM=3D grayscale image.')
end

if sum(isnan(IM(:)))>0 || sum(isinf(IM(:)))>0
error('Input image contains NaN and/or Inf elements. Remove them and try again.')
end

if ~isnumeric(vr) || ~isvector(vr) || numel(vr)~=3 || sum(vr<eps)>0 || sum(isnan(vr))>0 || sum(isinf(vr))>0
error('Incorrect entry for 2nd input argument. VoxRes=[dx dy dz], where dx, dy and dz are voxel dimensions.')
end

if ~isvector(opt) || numel(opt)<3
error('Incorrect entry for 3rd input argument. opt=[Wd Wr R pad n ref]. See function documentation for more info.')
end

Wd=opt(1); Wr=opt(2); R=opt(3);
if Wd<eps || isnan(Wd),error('Revise entry for Wd option. Wd>0'); end
if Wr<eps || isnan(Wr), error('Revise entry for Wr option. Wr>0'); end
if R<1 || isnan(R), error('Revise entry for R option. R>=1'); end

if numel(opt)<4
opt(4)=0;
else
pad=opt(4);
if ~(pad==0 || pad==1 || pad==2), error('Revise entry for pad option. See function documentation for more info.'); end
end

if numel(opt)<5
opt(5)=1;
else
n=opt(5);
if n<1 || round(n)~=n, error('Revise entry for n option. n must be a postive integer >=1'); end
end

if numel(opt)<6
opt(6)=1;
else
ref=opt(6);
if ~(ref==1 || ref==2), error('Revise entry for ref option. ref=1 or ref=2'); end
end

.