Filters and Denoising Introduction

This tutorial demonstrates convolution using the scikit-image library. The tutorial covers:

  • Convolution/correlation with manual and predefined filters.
  • Explicitly construct filters to perform specific operations.

As a first step, it is necessary to import the libraries:

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import convolve
from skimage import data, filters
from scipy import signal
import scipy
import scipy.ndimage as sp
import seaborn as sns
import pandas as pd
from skimage import data, color
from scipy.signal import convolve2d

Scikit-Image has a set of example images (see the documentation for a complete list). The cell below prints the name of all images and displays the rocket image:

print(data.__all__)
plt.imshow(data.rocket(), cmap="gray")
plt.axis(False)
['astronaut', 'binary_blobs', 'brain', 'brick', 'camera', 'cat', 'cell', 'cells3d', 'checkerboard', 'chelsea', 'clock', 'coffee', 'coins', 'colorwheel', 'create_image_fetcher', 'data_dir', 'download_all', 'eagle', 'file_hash', 'grass', 'gravel', 'horse', 'hubble_deep_field', 'human_mitosis', 'image_fetcher', 'immunohistochemistry', 'kidney', 'lbp_frontal_face_cascade_filename', 'lfw_subset', 'lily', 'logo', 'microaneurysms', 'moon', 'nickel_solidification', 'page', 'protein_transport', 'retina', 'rocket', 'shepp_logan_phantom', 'skin', 'stereo_motorcycle', 'text', 'vortex']

Linear filters are implemented as inner products and have many applications, such as blurring, edge detection or finding gradients and are essential in certain types of neural networks. Several of the filters presented at the lecture can be obtained and applied using the skimage.filters module. For example, a gaussian filter can be applied to an image using the filters.gaussian function as shown below. The documentation contains a complete list of filters. You can expermient with other filters and observe the results. 

blurred = filters.gaussian(data.rocket(), channel_axis=2, sigma=5)
plt.imshow(blurred, cmap="gray")
plt.axis(False)

Filters from scratch

In image processing a filter is an $N \times N$ matrix with values that are used as weights applied to the corresponding pixel and its neighboring pixels (dot product). The mean filter (matrix) has filter values of ones and calculates the average over a region specified by the filter size. The size of the filter determines the extent of smoothing/blurring applied to the signal. A larger filter will lead to more smoothing by averaging over a larger area. The function mean_filter in the cell below produces a mean filter:

def create_mean_filter(size):
    return np.ones((size, size, 1))/(size**2)

mean_filter =  create_mean_filter(10)

The cell below plots the mean filter:

# Create a 10x10 grid to display the filter values
plt.figure(figsize=(8, 8))
plt.imshow(mean_filter[:, :, 0], cmap='gray', interpolation='none', vmin=0, vmax=0.01, extent=[0, 10, 0, 10])  # Adjust vmin and vmax for your specific filter values

# Display the values inside each cell
for i in range(mean_filter.shape[0]):
    for j in range(mean_filter.shape[1]):
        plt.text(j + 0.5, i + 0.5, f'{mean_filter[i, j, 0]:.2f}', color='black', ha='center', va='center')

plt.title('Mean Filter (10x10)')
plt.xticks(np.arange(0, 11, 1))
plt.yticks(np.arange(0, 11, 1))
plt.grid(color='black', linestyle='-', linewidth=0.5)
plt.show()

When the filter is convolved with an image (applying the dot product iteratively on each color channel), each pixel in the output image represents the average of its surrounding pixels, producing a smoothed image. The cell below visualizes the image of the rocket ship. The pixel values have been normalized and downsampled to 20x20 pixels for display purposes: 

# Load the "camera" image from scikit-learn
data1 = data.rocket()
from skimage.transform import resize

# Resize the image to a smaller size
target_height, target_width = 20, 20  # Set your desired dimensions
resized_image = resize(data1, (target_height, target_width))

# Create a figure with two subplots
fig, axes = plt.subplots(1, 2, figsize=(15, 8))

# Plot the grid of values in the first subplot
axes[0].imshow(np.ones_like(resized_image), cmap='gray', extent=[0, 20, 0, 20], vmin=0, vmax=1)  # Display white cells

# Display the values strictly inside each cell
for i in range(target_height):
    for j in range(target_width):
        pixel_value = resized_image[i, j, 0]
        axes[0].text(j + 0.5, target_height - i - 0.5, f'{pixel_value:.2f}', color='black', ha='center', va='center')

axes[0].set_title('Grid of Values (20x20)')
axes[0].set_xticks(np.arange(0, 21, 1))
axes[0].set_yticks(np.arange(0, 21, 1))
axes[0].grid(color='black', linestyle='-', linewidth=0.5)

# Plot the downsampled image in the second subplot without pixel values
axes[1].imshow(resized_image, cmap='hot', vmin=0, vmax=1)
axes[1].set_title('Downsampled Image (20x20)')
axes[1].axis('off')  # Hide axis for the image plot

plt.tight_layout()
plt.show()

The cell below convolves the mean filter to the rocket image using scipy.ndimage.convolve :

blurred = scipy.ndimage.convolve(data.rocket(), mean_filter)
plt.imshow(blurred)
plt.axis(False)

Edge Filter by convolution

Edges are places in the images with large differences between intensity/color values and correspond to intensity gradients. By constructing a filter that calculates the partial derivatives of an image in either the $x$ or $y$ direction linear filters can be used to find edges in an image in a certain direction.

The definition of the partial derivative of a function $f$ with respect to $x$ is defined as:

$$ \begin{align} \frac{\partial f(x, y)}{\partial x} = \lim_{\Delta \rightarrow 0} \frac{f(x+\Delta, y) - f(x,y)}{\Delta} \end{align} $$

An image is a discrete approximation of the light distribution hitting the camera sensor with integer indices. The partial derivative can be approximated by setting $\Delta$ to $1$ (smallest possible step). This method is known as the finite differences method:

$$ \begin{align} \frac{\partial f(x, y)}{\partial x} \approx \frac{f(x+1, y) - f(x,y)}{1} = f(x+1, y) - f(x, y) \end{align} $$

This expression is the difference between neighboring pixels. The corresponding filter to calculate the derivative with respect to $x$ is:

$$ k_{x} = \begin{bmatrix}1&-1\end{bmatrix} $$

Equivalently, the filter for calculating the derivative in the $y$-direction is:

$$ \begin{align*} k_{y} &= \begin{bmatrix}1 \\ -1\end{bmatrix} \end{align*} $$

You may have noticed that the $1$ and $-1$ seem reversed in the filter. This is due the definition of convolution. The cell below implements the filter and convolves it with the image . The image is grayscale but displayed using brown colors for illustrative purposes:

k_x = np.array([[1, -1]])

image = data.camera().astype(np.float32)

dx = convolve(image, k_x)

# Create a 1x2 grid of subplots
fig, axs = plt.subplots(1, 2, figsize=(10, 5))

# Plot the first image in the first subplot
axs[0].imshow(data.camera(), cmap='copper')
axs[0].set_title('Original Image')
axs[0].axis('off')

# Plot the result of convolution in the second subplot
axs[1].imshow(dx, cmap='copper')
axs[1].set_title('Convolution Result')
axs[1].axis('off')

# Display the subplots
plt.show()