Simple image denoising example using 2-dimensional FFT

In this example we illustrate a very elementary type of image denoising with the two-dimensional Fourier Transform.



In [3]:

    
%matplotlib inline

import sys
import numpy as np
from matplotlib import pyplot as plt



In [4]:

    
def plot_spectrum(F, amplify=1000, ax=None):
    """Normalise, amplify and plot an amplitude spectrum."""

    # Note: the problem here is that we have a spectrum whose histogram is
    # *very* sharply peaked at small values.  To get a meaningful display, a
    # simple strategy to improve the display quality consists of simply
    # amplifying the values in the array and then clipping.

    # Compute the magnitude of the input F (call it mag).  Then, rescale mag by
    # amplify/maximum_of_mag.
    mag = abs(F) 
    mag *= amplify/mag.max() 
    
    # Next, clip all values larger than one to one.
    mag[mag > 1] = 1 

    if ax is None: ax = plt.gca()
    ax.imshow(mag, plt.cm.Blues)

Read in original image, convert to floating point for further manipulation; imread returns a MxNx4 RGBA image. Since the image is grayscale, just extract the 1st channel

Hints:

use plt.imread() to load the file
convert to a float array with the .astype() method
extract all rows, all columns, 0-th plane to get the first channel
the resulting array should have 2 dimensions only



In [7]:

    
fname = 'images/moonlanding.png'
im = plt.imread(fname).astype(float) 
print("Image shape:", im.shape)









    



Image shape: (474, 630)

Compute the 2d FFT of the input image Hint: Look for a 2-d FFT in np.fft. Note: call this variable 'F', which is the name we'll be using below. In the lines following, we'll make a copy of the original spectrum and truncate coefficients.
In the lines following, we'll make a copy of the original spectrum and truncate coefficients.



In [8]:

    
F = np.fft.fft2(im)  

# Define the fraction of coefficients (in each direction) we keep
keep_fraction = 0.1

# Call ff a copy of the original transform.  Numpy arrays have a copy
# method for this purpose.
ff = F.copy() 

# Set r and c to be the number of rows and columns of the array.
r,c = ff.shape 

# Set to zero all rows with indices between r*keep_fraction and
# r*(1-keep_fraction):
ff[r*keep_fraction:r*(1-keep_fraction)] = 0  

# Similarly with the columns:
ff[:, c*keep_fraction:c*(1-keep_fraction)] = 0

Reconstruct the denoised image from the filtered spectrum, keep only the real part for display. Hint: There's an inverse 2d fft in the np.fft module as well (don't forget that you only want the real part). Call the result im_new and plot the results



In [9]:

    
im_new = np.fft.ifft2(ff).real  

fig, ax = plt.subplots(2, 2, figsize=(10,7))

ax[0,0].set_title('Original image')
ax[0,0].imshow(im, plt.cm.gray)

ax[0,1].set_title('Fourier transform')
plot_spectrum(F, ax=ax[0,1])

ax[1,1].set_title('Filtered Spectrum')
plot_spectrum(ff, ax=ax[1,1])

ax[1,0].set_title('Reconstructed Image')
ax[1,0].imshow(im_new, plt.cm.gray);