Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Topic: convn and (i)fftn
Replies: 10   Last Post: May 16, 2016 4:51 PM

 Search Thread: Advanced Search

 Messages: [ Previous | Next ]
 Alle Mejie Posts: 12 Registered: 2/24/11
Re: convn and (i)fftn
Posted: Jun 23, 2011 5:49 AM
 Plain Text Reply

"Matt J" wrote in message <itsqn1\$ebj\$1@newscl01ah.mathworks.com>...
> There's no reason you should be using FFT's with such a small kernel. Also,
> you're not exploiting the separability of the Gaussian. The separability allows you to
> convolve with a 1D Gaussian first along x, then along y, then along z, which is more
> efficient than consolidating into a 3D convolution kernel.

Got it! I needed to download a substitute for bsx_func() which is in the latest matlab version only, and change a line where it says
[~, output] = function(inputs)
which is also new-matlab syntax?

Also had to apply a circular shift of [-ksize -ksize -ksize] to the volume (because the kernel was not centred), but still it was *A LOT* faster than the other two.

>> [a b c]=gauss_filter_test;
gave me
KronX - based convolution -> Elapsed time is 0.078794 seconds.
fftn() - based convolution -> Elapsed time is 0.298942 seconds.
convn() - based convolution -> Elapsed time is 6.081766 seconds.

>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>

function [ filtered, filtered2, filtered3 ] = gauss_filter_test;

addpath('kronProd') % http://www.mathworks.com/matlabcentral/fileexchange/25969
addpath('interpMatrix') % http://www.mathworks.com/matlabcentral/fileexchange/26292
addpath('bsxfun_substitute') % http://www.mathworks.com/matlabcentral/fileexchange/23005

% volume of [ 144 x 144 x 44 ] zeros
% with a white 'box' of thickness 5,
% at 10 points % away from the edge
%
% This idea :
% 0 0 0 0 0 0 0 0
% 0 1 1 1 1 1 1 0
% 0 1 0 0 0 0 1 0
% 0 1 0 0 0 0 1 0
% 0 1 0 0 0 0 1 0
% 0 1 0 0 0 0 1 0
% 0 1 1 1 1 1 1 0
% 0 0 0 0 0 0 0 0
% but then in 3D
%
edge = 10; width = 5; inside = edge + width;
volume = zeros ( 144, 144, 44 ) ;
volume ( edge : ( size ( volume, 1 ) -edge + 1 ) , ...
edge : ( size ( volume, 2 ) -edge + 1 ) , ...
edge : ( size ( volume, 3 ) -edge + 1 ) ) = 1;
volume ( inside : ( size ( volume, 1 ) -inside + 1 ) , ...
inside : ( size ( volume, 2 ) -inside + 1 ) , ...
inside : ( size ( volume, 3 ) -inside + 1 ) ) = 0;

% settings for PSF
pixelmmsize = 4; fwhm = 6.5; width = fwhm / pixelmmsize; sigma = width / 2.35482; ksize=3;

%% KronProd based convolution %%
fprintf('KronX - based convolution -> '); warning off; tic
%%
dist = sqrt ( (-ksize:ksize).^ 2) ;
gauss = 1 / ( sigma * sqrt ( 2 * pi ) ) * exp ( -0.5 * dist .^ 2 / ( sigma ) ^2 ) ;
gauss=gauss/sum(gauss(:));
%%
Hxy=interpMatrix(gauss,'max',144);
Hz=interpMatrix(gauss,'max',44);
K=KronProd({Hxy,Hz},[1 1 2]);
filtered3=K*circshift(volume,[-3 -3 -3]);
%%
toc; warning on

%% frequency-domain convolution %%
fprintf('fftn() - based convolution -> '); tic
%%
[ xg, yg, zg ] = ndgrid ( -ksize : ksize, -ksize : ksize, -ksize : ksize ) ;
dist = sqrt ( xg .^ 2 + yg .^ 2 + zg .^ 2 ) ;
gauss = 1 / ( sigma * sqrt ( 2 * pi ) ) * exp ( -0.5 * dist .^ 2 / ( sigma ) ^2 ) ;
gauss = gauss / sum ( gauss ( : ) ) ;
%%
filtered2 = conv3Dfreq ( volume, gauss ) ;
%%
toc

fprintf('convn() - based convolution -> '); tic
[ xg, yg, zg ] = ndgrid ( -ksize : ksize, -ksize : ksize, -ksize : ksize ) ;
dist = sqrt ( xg .^ 2 + yg .^ 2 + zg .^ 2 ) ;
gauss = 1 / ( sigma * sqrt ( 2 * pi ) ) * exp ( -0.5 * dist .^ 2 / ( sigma ) ^2 ) ;
gauss = gauss / sum ( gauss ( : ) ) ;
%%
filtered = convn ( volume, gauss, 'same' ) ;
%%
toc

% show the differences
close all;
for i=1:size(volume,3)
imagesc([filtered(:,:,i) filtered2(:,:,i) filtered3(:,:,i)]);
pause(.1);
end

rmpath('bsxfun_substitute')
rmpath('interpMatrix')
rmpath('kronProd')

return;

function MFout=conv3Dfreq(cprocl,vker)
% function from http://www.mathworks.com/matlabcentral/newsreader/view_thread/156782

cprocl=double(cprocl);
smoothcell=zeros(size(cprocl));
centerpix=floor(size(cprocl)/2);
centerker=floor(size(vker)/2);

embed1i=centerpix(1)-centerker(1)+[1:size(vker,1)];
embed2i=centerpix(2)-centerker(2)+[1:size(vker,2)];
embed3i=centerpix(3)-centerker(3)+[1:size(vker,3)];

smoothcell(embed1i,embed2i,embed3i)=vker;

% disp('matched filter correlating observed volume');
% Correlation through FFT needs:
% a. Finding FFT(S) and FFT(C)
FFT1=fftn(cprocl);
FFT2=fftn(smoothcell);
% b. ComplexArray=FFT(S) * Conj{FFT(C)}
CPXARR2=FFT1.*conj(FFT2);
% c. IFFT{ComplexArray}
MFout=abs(fftshift(ifftn(CPXARR2)));

return;

<<<<<<<<<<<
<<<<<<<<<<<
<<<<<<<<<<<

Date Subject Author
6/22/11 Alle Mejie
6/22/11 Alle Mejie
6/22/11 Matt J
6/22/11 Alle Mejie
6/22/11 Matt J
6/23/11 Alle Mejie
6/23/11 Alle Mejie
6/23/11 Alle Mejie
6/23/11 Matt J
6/23/11 Alle Mejie
5/16/16 mathango

© The Math Forum at NCTM 1994-2017. All Rights Reserved.