MATLAB Language Fourier Transforms and Inverse Fourier Transforms Images and multidimensional FTs

Example

In medical imaging, spectroscopy, image processing, cryptography and other areas of science and engineering it is often the case that one wishes to compute multidimensional Fourier transforms of images. This is quite straightforward in Matlab: (multidimensional) images are just n-dimensional matrices, after all, and Fourier transforms are linear operators: one just iteratively Fourier transforms along other dimensions. Matlab provides fft2 and ifft2 to do this in 2-d, or fftn in n-dimensions.

One potential pitfall is that the Fourier transform of images are usually shown "centric ordered", i.e. with the origin of k-space in the middle of the picture. Matlab provides the fftshift command to swap the location of the DC components of the Fourier transform appropriately. This ordering notation makes it substantially easier to perform common image processing techniques, one of which is illustrated below.

Zero filling

One "quick and dirty" way to interpolate a small image to a larger size is to Fourier transform it, pad the Fourier transform with zeros, and then take the inverse transform. This effectively interpolates between each pixel with a sinc shaped basis function, and is commonly used to up-scale low resolution medical imaging data. Let's start by loading a built-in image example

%Load example image 
I=imread('coins.png'); %Load example data -- coins.png is builtin to Matlab 
I=double(I); %Convert to double precision -- imread returns integers 
imageSize = size(I); % I is a 246 x 300 2D image

%Display it
imagesc(I); colormap gray; axis equal; 
%imagesc displays images scaled to maximum intensity

We can now obtain the Fourier transform of I. To illustrate what fftshift does, let's compare the two methods:

% Fourier transform 
%Obtain the centric- and non-centric ordered Fourier transform of I
k=fftshift(fft2(fftshift(I))); 
kwrong=fft2(I); 

%Just for the sake of comparison, show the magnitude of both transforms: 
figure; subplot(2,1,1); 
imagesc(abs(k),[0 1e4]); colormap gray; axis equal; 
subplot(2,1,2); 
imagesc(abs(kwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear). 

We now have obtained the 2D FT of an example image. To zero-fill it, we want to take each k-space, pad the edges with zeros, and then take the back transform:

%Zero fill
kzf = zeros(imageSize .* 2); %Generate a  492x600 empty array to put the result in 
kzf(end/4:3*end/4-1,end/4:3*end/4-1) = k; %Put k in the middle 
kzfwrong = zeros(imageSize .* 2); %Generate a  492x600 empty array to put the result in 
kzfwrong(end/4:3*end/4-1,end/4:3*end/4-1) = kwrong; %Put k in the middle 

%Show the differences again 
%Just for the sake of comparison, show the magnitude of both transforms: 
figure; subplot(2,1,1); 
imagesc(abs(kzf),[0 1e4]); colormap gray; axis equal; 
subplot(2,1,2); 
imagesc(abs(kzfwrong),[0 1e4]); colormap gray; axis equal;
%(The second argument to imagesc sets the colour axis to make the difference clear). 

At this point, the result fairly unremarkable:

Once we then take the back-transforms, we can see that (correctly!) zero-filling data provides a sensible method for interpolation:

% Take the back transform and view 
Izf = fftshift(ifft2(ifftshift(kzf))); 
Izfwrong = ifft2(kzfwrong);

figure; subplot(1,3,1); 
imagesc(abs(Izf)); colormap gray; axis equal; 
title('Zero-filled image'); 
subplot(1,3,2); 
imagesc(abs(Izfwrong)); colormap gray; axis equal;
title('Incorrectly zero-filled image'); 
subplot(1,3,3); 
imagesc(I); colormap gray; axis equal; 
title('Original image'); 
set(gcf,'color','w'); 

Note that the zero-filled image size is double that of the original. One can zero fill by more than a factor of two in each dimension, although obviously doing so does not arbitrarily increase the size of an image.

Hints, tips, 3D and beyond

The above example holds for 3D images (as are often generated by medical imaging techniques or confocal microscopy, for example), but require fft2 to be replaced by fftn(I, 3), for example. Due to the somewhat cumbersome nature of writing fftshift(fft(fftshift(... several times, it is quite common to define functions such as fft2c locally to provide easier syntax locally -- such as:

function y = fft2c(x)

y = fftshift(fft2(fftshift(x)));

Note that the FFT is fast, but large, multidimensional Fourier transforms will still take time on a modern computer. It is additionally inherently complex: the magnitude of k-space was shown above, but the phase is absolutely vital; translations in the image domain are equivalent to a phase ramp in the Fourier domain. There are several far more complex operations that one may wish to do in the Fourier domain, such as filtering high or low spatial frequencies (by multiplying it with a filter), or masking out discrete points corresponding to noise. There is correspondingly a large quantity of community generated code for handling common Fourier operations available on Matlab's main community repository site, the File Exchange.