Properties of the Fourier Transform

import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
%matplotlib inline
plt.ion()
from scipy import ndimage, misc
from scipy.signal import gaussian

def Gaussian2D(dims, std):
    '''
        Creates an array containing a 2-D Gaussian function, with mean
        at mu, and standard deviation std.

        Inputs:
          dims is a 2-vector containing the number of rows and columns
          std is the standard deviation
    '''
    nr, nc = dims
    rr, cc = np.mgrid[0:nr,0:nc]
    ctr = GetCentre(rr)
    temp1 = rr - ctr[0]
    temp2 = cc - ctr[1]
    temp1 = np.exp(-temp1**2/std)
    temp2 = np.exp(-temp2**2/std)
    blah = temp1*temp2
    return blah/sum(blah.flatten())

def GetCentre(f):
    return np.array(np.floor(np.array(np.shape(f))/2.), dtype=int)

f = plt.imread('pd.jpg')[:,:,0]
plt.imshow(f,cmap='gray');

plt.clf()
F = np.zeros(128)
F[3] = 1.
#F[10] = 0.5
f = ifft(F)
plt.plot(np.real(f));

2D DFT

F = np.zeros([256,256])
F1 = F.copy()
#F1[0,2] = 1.
F1[5,10] = 1.0
#F[7,19] = 2.0
f1 = ifft2(F1)
plt.figure(figsize=[18,9])
plt.subplot(1,2,1)
plt.imshow(F1, cmap='gray')
plt.subplot(1,2,2)
plt.imshow(np.real(f1), cmap='gray');

F2 = F.copy()
F2[1,-2] = 1.
f2 = ifft2(F2)
plt.figure(figsize=[18,9])
plt.subplot(1,2,1)
plt.imshow(np.real(f2), cmap='gray')
plt.subplot(1,2,2)
plt.imshow(np.real(f1+f2), cmap='gray');

Convolution

rr, cc = np.mgrid[-128:128, -128:128]

Try a circular kernel

d = np.sqrt(rr**2 + cc**2)
circle = np.zeros([256,256])
circle_mask = d<15
circle[circle_mask] = 1.
circle = ifftshift(circle)
circle[0:22,0:5] = 1.
plt.figure(figsize=(8,4))
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('Shifted');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.title('Unshifted');

# Create an image with just a few non-zero pixels.
f = np.zeros([256,256])
f[50,100] = 1.0
f[48,110] = 2.
plt.figure(figsize=(10,10))
plt.imshow(f, cmap='gray');

g = ifft2( fft2(f) * fft2(circle))

plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('Circle convolved with impulses');

Or, try a Gaussian edge-detector kernel

# Make and edge detector
g = Gaussian2D([256,256], 20)
G = ifftshift( fft2(fftshift(g)) )
ramp = np.outer(range(256), range(256)) - 128**2
H = G * ramp*1.j
circle = np.real( ifft2( fftshift(H) ) )
plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('Shifted');

plt.imshow(circle, cmap='gray'); plt.title('Unshifted');

f = plt.imread('pd.jpg')
f = f[:,:,0]

plt.imshow(f, cmap='gray');

F = fftshift( fft2( f ) )

plt.figure(figsize=[10,20])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray')
plt.subplot(1,2,2)
plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.title('Shifted DFT (log-modulus)');

g = fftshift( ifft2( fft2(ifftshift(f)) * fft2(circle)) )

plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('Edge detector convolved with image');

G = fftshift(fft2(ifftshift(g)))
plt.imshow(np.abs(G), cmap='gray');

Filtering in the frequency domain

thresh = 10
rows, cols = np.shape(f)
ctr = np.floor(np.array(np.shape(f))/2) +1
rr, cc = np.mgrid[-ctr[0]:(rows-ctr[0]),
               -ctr[1]:(cols-ctr[1])]

# Choose a filter
f_type = 2
if f_type == 0:
    filt_mask = abs(cc)<thresh
elif f_type == 1:
    filt_mask = abs(rr)<thresh
elif f_type == 2:
    filt_mask = np.sqrt(rr**2 + cc**2)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0

F_filt = filt*F
f_filt = ifft2( fftshift(F_filt) )
plt.figure(figsize=[15,15])
plt.subplot(2,2,1); plt.imshow(f, cmap='gray');
plt.subplot(2,2,3); plt.imshow(np.log(abs(F)+1), cmap='gray')
plt.subplot(2,2,4); plt.imshow(np.log(abs(F_filt)+1), cmap='gray')
plt.subplot(2,2,2); plt.imshow(np.real(f_filt), cmap='gray');

Applying that filter in the frequency domain is the same as convolving the image with the kernel below, attained by taking the IDFT of the filter kernel.

Filt = fftshift(ifft2(ifftshift(filt)))
plt.figure()
plt.imshow(np.log(abs(Filt)+1), cmap='gray');

What happens if we take the DFT twice in a row?

plt.imshow(f, cmap='gray');

F = fft2(f)
plt.imshow(np.log(abs(F)), cmap='gray');

G = fft2(F)
plt.imshow(np.real(G), cmap='gray');

Scaling

You can scale an image using the DFT by padding in cropping in the two domains (time/space, and frequency).

# We will magnify an image by padding its DFT coefficients, and then cropping the image
# that results for the IDFT of those padded coefs.
scale_factor = 2
padamt = int(rows*(scale_factor-1.)/2.)
G = np.pad(ifftshift(F), (padamt,), mode='constant')
g512 = ifft2( ifftshift(G) )
if padamt!=0:
    g = g512[padamt:-padamt,padamt:-padamt]
