Subscribe to DSC Newsletter

Solving Some Image Processing Problems with Python libraries - Part 2

In this article a few more popular image processing problems along with their solutions are going to be discussed. Python image processing libraries are going to be used to solve these problems.

Some Applications of DFT 

0. Fourier Transform of a Gaussian Kernel is another Gaussian Kernel

Also, the spread in the frequency domain  inversely proportional to the spread in the spatial domain. Here is the proof:

The following animation shows an example visualizing the Gaussian contours in spatial and corresponding frequency domains:

gaussian.gif

 

1. Using DFT to up-sample an image

  1. Let’s use the lena gray-scale image.
  2. First double the size of the by padding zero rows/columns at every alternate positions.
  3. Use FFT followed by an LPF.
  4. Finally use IFFT to get the output image.

 

The following code block shows the python code for implementing the steps listed above:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
import numpy as np
import numpy.fft as fp
import matplotlib.pyplot as plt
 
im = np.mean(imread('images/lena.jpg'), axis=2)
im1 = np.zeros((2*im.shape[0], 2*im.shape[1]))
print(im.shape, im1.shape)
for i in range(im.shape[0]):
    for j in range(im.shape[1]):
        im1[2*i,2*j] = im[i,j]
 
def padwithzeros(vector, pad_width, iaxis, kwargs):
    vector[:pad_width[0]] = 0
    vector[-pad_width[1]:] = 0
    return vector
 
