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!
