Rank: 2292 based on 54 downloads (last 30 days) and 2 files submitted
photo

Igor Solovey

E-mail

Personal Profile:

 

Watch this Author's files

 

Files Posted by Igor Solovey View all
Updated   File Tags Downloads
(last 30 days)
Comments Rating
07 Apr 2011 Screenshot Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: Igor Solovey bilateral, 3d, volume, fast bilateral filter, edgepreserving smooth... 41 7
  • 5.0
5.0 | 2 ratings
26 Feb 2011 Separable N-Dimensional Convolution N-dimensional convolution for separable kernels, similar to functionality of "conv2(hcol, hrow, A)" Author: Igor Solovey convolution, separability, nd, 3d, 4d, signal processing 13 5
  • 4.0
4.0 | 1 rating
Comments and Ratings by Igor Solovey View all
Updated File Comments Rating
28 Nov 2011 Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: 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.

28 Nov 2011 Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: Igor Solovey

Anton: it would be helpful if you provided a bit more information, e.g. what specific step caused the memory issues (line number would suffice), and more importantly what your results were when you tried a higher downsampling rate and a smaller subset of the image, i.e. 512x512x34. Without some effort on your part, your one-star rating is a little impolite and not helpful to others. Bilateral filtering is a resource-intensive algorithm, so my implementation is of course not infinitely scalable. I understand your frustration at having to deal with a problem, and I'd understand the one-star rating if the code didn't do what it advertises, or crashed on a smaller dataset.

29 May 2011 Separable N-Dimensional Convolution N-dimensional convolution for separable kernels, similar to functionality of "conv2(hcol, hrow, A)" Author: Igor Solovey

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

26 Feb 2011 Separable N-Dimensional Convolution N-dimensional convolution for separable kernels, similar to functionality of "conv2(hcol, hrow, A)" Author: Igor Solovey

Thank you, Maxime. Sorry for not including n2s - it's just an alias function I made for num2str. I submitted an update which simply replaces n2s with num2str.

As for trimarray, perhaps the way I wrote it is a little weird, but the goal is to trim pixels away from the edges of the result, to replicate the behaviour of convn with shape set to 'valid' or 'same'.

Comments and Ratings on Igor Solovey's Files View all
Updated File Comment by Comments Rating
29 Nov 2011 Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: Igor Solovey Anton Semechko

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

28 Nov 2011 Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: Igor Solovey 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.

28 Nov 2011 Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: Igor Solovey Anton Semechko

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

28 Nov 2011 Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: Igor Solovey Anton Semechko

Hi Igor,
I believe that my rating is entirely reflective of the quality of the code and its functionality. Your implementation of the 3D bilateral filter is simply not suitable for processing images of moderate complexity, even on a machine with ample memory resources.
In my previous comment, I pointed out that your function required too much memory to process a 1024x1024x33 image. I admit that images of this size are by no means small, so I tried your function on 256x256x33 image. It still run into exactly the same problem. Here is the error:
??? Error using ==> zeros
Out of memory. Type HELP MEMORY for your options.
Error in ==> bilateral3>bilateral3i at 148
GData =zeros(Xbins,Ybins,Zbins,1/samR);
Error in ==> bilateral3 at 81
Ibf=bilateral3i(I,sigmaSx,sigmaSy,sigmaSz,sigmaR,samSx,samSy,samSz,samR);

Determined to find the appropriate size of an image where your function would actually work, I downsampled the image once again to 128x128x33. This time, however, there was no error, but it gobbled up so much memory that my machine froze.

Below is my own, somewhat brute implementation of a 3D bilateral filter. Let me know how your advanced implementation compares to mine on the MRI image of the head included in the IP Toolbox.

% Compare performance (if you can get your function to run)
%-------------
load mri
D=squeeze(D);

t0=clock;
Da=Bilateral3D(D,[1 1 1],[3 5 6]);
t=etime(clock,t0);

t0=clock;
Db= bilateral3(D,3,3,5,1,1);
Bilateral3D(D,[1 1 1],[3 5 6]);
t(2)=etime(clock,t0);

28 Nov 2011 Fast Bilateral Filter for 3D Images Fast bilateral filter for 3D images Author: Igor Solovey Anton Semechko

.

Contact us