Code covered by the BSD License  

Highlights from
3D convolution in the FFT domain

image thumbnail

3D convolution in the FFT domain

by

 

13 Mar 2012 (Updated )

Achieve a 3D convolution in the fourrier domain.

convolution3D_FFTdomain(inVol,inKer)
function [outVol] = convolution3D_FFTdomain(inVol,inKer)
%% convolution3D_FFTdomain - Performs a fast 3D convolution between a volume and a kernel using mutliplication in the Fourrier space. 
%
% Syntax:  [outVol] = convolutionInFFTdomain(inVol,inKer,inMsg)
% Inputs:
%       inVol - input volume (real / complex)
%       inKer - input kernel (real / complex)
% Outputs:
%       outVol - output convolved volume (real / complex) - precision of the output format is the same as the input
%       volume. The output volume is the central part of the convolution with same size as inVol.
%       size(outVol)=size(inVol) ('same' option of convn).
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none

% Author: Christopher Coello
% Work address: Preclinical PET/CT Unit
% Email address: s.c.coello@medisin.uio.no
% Website: http://www.med.uio.no/imb/english/services/public/pet/
% 2012/03/13; Last revision: 2012/03/16
% Created with Matlab version: 7.13.0.564 (R2011b)
%%

% Check if both inputs are real numbers
realInput = isreal(inVol) && isreal(inKer);

% Check the precision of the input volume
if strcmp(class(inVol),'double'),
    indDb=1;
else
    indDb=0;
end

% Size of the input volumes
inVolSize=size(inVol);
inVolSide=max(inVolSize);
inKerSize=size(inKer);
inKerSide=max(size(inKer));

% Fourrier tranform of the volume and inKer.
extr(1:3)={};
for iDim=(1:3),
    inVol=fft(inVol,inVolSide+inKerSide-1,iDim);
    inKer=fft(inKer,inVolSide+inKerSide-1,iDim);
    extr{iDim}=ceil((inKerSize(iDim)-1)/2)+(1:inVolSize(iDim));
end

% Multiplication of the Fourrier tranforms
conv_FFT=inVol.*inKer;

% Inverse Fourrier Transform of the convolution
for iDim=(1:3),
    conv_FFT=ifft(conv_FFT,[],iDim);
end

% Crop the side of the image in relation to the size of the kernel
convinVol=conv_FFT(extr{:});

% limit the results
if realInput,
    convinVol=real(convinVol);
end

% Format the output into single precision if needed
if indDb,
    outVol=convinVol;
else
    outVol=single(convinVol);
end

Contact us