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  

 image=double(rgb2gray(imread('myimage.jpg')));  

 % crop to smallest dimension  

 [m,n]=size(im);  

 im=im(1:min(m,n),1:min(m,n));  

 % find image mean  

 imMean=mean(im(:));  

 % and subtract mean from the image  

 demean=im-imMean;  

 % find the 2D fourier transform  

 ft=fft2(demean);  

 % power spectrum is the element-by-element square of the fourier transform  

 spectrum=abs(ft).^2;  

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  

 im=Image.open('myimage.jpg').convert("L")  

 # crop to smallest dimension  

 box=(0,0,im.ny,im.nx)  

 region=im.crop(box)  

 # find mean  

 imMean=mean(asarray(region))  

 # 2D fourier transform of demeaned image  

 ft=fft2(region-imMean)  

 # power spectrum  

 spectrum=pow(abs(ft),2)  

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:

 ./configure  

 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  

 #-------------------------------------------------  

 #!/bin/bash  

 # 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!