Properties of the Fourier Transform¶

These demonstrations go beyond the content of the course. I don't necessarily expect you to understand all that is going on in these snippets of code. But I show them to you so you see what sorts of things can be implemented by manipulating the Fourier coefficients.

Preliminaries¶

In [1]:
import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
from scipy import ndimage, misc
from scipy.signal import gaussian
In [2]:
from IPython.display import display, HTML
display(HTML("<style>.container { width:98% !important; }</style>"))

np.set_printoptions(edgeitems=30, linewidth=250, 
    formatter=dict(float=lambda x: "% .4f" % x))
In [3]:
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)

Conjugate Symmetry¶

In [4]:
# If f is real-valued...
f = np.random.random(8)
print(f)

[ 0.4483  0.8407  0.0003  0.7059  0.0499  0.9091  0.3707  0.2379]
In [5]:
# Fourier coefs will be cojugate symmetric
F = np.fft.fft(f)
print(F)
print(F[2], F[-2])

[ 3.56264689+0.j          0.01917726+0.08782488j  0.12722956-0.8059954j   0.77776864-0.65306241j -1.8244074 +0.j          0.77776864+0.65306241j  0.12722956+0.8059954j   0.01917726-0.08782488j]
(0.12722956459780244-0.8059953973209955j) (0.12722956459780244+0.8059953973209955j)

Convolution¶

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

Try a circular kernel¶

In [7]:
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('on'); plt.title('g, Shifted');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('on'); plt.title('g, Unshifted');
In [8]:
# 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=(6,6))
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('f');
In [9]:
g = ifft2( fft2(f) * fft2(circle))

plt.figure(figsize=(6,6))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title(r'$(f \ast g)$');
In [10]:
f = plt.imread('images/pd.jpg')
f = f[:,:,0]
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('An image f');
In [11]:
g = ifft2( fft2(f) * fft2(circle))

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

Or, try a Gaussian edge-detector kernel¶

In [12]:
# 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.figure(figsize=[8,5]);
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Shifted)');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Unshifted)');
In [13]:
plt.imshow(np.imag(H), cmap='gray');
plt.title('Imaginary part of the DFT of the kernel (shifted)');
In [14]:
f = plt.imread('images/pd.jpg')
f = f[:,:,0]
In [15]:
F = fftshift( fft2( f ) )
In [16]:
plt.figure(figsize=[10,20])
plt.subplot(1,2,1)
plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off'); plt.title('Shifted DFT (log-modulus)');
plt.subplot(1,2,2)
plt.imshow(f, cmap='gray'); plt.axis('off');
In [17]:
g = fftshift( ifft2( fft2(ifftshift(f)) * fft2(circle)) )

plt.figure(figsize=(16,8))
plt.subplot(1,2,1);
G = fftshift(fft2(ifftshift(g)))
plt.imshow(np.abs(G), cmap='gray'); plt.axis('off')
plt.title('DFT(Edge Kernel) x DFT(Image)');
plt.subplot(1,2,2); plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('Edge Kernel convolved with image');

Filtering in the frequency domain¶

In [18]:
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])]
In [19]:
def PlotFiltImg(F, filt):
    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.axis('off');
    plt.subplot(2,2,3); plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off');
    plt.subplot(2,2,4); plt.imshow(np.log(abs(F_filt)+1), cmap='gray'); plt.axis('off');
    plt.subplot(2,2,2); plt.imshow(np.real(f_filt), cmap='gray'); plt.axis('off');
In [20]:
filt_mask = abs(cc)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [21]:
filt_mask = abs(rr)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [22]:
filt_mask = np.sqrt(rr**2 + cc**2)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [23]:
# Make and edge detector
g = Gaussian2D([256,256], 5)
G = ifftshift( fft2(fftshift(g)) )
# ramp = np.outer(range(256), np.ones(256)) - 128**2
# Hv = G * ramp*1.j
# ramp = np.outer(np.ones(256), range(256)) - 128**2
# Hh = G * ramp*1.j

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.

In [24]:
Filt = fftshift(ifft2(ifftshift(filt)))
plt.figure()
plt.subplot(1,2,1); plt.imshow(filt, cmap='gray'); plt.title('DFT of filter'); plt.axis('off')
plt.subplot(1,2,2); plt.imshow(np.log(abs(Filt)+1), cmap='gray'); plt.axis('off'); plt.title('Spatial filter');

Scaling¶

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

In [25]:
f = plt.imread('images/pd.jpg')
f = f[:,:,0]

# We will magnify an image by padding its DFT coefficients, and then cropping the image
# that results from the IDFT of those padded coefs.
scale_factor = 1.5
padamt = int(rows*(scale_factor-1.)/2.)
G = np.pad(fftshift(fft2(f)), (padamt,), mode='constant')
g512 = ifft2( ifftshift(G) )
if padamt!=0:
    g = g512[padamt:-padamt,padamt:-padamt]
else:
    g = g512
In [26]:
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.

In [27]:
dr = 44
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(figsize=(10,5))
plt.subplot(1,2,1); plt.imshow(np.real(wave))
plt.subplot(1,2,2); plt.imshow(np.imag(ramp))
plt.title('Ramp in imaginary part');
In [28]:
# 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 -> (2491737+0j)
In [29]:
G = wave * fftshift(fft2(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)+')');

Pendulum synchronization video (YouTube)

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.

In [30]:
plt.figure(figsize=(8,8));
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);
In [31]:
plt.figure(figsize=(8,8));
h = ifft2(H)
plt.imshow(np.real(h), cmap='gray');
In [32]:
# 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 (44, 12)
True shift was (44, 12)

Rotation¶

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

In [33]:
theta = 30
f = plt.imread('images/pd.jpg')
f = np.array(f[:,:,0], dtype=float)
In [34]:
padamt = int(rows/2)
#padamt = 0
In [41]:
# If we pad the image first, it makes the frequency domain smoother.
ff = f
ff = np.pad(f, (padamt,), mode='constant')
FF = fftshift( fft2( ifftshift(ff) ) )
In [36]:
# 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
In [37]:
g = fftshift( ifft2( ifftshift(G) ) )
if padamt!=0:
    g = g[padamt:-padamt,padamt:-padamt]
In [38]:
plt.figure(figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original'); plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); 
plt.title('Rotated by '+str(theta)+' degrees'); plt.axis('off');

Fourier Projection Theorem¶

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

In [39]:
f = plt.imread('images/ct.jpg')
f1 = np.array(f[:,:,0], dtype=float)
In [42]:
plt.figure(1, figsize=(10,10))
plt.clf()
plt.subplot(2,2,1)
plt.imshow(f1, cmap=plt.cm.gray);
plt.title('CT image')

px = np.sum(f1,axis=0)
plt.subplot(2,2,3)
plt.plot(px); plt.title('Projection')
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);
plt.title('2D-DFT, and slice')

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'); plt.ylabel('log modulus')

Px = fftshift( np.fft.fft(px) )
#plt.subplot(2,2,4)
plt.plot(np.log(abs(Px)+1), 'b--');
plt.legend(['Slice through 2D-DFT', '1D-DFT of Projection'])
Out[42]:
<matplotlib.legend.Legend at 0x1907d49d0>
In [ ]: