.

# 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:

### 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.

### 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)

The temple image (frequency domain)

#### The Gaussian Kernel LPF in 2D (frequency domain)

https://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 smoothed temple image with the LPF (frequency domain)

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)

The output image after convolution (spatial / time domain)

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 &amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp; original image'``, size``=``20``)` `plt.axis(``'off'``)` `plt.show()`   `# Plot the surface of the frequency responses here`

Frequency response of the input image

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

https://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)

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

https://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)

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:

https://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 imagehttps://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)https://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

https://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)

https://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 imagehttps://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.

https://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.

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
https://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

Output image

https://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
https://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 imagehttps://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: 5762

Comment