# the LPF kernel
kernel = [[0.25, 0.5, 0.25], [0.5, 1, 0.5], [0.25, 0.5, 0.25]]
# enlarge the kernel to the shape of the image
kernel = np.pad(kernel, (((im1.shape[0]-3)//2,(im1.shape[0]-3)//2+1), ((im1.shape[1]-3)//2,(im1.shape[1]-3)//2+1)), padwithzeros)
 
plt.figure(figsize=(15,10))
plt.gray() # show the filtered result in grayscale
 
freq = fp.fft2(im1)
freq_kernel = fp.fft2(fp.ifftshift(kernel))
freq_LPF = freq*freq_kernel # by the Convolution theorem
im2 = fp.ifft2(freq_LPF)
freq_im2 = fp.fft2(im2)
 
plt.subplot(2,3,1)
plt.imshow(im)
plt.title('Original Image', size=20)
plt.subplot(2,3,2)
plt.imshow(im1)
plt.title('Padded Image', size=20)
plt.subplot(2,3,3)
plt.imshow( (20*np.log10( 0.1 + fp.fftshift(freq))).astype(int), cmap='jet')
plt.title('Original Image Spectrum', size=20)
plt.subplot(2,3,4)
plt.imshow( (20*np.log10( 0.1 + fp.fftshift(freq_kernel))).astype(int), cmap='jet')
plt.title('Image Spectrum of the LPF', size=20)
plt.subplot(2,3,5)
plt.imshow( (20*np.log10( 0.1 + fp.fftshift(freq_im2))).astype(int), cmap='jet')
plt.title('Image Spectrum after LPF', size=20)
plt.subplot(2,3,6)
plt.imshow(im2.astype(np.uint8)) # the imaginary part is an artifact
plt.title('Output Image', size=20)

The next figure shows the output. As can be seen from the next figure, the LPF removed the high frequency components from the Fourier spectrum of the padded image and with a subsequent inverse Fourier transform  we get a decent enlarged image.

lena_LPF_fft.pnghttps://sandipanweb.files.wordpress.com/2018/07/lena_lpf_fft.png?w=150 150w, https://sandipanweb.files.wordpress.com/2018/07/lena_lpf_fft.png?w=300 300w, https://sandipanweb.files.wordpress.com/2018/07/lena_lpf_fft.png?w=768 768w, https://sandipanweb.files.wordpress.com/2018/07/lena_lpf_fft.png 881w" sizes="(max-width: 676px) 100vw, 676px" />

 

2. Frequency Domain Gaussian Filter

  1. Use an input image and use DFT to create the frequency 2D-array.
  2. Create a small Gaussian 2D Kernel (to be used as an LPF) in the spatial domain and pad it to enlarge it to the image dimensions.
  3. Use DFT to obtain the Gaussian Kernel in the frequency domain.
  4. Use the Convolution theorem to convolve the LPF with the input image in the frequency domain.
  5. Use IDFT to obtain the output image.
  6. Plot the frequency spectrum of the image, the gaussian kernel and the image obtained after convolution in the frequency domain, in 3D.

The following code block shows the python code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import matplotlib.pyplot as plt
from matplotlib import cm
from skimage.color import rgb2gray
from skimage.io import imread
import scipy.fftpack as fp
 
im = rgb2gray(imread('images/temple.jpg'))
kernel = np.outer(signal.gaussian(im.shape[0], 10), signal.gaussian(im.shape[1], 10))
freq = fp.fft2(im)
assert(freq.shape == kernel.shape)
freq_kernel = fp.fft2(fp.ifftshift(kernel))
convolved = freq*freq_kernel # by the Convolution theorem
im_blur = fp.ifft2(convolved).real
im_blur = 255 * im_blur / np.max(im_blur)
 
# center the frequency response
plt.imshow( (20*np.log10( 0.01 + fp.fftshift(freq_kernel))).real.astype(int), cmap='coolwarm')
plt.colorbar()
plt.show()
 
plt.figure(figsize=(20,20))
plt.imshow(im, cmap='gray')
plt.show()
 
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# ... code for 3D visualization of the spectrums

The original color temple image (time / spatial domain)

temple.JPGhttps://sandipanweb.files.wordpress.com/2018/07/temple.jpg?w=1352 1352w, https://sandipanweb.files.wordpress.com/2018/07/temple.jpg?w=150 150w, https://sandipanweb.files.wordpress.com/2018/07/temple.jpg?w=300 300w, https://sandipanweb.files.wordpress.com/2018/07/temple.jpg?w=768 768w, https://sandipanweb.files.wordpress.com/2018/07/temple.jpg?w=1024 1024w" sizes="(max-width: 676px) 100vw, 676px" />

The temple image (frequency domain)

temple_freqhttps://sandipanweb.files.wordpress.com/2018/07/temple_freq.png?w=150 150w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq.png?w=300 300w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq.png?w=768 768w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq.png?w=1024 1024w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq.png 1137w" sizes="(max-width: 676px) 100vw, 676px" />

The Gaussian Kernel LPF in 2D (frequency domain)

gaussian_freq2dhttps://sandipanweb.files.wordpress.com/2018/07/gaussian_freq2d.png... 150w, https://sandipanweb.files.wordpress.com/2018/07/gaussian_freq2d.png... 300w" sizes="(max-width: 801px) 100vw, 801px" />

The Gaussian Kernel LPF (frequency domain)
gaussian_freqhttps://sandipanweb.files.wordpress.com/2018/07/gaussian_freq.png?w... 150w, https://sandipanweb.files.wordpress.com/2018/07/gaussian_freq.png?w... 300w, https://sandipanweb.files.wordpress.com/2018/07/gaussian_freq.png?w... 768w, https://sandipanweb.files.wordpress.com/2018/07/gaussian_freq.png?w... 1024w, https://sandipanweb.files.wordpress.com/2018/07/gaussian_freq.png 1137w" sizes="(max-width: 676px) 100vw, 676px" />

The smoothed temple image with the LPF (frequency domain)

temple_freq_LPFhttps://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf.png... 150w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf.png... 300w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf.png... 768w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf.png... 1024w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf.png 1137w" sizes="(max-width: 676px) 100vw, 676px" />

If we set the standard deviation of the LPF Gaussian kernel to be 10 we get the following output as shown in the next figures. As can be seen, the frequency response value drops much quicker from the center.

The smoothed temple image with the LPF with higher s.d. (frequency domain)

temple_freq_LPF2https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf2.pn... 150w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf2.pn... 300w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf2.pn... 768w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf2.pn... 1024w, https://sandipanweb.files.wordpress.com/2018/07/temple_freq_lpf2.png 1137w" sizes="(max-width: 676px) 100vw, 676px" />

The output image after convolution (spatial / time domain)

temple_blurredhttps://sandipanweb.files.wordpress.com/2018/07/temple_blurred.png?... 150w, https://sandipanweb.files.wordpress.com/2018/07/temple_blurred.png?... 300w, https://sandipanweb.files.wordpress.com/2018/07/temple_blurred.png?... 768w, https://sandipanweb.files.wordpress.com/2018/07/temple_blurred.png?... 1024w, https://sandipanweb.files.wordpress.com/2018/07/temple_blurred.png 1166w" sizes="(max-width: 676px) 100vw, 676px" />

3. Using the inverse filter to restore a motion-blurred image

  1. First create a motion blur kernel of a given shape.
  2. Convolve the kernel with an input image in the frequency domain.
  3. Get the motion-blurred image in the spatial domain with IDFT.
  4. Compute the inverse filter kernel and convolve with the blurred image in the frequency domain.
  5. Get the convolved image back in the spatial domain.
  6. Plot all the images and kernels in the frequency domain.

The following code block shows the python code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
im = rgb2gray(imread('../my images/madurga.jpg'))
# create the motion blur kernel
size = 21
kernel = np.zeros((size, size))
kernel[int((size-1)/2), :] = np.ones(size)
kernel = kernel / size
kernel = np.pad(kernel, (((im.shape[0]-size)//2,(im.shape[0]-size)//2+1), ((im.shape[1]-size)//2,(im.shape[1]-size)//2+1)), padwithzeros)
 
freq = fp.fft2(im)
freq_kernel = fp.fft2(fp.ifftshift(kernel))
convolved1 = freq1*freq_kernel1
im_blur = fp.ifft2(convolved1).real
im_blur = im_blur / np.max(im_blur)
 
epsilon = 10**-6
 
freq = fp.fft2(im_blur)
freq_kernel = 1 / (epsilon + freq_kernel1)
 
convolved = freq*freq_kernel
im_restored = fp.ifft2(convolved).real
im_restored = im_restored / np.max(im_restored)
 
plt.figure(figsize=(18,12))
plt.subplot(221)
plt.imshow(im)
plt.title('Original image', size=20)
plt.axis('off')
plt.subplot(222)
plt.imshow(im_blur)
plt.title('Blurred image with motion blur kernel', size=20)
plt.axis('off')
plt.subplot(223)
plt.imshow(im_restored)
plt.title('Restored image with inverse filter', size=20)
plt.axis('off')
plt.subplot(224)
plt.imshow(im_restored - im)
plt.title('Diff restored & original image', size=20)
plt.axis('off')
plt.show()
 
# Plot the surface of the frequency responses here

madurga_inversehttps://sandipanweb.files.wordpress.com/2018/07/madurga_inverse.png... 150w, https://sandipanweb.files.wordpress.com/2018/07/madurga_inverse.png... 300w, https://sandipanweb.files.wordpress.com/2018/07/madurga_inverse.png... 768w, https://sandipanweb.files.wordpress.com/2018/07/madurga_inverse.png 931w" sizes="(max-width: 676px) 100vw, 676px" />

Frequency response of the input image

madurga_freq.pnghttps://sandipanweb.files.wordpress.com/2018/07/madurga_freq.png?w=... 150w, https://sandipanweb.files.wordpress.com/2018/07/madurga_freq.png?w=... 300w" sizes="(max-width: 630px) 100vw, 630px" />

(log) Frequency response of the motion blur kernel (LPF)

motion_blur_freq.pnghttps://sandipanweb.files.wordpress.com/2018/07/motion_blur_freq1.p... 150w, https://sandipanweb.files.wordpress.com/2018/07/motion_blur_freq1.p... 300w" sizes="(max-width: 693px) 100vw, 693px" />

Input image convolved with the motion blur kernel (frequency domain)

madurga_convolved1https://sandipanweb.files.wordpress.com/2018/07/madurga_convolved11... 150w, https://sandipanweb.files.wordpress.com/2018/07/madurga_convolved11... 300w" sizes="(max-width: 637px) 100vw, 637px" />

(log) Frequency response of the inverse frequency filter kernel (HPF)

motion_blur_freqhttps://sandipanweb.files.wordpress.com/2018/07/motion_blur_freq.pn... 150w, https://sandipanweb.files.wordpress.com/2018/07/motion_blur_freq.pn... 300w" sizes="(max-width: 636px) 100vw, 636px" />

Motion-blurred image convolved with the inverse frequency filter kernel (frequency domain)

madurga_convolved.pnghttps://sandipanweb.files.wordpress.com/2018/07/madurga_convolved.p... 150w, https://sandipanweb.files.wordpress.com/2018/07/madurga_convolved.p... 300w" sizes="(max-width: 583px) 100vw, 583px" />

 

4. Impact of noise on the inverse filter

  1. Add some random noise to the Lena image.
  2. Blur the image with a Gaussian kernel.
  3. Restore the image using inverse filter.

With the original image

Let’s first blur and apply the inverse filter on the noiseless blurred image. The following figures show the outputs:

lena_inversehttps://sandipanweb.files.wordpress.com/2018/07/lena_inverse.png?w=150 150w, https://sandipanweb.files.wordpress.com/2018/07/lena_inverse.png?w=300 300w, https://sandipanweb.files.wordpress.com/2018/07/lena_inverse.png 745w" sizes="(max-width: 676px) 100vw, 676px" />

(log) Frequency response of the input imagelena_freq.pnghttps://sandipanweb.files.wordpress.com/2018/07/lena_freq.png?w=150 150w, https://sandipanweb.files.wordpress.com/2018/07/lena_freq.png?w=300 300w" sizes="(max-width: 579px) 100vw, 579px" />

(log) Frequency response of the Gaussian blur kernel (LPF)lena_gaussian_freqhttps://sandipanweb.files.wordpress.com/2018/07/lena_gaussian_freq.... 150w, https://sandipanweb.files.wordpress.com/2018/07/lena_gaussian_freq.... 300w" sizes="(max-width: 635px) 100vw, 635px" />

(log) Frequency response of the blurred image

lena_blurred_spectrum.pnghttps://sandipanweb.files.wordpress.com/2018/07/lena_blurred_spectr... 150w, https://sandipanweb.files.wordpress.com/2018/07/lena_blurred_spectr... 300w" sizes="(max-width: 621px) 100vw, 621px" />

(log) Frequency response of the inverse kernel (HPF)

gauss_inverse_kernel.pnghttps://sandipanweb.files.wordpress.com/2018/07/gauss_inverse_kerne... 150w, https://sandipanweb.files.wordpress.com/2018/07/gauss_inverse_kerne... 300w" sizes="(max-width: 639px) 100vw, 639px" />

Frequency response of the output imagelena_inverse_spectrumhttps://sandipanweb.files.wordpress.com/2018/07/lena_inverse_spectr... 150w, https://sandipanweb.files.wordpress.com/2018/07/lena_inverse_spectr... 300w" sizes="(max-width: 563px) 100vw, 563px" />

Adding noise to the original image

The following python code can be used to add Gaussian noise to an image:

1
2
from skimage.util import random_noise
im = random_noise(im, var=0.1)


The next figures show the noisy lena image, the blurred image with a Gaussian Kernel and the restored image with the inverse filter. As can be seen, being a high-pass filter, the inverse filter enhances the noise, typically corresponding to high frequencies.

 

lena_noisy_inversehttps://sandipanweb.files.wordpress.com/2018/07/lena_noisy_inverse1... 150w, https://sandipanweb.files.wordpress.com/2018/07/lena_noisy_inverse1... 300w, https://sandipanweb.files.wordpress.com/2018/07/lena_noisy_inverse1... 745w" sizes="(max-width: 676px) 100vw, 676px" />

5. Use a notch filter to remove periodic noise from the following half-toned car image.

halftone.pnghttps://sandipanweb.files.wordpress.com/2018/07/halftone.png?w=150 150w, https://sandipanweb.files.wordpress.com/2018/07/halftone.png?w=300 300w, https://sandipanweb.files.wordpress.com/2018/07/halftone.png?w=768 768w, https://sandipanweb.files.wordpress.com/2018/07/halftone.png 1000w" sizes="(max-width: 676px) 100vw, 676px" />

  1. Use DFT to obtain the frequency spectrum of the image.
  2. Block the high frequency components that are most likely responsible fro noise.
  3. Use IDFT to come back to the spatial domain.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
from scipy import fftpack
im = imread('images/halftone.png')
F1 = fftpack.fft2((im).astype(float))
F2 = fftpack.fftshift(F1)
for i in range(60, w, 135):
for j in range(100, h, 200):
if not (i == 330 and j == 500):
F2[i-10:i+10, j-10:j+10] = 0
for i in range(0, w, 135):
for j in range(200, h, 200):
if not (i == 330 and j == 500):
F2[max(0,i-15):min(w,i+15), max(0,j-15):min(h,j+15)] = 0
plt.figure(figsize=(6.66,10))
plt.imshow( (20*np.log10( 0.1 + F2)).astype(int), cmap=plt.cm.gray)
plt.show()
im1 = fp.ifft2(fftpack.ifftshift(F2)).real
plt.figure(figsize=(10,10))
plt.imshow(im1, cmap='gray')
plt.axis('off')
plt.show()


Frequency response of the input image
car_freqhttps://sandipanweb.files.wordpress.com/2018/07/car_freq.png?w=150 150w, https://sandipanweb.files.wordpress.com/2018/07/car_freq.png?w=300 300w, https://sandipanweb.files.wordpress.com/2018/07/car_freq.png?w=768 768w, https://sandipanweb.files.wordpress.com/2018/07/car_freq.png?w=1024 1024w, https://sandipanweb.files.wordpress.com/2018/07/car_freq.png 1160w" sizes="(max-width: 676px) 100vw, 676px" />

Frequency response of the input image with blocked frequencies with notch

car_notch_blockedhttps://sandipanweb.files.wordpress.com/2018/07/car_notch_blocked.p... 150w, https://sandipanweb.files.wordpress.com/2018/07/car_notch_blocked.p... 300w, https://sandipanweb.files.wordpress.com/2018/07/car_notch_blocked.p... 768w, https://sandipanweb.files.wordpress.com/2018/07/car_notch_blocked.p... 1024w, https://sandipanweb.files.wordpress.com/2018/07/car_notch_blocked.png 1160w" sizes="(max-width: 676px) 100vw, 676px" />

Output image

car_notchhttps://sandipanweb.files.wordpress.com/2018/07/car_notch1.png?w=15... 150w, https://sandipanweb.files.wordpress.com/2018/07/car_notch1.png?w=30... 300w" sizes="(max-width: 718px) 100vw, 718px" />

With a low-pass-filter (LPF):

Frequency response of the input image with blocked frequencies with LPF
car_freq_notchhttps://sandipanweb.files.wordpress.com/2018/07/car_freq_notch.png?... 150w, https://sandipanweb.files.wordpress.com/2018/07/car_freq_notch.png?... 300w" sizes="(max-width: 696px) 100vw, 696px" />

Output imagecar_notchhttps://sandipanweb.files.wordpress.com/2018/07/car_notch.png?w=150... 150w, https://sandipanweb.files.wordpress.com/2018/07/car_notch.png?w=300... 300w" sizes="(max-width: 775px) 100vw, 775px" />

Views: 1619

Comment

You need to be a member of Data Science Central to add comments!

Join Data Science Central

Videos

  • Add Videos
  • View All

© 2019   Data Science Central ®   Powered by

Badges  |  Report an Issue  |  Privacy Policy  |  Terms of Service