Search All of the Math Forum:

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

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

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

 Messages: [ Previous | Next ]
 Alle Mejie Posts: 12 Registered: 2/24/11
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>

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