ND Fourier Transform in Python

[EnSlavingBlair]

Hi,

My aim is to get a series of images in 2D space that run over different timestamps and put them through a 3D Fourier Transform. So my 3D FT has 2 spatial axes and one temporal axis. However I have never done anything like this before, and I have a very basic knowledge of Python.

So far, I can do the FFT for a list (or 1D array) of point sources as follows:

$# Import libs import matplotlib.pyplot as plt import numpy as np # create point sources tmp = range(100) source = [0 for x in tmp] source.insert(50,1) source.insert(5,1) source.insert(60,0.5) # make t useable later t = np.array(source) # equations for later f = np.fft.fft(t) g = np.sqrt(np.abs(f)**2) # set up plot fig = plt.figure() # add sources to plot ax1 = fig.add_subplot(211) plt.plot(source) ax1.set_title('Source') ax1.xaxis.set_visible(False) # add FT of sources to plot ax2 = fig.add_subplot(212) plt.plot(f) ax2.set_title('Fourier Transform') ax2.xaxis.set_visible(False) # show plot plt.show()$

Now I would like to turn my 1D image into a 2D image, and I just can't work out how to do this. I've tried by doing something like this for my point sources:

$# Array of 3 source lists # Create list of lists - seems dodgy tmp_array = range(3) array_list = [] # create source lists to go in array tmp = range(10) source1 = [0 for x in tmp] source1.insert(5,1) source2 = [0 for x in tmp] source2.insert(4,1) source3 = [0 for x in tmp] source3.insert(2,1) # put lists in array array_list.insert(1,source1) array_list.insert(2,source2) array_list.insert(3,source3)$

Though calling it "source" instead of "array_list" would fit better with previous code.
However it is not working and I cannot figure out why.

I was also wandering if I need to bother with getting the 2D FT working before trying the 3D, or if I can just jump forward? Not that I have any idea how to do that yet.

Thank you for your help

 

[Daverz]

I would suggest creating arrays as numpy arrays first, not as Python lists:

For the 1d array:

t = np.zeros(103)
t[5] = 1
t[51] = 1
t[60] = 0.5

(if that was your intention)

For the 2d array:

array_list = np.zeros((3, 11))
array_list[0,5] = 1
# etc...

For a 2d fft of with real-valued input, use rfft2 or rfftn.

Note that for large FFT sizes, try to avoid a size with large prime factors, or pad out to the next largest power of 2 if that is unavoidable (using the optional size parameter). Fftpack seems to handle small prime factors like 3 or 5 OK, though.

 

[Daverz]

One more comment:

f = np.fft.fft(t)
g = np.sqrt(np.abs(f)**2)

The sqrt and square here are unnecessary: np.abs(f) is already sqrt(f.real**2 + f.imag**2).

 

[EnSlavingBlair]

Thank you :)

I've gone for a slightly different approach now than above, and hitting different problems.

# A lot of this is taken from http://matplotlib.org/users/image_tutorial.html
# as well as http://opencv-python-tutroals.readthedocs.org/en/latest/py_tutorials/py_imgproc/py_transforms/py_fourier_transform/py_fourier_transform.html

$img=mpimg.imread('stinkbug.png') img_red = img[:,:,0] img_green = img[:,:,1] img_blue = img[:,:,2] # I assume this is taking the data in each of the red, green and blue columns? array_list = [] array_list.append(img_red) array_list.append(img_green) array_list.append(img_blue) t = np.array(array_list) f = np.fft.fftn(t) fshift = np.fft.fftshift(f) magnitude_spectrum = np.log(np.abs(fshift)) fig = plt.figure() ax1 = fig.add_subplot(211) imgplot=plt.imshow(img) ax1.set_title('FT this Bug') imgplot.set_cmap('gray') plt.colorbar() ax2 = fig.add_subplot(212) plt.imshow(magnitude_spectrum) plt.colorbar() plt.show()$

I get the following error

$# Traceback (most recent call last): File "<stdin>", line 1, in <module> File "image_test2", line 51, in <module> plt.imshow(magnitude_spectrum) File "/usr/lib/pymodules/python2.7/matplotlib/pyplot.py", line 2892, in imshow imlim=imlim, resample=resample, url=url, **kwargs) File "/usr/lib/pymodules/python2.7/matplotlib/axes.py", line 7300, in imshow im.set_data(X) File "/usr/lib/pymodules/python2.7/matplotlib/image.py", line 429, in set_data raise TypeError("Invalid dimensions for image data") TypeError: Invalid dimensions for image data$

Which I believe is because I'm trying to plot 3 images in 1 2D image, where I should be making a very thin rectangular cube. Because the images are basically the same, I expect that all the information will be in the first part and only left-over stuff in the others as I am sampling a limited sized image. I just don't know how to do that!

 

[Daverz]

Try plotting the red, green, and blue sections separately.

A few other comments on the code

Thank you :)
img=mpimg.imread('stinkbug.png')
img_red = img[:,:,0]
img_green = img[:,:,1]
img_blue = img[:,:,2]
# I assume this is taking the data in each of the red, green and blue columns?

array_list = []
array_list.append(img_red)
array_list.append(img_green)
array_list.append(img_blue)

t = np.array(array_list)

This could more easily be achieved using rollaxis(img, 2, 0):

img_rgba = img.rollaxis(img, 2, 0)
img_rgb = img_rgba[:3]

f = np.fft.fftn(t)
fshift = np.fft.fftshift(f)
magnitude_spectrum = np.log(np.abs(fshift))

Again, for real-valued input, the fft will have a second set of redundant conjugate negative-frequency values (the shift just moves them to to the beginning of the array, IIRC). I suggest using rfftn.