2-dimensional Fourier transform interpretation

fourier analysisfourier transform

I have this image of a duck:
enter image description here

Taking the FFT of this image, I can generate two new images. A magnitude spectrum (right) and phase spectrum (bottom left).

I am trying to re-construct the image of the duck MANUALLY by adding together the 65,535 sine waves encoded in this information. I am running into a problem where my duck is reconstructed but looks a little funny.

If the 2D Fourier transform is as follows:

$f(x,y) = \sum\sum F[h,k]e^{-2\pi i(hx+ky+\phi)}$

I was hoping to be schooled in my interpretation of this. My sine wave amplitude is $F[h,k]$. My phase is $\phi$. My indices for the wave (the number of time it 'cuts' each respective axis are $h, k$.

So, what the magnitude spectrum is showing is:

  1. The center is $h,k = 0,0$. The zero frequency component intensity represents the average intensity of the entire image.
  2. Each coordinate, for example $(2, 2)$ from the center represents a wave with $(h, k) = (2, 2)$. The intensity of this pixel is an indication of the amplitude of the wave.
  3. Similarly, in the phase spectrum, the values range between $(0, \pi)$ and represent the offset of this sine wave from the origin.

The code I am using to doing the following is below. I SHOULD be able to perfectly reconstruct the image, but I keep getting some weird looking duck. The following is the duck I get after adding together 65,535 waves:

enter image description here

I was hoping somebody could point out my error. Thank you!

EDIT

If I only add the first 128 rows of my 2D FFT arrays, I get a duck that looks closer to the original. I feel like it has something to do with the fact that you only need 1/2 the information from these spectra, but I am adding together some of the wrong info. Like I need to cut each spectrum diagonally (about the symmetrical axis) and add only those waves…

EDIT

Final product…

enter image description here

import numpy as np
import matplotlib.pyplot as plt
import cv2

#Wave generator
def waveGenerator(h, k, a, phi, res):
    mesh = np.fromfunction(lambda x, y: a*np.sin(2*np.pi*h*x/res + 2*np.pi*k*y/res + phi), shape = (res,res))
    return mesh

# Import some image...
duck = cv2.imread("data/fourier_duck.png", 0)

#Take the fft of the duck
duck_fft = np.fft.fft2(duck)
duck_fshift = np.fft.fftshift(duck_fft)
magnitude_duck = np.log(np.abs(duck_fshift) + 1)
unshifted = np.abs(duck_fft)
phases = np.angle(duck_fft)

#THE RECONSTRUCTED IMAGE (MANUALLY)
recon_img = np.full((256, 256), unshifted[0][0])

for h in range(len(unshifted)):
    for k in range(len(unshifted[h])):
        recon_img = np.add(recon_img, waveGenerator(h, k, unshifted[h][k], phases[h][k], 256))
        print(str(k) + ',' + str(h))

plt.imshow(recon_img)
plt.show()

Best Answer

Try changing:

a*np.sin

to

a*np.cos

in your waveGenerator function. When the phase is 0, the real part of the wavefunction should be maximum which corresponds to a cosine: i.e. cos(0)=1. The rest seems to be okay, other than scaling issues I believe. I hope this helps.

Simplified test program below adapted from the question:

import numpy as np
import matplotlib.pyplot as plt
#import cv2

#Wave generator
def waveGenerator(h, k, a, phi, res):
    mesh = np.fromfunction(lambda x, y: a*np.cos(2*np.pi*h*x/res + 2*np.pi*k*y/res + phi), shape = (res,res))
    return mesh

# Import some image...
#duck = cv2.imread("data/fourier_duck.png", 0)
n_size = 64
duck = np.zeros((n_size,n_size))
duck[n_size//2][n_size//2]=1.0
duck[(n_size+10)//2][(n_size+5)//4]=2.0
plt.imshow(duck)
plt.show()

#Take the fft of the duck
duck_fft = np.fft.fft2(duck)
duck_fshift = np.fft.fftshift(duck_fft)
magnitude_duck = np.log(np.abs(duck_fshift) + 1)
unshifted = np.abs(duck_fft)
phases = np.angle(duck_fft)

#THE RECONSTRUCTED IMAGE (MANUALLY)
recon_img = np.full((n_size, n_size), unshifted[0][0])

for h in range(len(unshifted)):
    for k in range(len(unshifted[h])):
        recon_img = np.add(recon_img, waveGenerator(h, k, unshifted[h][k], phases[h][k], n_size))

plt.imshow(recon_img)
plt.show()
Related Question