Thread Subject: Out of memory when I am doing 3D FFT

Hi, there

I am working on CT image processing. I need to calculate the Hessian matrix of the volumetric dataset. It can be accomplished by convolution the image with the 2nd order partial derivatives of the Gaussian kernel. Here is the way I am doing the convolution.


function FilterOut=perform_convolution_3d(Vol,vker)

%% Vol is the input volumetric image data whose size is 341 by 341 by 341, uint16,
%% vker is an isotropic Gaussian vker whose size is 30 by 30 by 30

SizeV=size(Vol);

smoothcell=zeros(SizeV);
centerpix=floor(SizeV/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;

FFT1=fftn(double(Vol));
FFT2=fftn(double(smoothcell));
CPXARR2=FFT1.*conj(FFT2);
MFout=fftshift(ifftn(CPXARR2));

But the 'out of memory' problem always occurs when I implement the code. I tried to save some of the un-used variables to the disk, only leaving the essential ones in the workspace, but the 'out of memory' problem still happens.

I use window xp 32 OP system with 2GB RAM.

Thanks for the help,
David

"David " <david2008.w@gmail.com> wrote in message <gk7vku$naj$1@fred.mathworks.com>...
> Hi, there
>
> I am working on CT image processing. I need to calculate the Hessian matrix of the volumetric dataset. It can be accomplished by convolution the image with the 2nd order partial derivatives of the Gaussian kernel. Here is the way I am doing the convolution.
>
>
> function FilterOut=perform_convolution_3d(Vol,vker)
>
> %% Vol is the input volumetric image data whose size is 341 by 341 by 341, uint16,
> %% vker is an isotropic Gaussian vker whose size is 30 by 30 by 30
>
> SizeV=size(Vol);
>
> smoothcell=zeros(SizeV);
> centerpix=floor(SizeV/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;
>
> FFT1=fftn(double(Vol));
> FFT2=fftn(double(smoothcell));
> CPXARR2=FFT1.*conj(FFT2);
> MFout=fftshift(ifftn(CPXARR2));
>
> But the 'out of memory' problem always occurs when I implement the code. I tried to save some of the un-used variables to the disk, only leaving the essential ones in the workspace, but the 'out of memory' problem still happens.
>
> I use window xp 32 OP system with 2GB RAM.

How big are your arrays? Since they seem to be pretty big, maybe you'd settle for representing this as singles rather than doubles.

Other than that, you probably don't need to compute FFT2, or at least not as a 3D array whose size is the same as Vol.

Since your kernel is Gaussian, it is separable, and so are its derivatives. So, instead, you should be computing 1D kernels and using 1D fft's to do the convolution on each row, column, etc... of Vol.

Follow how memory is used on your computer in windows task manager
(ctrl+alt+del). See MATLAB.exe. Add virtual memory colomn. Find info
about 3GB mode in Windows. In 64 bit computer there are much more
virtual memory is avaliable.

-----------------------------------------
Maxim Vedenev, MATLAB Custom Programming
http://simulations.narod.ru/

"David " <david2008.w@gmail.com> wrote in message <gk7vku$naj$1@fred.mathworks.com>...

>
> FFT1=fftn(double(Vol));
> FFT2=fftn(double(smoothcell));
> CPXARR2=FFT1.*conj(FFT2);
> MFout=fftshift(ifftn(CPXARR2));
>

Few tips:
- use SINGLE instead of DOUBLE
- Inplace fft by slice:

  FFT1=zeros(size(Vol),class(Vol));
  for k=1:size(Vol,3)
     FFT1(:,:,k)=fftn(Vol(:,:,k));
  end
  for i=1:size(Vol,1)
     for j=1:size(Vol,2)
        FFT1(i,j,:)=fft(squeeze(Vol(i,j,:)));
     end
  end
  
  % Same for FFT2

- Clear intermediate variables or reuse the same variables at different steps
- As someone has pointed out: because the kernel is tensorial of 1d kernels, Do the convolution in 1D separately.

Bruno

:,:,k)=fftn(Vol(:,:,k));
> end
> for i=1:size(Vol,1)
> for j=1:size(Vol,2)
> FFT1(i,j,:)=fft(squeeze(Vol(i,j,:)));
> end
> end

The computation statement should be read as
FFT1(i,j,:)=fft(squeeze(FFT1(i,j,:)));

Bruno

Hi,
thanks for help, I belive that perform the 3D convolution by using three separate one dimensional kernel is good choice. As Bruno and Matt suggested, the use of 1D kernel avoid to calculating the FFT of the 3D smooth kernel which reduces the memory cost.