|
|
convn and (i)fftn
Posted:
Jun 22, 2011 5:45 AM
|
|
Dear all,
I'm trying to speed up a convolution of a 3D image with a Gaussian kernel by using fftn() and ifftn() as convn() appears to be unacceptably slow...
It sort of works, but I cannot get my new code to exactly reproduce the results of the current code that uses convn.
I have distilled the differences down to the test program below -- my gut feeling is that only the kernel needs adjusting and then the fftn convolution should produce the same result as the convn version -- and quite a bit faster!
The program does not need toolboxes. The matlab command "[a,b]=gauss_filter_test;" should give the convn-result in variable a and the fftn result in variable b.
Does anyone see what the problem is?
Am I mixing up shifts, conjugates?
Does it need twice the #samples?
Many thanks
Alle Meije
<code>
function [ filtered, filtered2 ] = gauss_filter_test;
% 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 ( 44, 44, 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;
% kernel size
ksize = 3;
[ xg, yg, zg ] = ndgrid ( -ksize : ksize, -ksize : ksize, -ksize : ksize ) ;
% making the gaussian kernel
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 ( : ) ) ;
% making a kernel for fft convolution
dist2 = sqrt ( ( xg + ksize ) .^ 2 + ( yg + ksize ) .^ 2 + ( zg + ksize ) .^ 2 ) ;
gauss2 = 1 / ( sigma * sqrt ( 2 * pi ) ) * exp ( -0.5 * dist2 .^ 2 / ( sigma ) ^2 ) ;
gauss2 = gauss2 / sum ( gauss2 ( : ) ) ;
% 3-dimensional filtering with convn
filtered = convn ( volume, gauss, 'same' ) ;
% 3-dimensional filtering with fftn
gf = fftn ( gauss2, size ( volume ) ) ;
gv = fftn ( volume, size ( volume ) ) ;
filtered2 = real ( ifftn ( conj ( gf ) .* gv ) ) ;
% show the differences
close all;
for i=1:size(volume,3)
imagesc([filtered(:,:,i) filtered2(:,:,i)])
pause(.1);
end
</code>
|
|
