Thursday, 16 April 2009

Frequency Transforms on Images

As part of my work, I have to do a lot of spectral analysis with images. Traditionally I have used MATLAB to do this - it is exceptionally easy and powerful, and works just fine. As a simple example, if I want a power-spectrum of an image, treated as a 2D matrix of 8-bit numbers, the sequence of commands is this:

 %------------------- MATLAB  

% read image in, convert to greyscale and double-precision
% crop to smallest dimension
% find image mean
% and subtract mean from the image
% find the 2D fourier transform
% power spectrum is the element-by-element square of the fourier transform

Seven lines of commands. Couldn't be simpler. That's the beauty of MATLAB. The snag is that it costs several hundreds to thousands of dollars for the program, which makes it difficult to share your code with other people. So recently I've thought about coding up some things I want to share with other people in Python.

Python is a cross-platform, open-source programming language and environment. It is not too dissimilar to MATLAB in many respects, especially when you install the right libraries. In order to do the same in Python as above, as easy as possible, I needed to install Python Imaging Library (PIL),  NumPy and SciPy. Then the sequence of operations is as follows:

 # ----------------PYTHON  

# declare libraries to use
import Image
import scipy.fftpack
from pylab import*
import cmath
from numpy import*
# load image, convert to greyscale'myimage.jpg').convert("L")
# crop to smallest dimension
# find mean
# 2D fourier transform of demeaned image
# power spectrum

Same answer. Not bad - not harder, just different. And completely free. Oh, and this'll work on UNIX, Linux, Mac, and PC when Python is installed (with a couple of extra libraries). Great! 

Hmm, I thought, I wonder if it is possible to do the same in a piece of non-numerical software. I tried ImageMagick because its free and ubiquitous on Linux and Mac, and also available for PC, and it has a massive user-community. It is incredibly powerful as a command-line graphics package - even more powerful than PhotoShop and Illustrator, and like all UNIX-based applications, infinitely flexible. It turns out it is a little more tricky to calculate the power spectrum in ImageMagick, but here's how.

First, you need the right tools - ImageMagick is usually already present on most Linux distributions - if not, if a debian/Ubuntu user use sudo apt-get install imagemagick, or download from the webiste hyperlinked above. Then you'll need an open-source C-program for calculation of FFTs: this is provided by FTW - download fftw-3.2.1.tar.gz, untar, cd to that directory, then install the binary by issuing the usual commands:


sudo make && make install

This C-program is no good on its own, so the last ingredient is a wrapper for ImageMagick. Fortunately - nay, crucially, this has been provided by Sean Burke and Fred Weinhaus - introducing IMFFT! The program can be downloaded and installed from here. And here is an in-depth tutorial of how to use it (I'm only just getting to grips with it myself). For some reason the compiled program is called 'demo'. I move it to my /bin directory so my computer can see it wherever I am working, and in the process change the name to something more suitable - fft:

 cp demo ./fft && sudo mv ./fft /bin  

Right, it's a little bit more involved, but here is the bash shell script for the power spectrum in ImageMagick:

 #----------------------------- IMAGEMAGICK  

# read image and convert to greyscale
convert myimage.tif -colorspace gray input_grey.tiff
# find min dimension for cropping
min=`convert input_grey.tiff -format "%[fx: min(w,h)]" info:`
# crop image
convert input_grey.tiff -background white -gravity center -extent ${min}x${min} input_grey_crop.tiff
# de-mean image
# get data
data=`convert input_grey_crop.tiff -verbose info:`
# find mean as proportion
mean=`echo "$data" | sed -n '/^.*[Mm]ean:.*[(]\([0-9.]*\).*$/{ s//\1/; p; q; }'`
# convert mean to percent
m=`echo "scale=1; $mean *100 / 1" | bc`
# de-mean image
convert input_grey_crop.tiff -evaluate subtract ${m}% ggc_dm.tiff
# pad to power of 2
p2=`convert ggc_dm.tiff -format "%[fx:2^(ceil(log(max(w,h))/log(2)))]" info:`
convert ggc_dm.tiff -background white -gravity center -extent ${p2}x${p2} ggc_dm_pad.tiff
# forward fast fourier transform
fft f ggc_dm_pad.tiff
# creates ggc_dm_F_Q16.tiff, on which get spectrum
convert ggc_dm_F_Q16.tiff -contrast-stretch 0% -gamma .5 -fill "rgb(1,1,1)" -opaque black ggc_dm_F_Q16_spec.tiff

Voila! The spectrum is the image 'ggc_dm_F_Q16_spec.tiff'. I haven't yet worked out how to keep all the data inside the program, to eliminate the need (and cumulative errors associated with) writing the outputs at every stage to file, and eading them back in again. This means I can't use the spectrum for any scientific application. But in theory it should be possible. Any suggestions welcome!