else:
    g = g512

plt.figure(1,figsize=[15,15])
plt.clf()
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); plt.title('Scaled x'+str(scale_factor));

Fourier Shift Theorem

You can shift an image by adding a linear 'ramp' to the phase portion of its DFT coefficients.

dr = 20
dc = 12
ramp = -2.*np.pi*1j*(cc*dc + rr*dr)/rows
ramp = np.roll(np.roll(ramp,-1, axis=0),-1, axis=1)
wave = np.exp(ramp)
plt.figure()
plt.imshow(np.imag(ramp))
plt.title('Ramp in imaginary part');

# Make sure the centre pixel of the ramp has a phase of zero.
ctr = int(np.floor(256/2))
print('This should be zero -> '+str(ramp[ctr,ctr]))
print('This is the DC of our image -> '+str(F[ctr,ctr]))  # This should be the DC of our image.

This should be zero -> -0j
This is the DC of our image -> (-19+0j)

G = wave * fftshift(F)
g = ifft2(ifftshift(G))

plt.figure(1,figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray');
plt.title('Shifted (down,right)=('+str(dr)+','+str(dc)+')');

Phase Correlation

You can detect a shift between two similar images by looking for a ramp in their phase difference. Or, you can take that phase difference and take the IDFT of it, and look for a spike. The offset of the spike from the origin gives the relative shift.

g = ndimage.shift(f,[dr,dc]) + np.random.normal(0.0, 5, size=np.shape(f))
G = fft2(g)
F = fft2(f)
H = G/F
plt.imshow(np.real(fftshift(H)), cmap='gray', vmin=-4, vmax=4);

h = ifft2(H)
plt.imshow(np.real(h), cmap='gray');

# Where does that lonely pixel occur?
coords = np.unravel_index(np.argmax(np.real(h)), h.shape)
print('Bright pixel at '+str(coords))
print('True shift was '+str((dr,dc)))

Bright pixel at (20, 12)
True shift was (20, 12)

Rotation

You can also rotate an image by simply applying the same rotation to the frequency domain.

theta = 30
f = plt.imread('pd.jpg')
f = np.array(f[:,:,0], dtype=float)

padamt = int(rows/2)
#padamt = 0

ff = f
ff = np.pad(f, (padamt,), mode='constant')
FF = fftshift( fft2( ifftshift(ff) ) )

# Apply rotation to the Fourier coefficients
Gr = ndimage.rotate(np.real(FF), theta, reshape=False, order=1)
Gi = ndimage.rotate(np.imag(FF), theta, reshape=False, order=1)
G = Gr + 1j*Gi

g = fftshift( ifft2( ifftshift(G) ) )
if padamt!=0:
    g = g[padamt:-padamt,padamt:-padamt]

plt.figure(figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); 
plt.title('Rotated by '+str(theta)+' degrees');

Fourier Projection Theorem

The Fourier projection theorem is the basis for reconstruction in Computed Tomography (CT) scans.

f = plt.imread('ct.jpg')
f1 = np.array(f[:,:,0], dtype=float)

plt.figure(1, figsize=(10,10))
plt.clf()
plt.subplot(2,2,1)
plt.imshow(f1, cmap=plt.cm.gray);

px = np.sum(f1,axis=0)
plt.subplot(2,2,3)
plt.plot(px)
plt.axis('tight');

F = fftshift(fft2(f1))
plt.subplot(2,2,2)
Fnice = np.log(abs(F)+1)
plt.imshow(Fnice, cmap=plt.cm.gray);

ctr = GetCentre(F)
Fx = F[ctr[0],:]
plt.plot([1,254],[128,128], 'r--')
plt.axis([0,255,255,0])

plt.subplot(2,2,4)
plt.plot(np.log(abs(Fx)+1), 'r')
plt.axis('tight');

Px = fftshift( np.fft.fft(px) )
#plt.subplot(2,2,4)
plt.plot(np.log(abs(Px)+1), 'b--');

Aliasing

x, y = np.mgrid[0:256,0:256]
f = np.sin(2.*np.pi*x/256*40)
#f = ndimage.imread('pd.jpg')
plt.figure(figsize=(10,5))
plt.clf()
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray');
c = plt.axis()
plt.title('Original');

g = f[::5,::5]  # every fifth sample
plt.subplot(1,2,2)
plt.imshow(g, cmap='gray'); plt.axis(c);
plt.title('Subsampled');

q = plt.imread('bricks.jpg')
q = np.array(q[:,:,0], dtype=float)
nr, nc = np.shape(q)

plt.figure(figsize=(15,5));
plt.subplot(1,3,1)
plt.imshow(q, cmap='gray', interpolation='none')
plt.title('Original');

q_small = q[::5,::5]
qsr, qsc = np.shape(q_small)
plt.subplot(1,3,2)
plt.imshow(q_small, cmap='gray', interpolation='nearest')
plt.title('Subsampled')

g = Gaussian2D(q.shape, 10)
Q = fftshift( fft2(q) )
G = fftshift( fft2(fftshift(g)) )
h = ifft2(ifftshift(Q*G))

plt.subplot(1,3,3)
plt.imshow(np.real(h[::5,::5]), cmap='gray', interpolation='nearest')
plt.title('Filtered ');