In this chapter, only image enhancement spatial filtering is done.
Spatial filtering is done by modifying all the pixels of an image by a function. There are linear and nonlinear spatial filtering.
The convolution of an image $f(x,y)$ of size $M \times N$ using a kernel $w$ of size $m \times n$ is defined as -
where
$ w(s,t) $ one element of the kernel and
$ f(x+s, y+t) $ is one pixel of the image
If an image $f$ is filterd with a kernel $w_1$, the result filtered with kernel $w_2$, and so on, for Q stages, then because of the cumulative property of convolution, this multistage filtering can be done in a single filtering operation $w*f$, where $w=w_1 * w_2 * w_3 * .... * w_q$
If all the kernels are of size $M \times N$, then the resulted image size will be $W_v\times W_h$ where
$W_v = Q \times (m-1) + m $
$W_h = Q \times (n-1) + n $
# padding of an image
# Zero padding
def my_zero_padding(f, w):
M, N = np.shape(f)
m, n = np.shape(w)
a = int((m-1)/2)
b = int((n-1)/2)
padded_f = np.zeros((M+2*a, N+2*b))
padded_f[a:M+a, b:N+b] = f
return padded_f
# Mirror padding
def my_mirror_padding(f, w):
M, N = np.shape(f)
m, n = np.shape(w)
a = int((m-1)/2)
b = int((n-1)/2)
padded_f = np.zeros((m+M-1, n+N-1))
padded_f[a:M+a, b:N+b] = f
for i in range(b):
padded_f[:, i] = padded_f[:, b+b-i-1]
padded_f[:, N+b+i] = padded_f[:, b+N-i-1]
for j in range(a):
padded_f[j,:] = padded_f[a+a-j-1,:]
padded_f[M+a+j,:] = padded_f[a+M-j-1,:]
return padded_f
# Convolution
def my_convolution2d(f, w, padding):
M, N = np.shape(f)
m, n = np.shape(w)
a = int((m-1)/2)
b = int((n-1)/2)
w_rotated = np.rot90(w, 2) # k = 2 for rotating by 90 degree twice
# padding the original image
if padding == 0:
pf = my_zero_padding(f, w)
elif padding == 1:
pf = my_mirror_padding(f, w)
M_pad, N_pad = np.shape(pf)
convoluted_pf = np.zeros((M_pad, N_pad))
for i in range(a,M+a):
for j in range(b,N+b):
f_st = pf[i-a:i+a+1, j-b:j+b+1] # form the m*n submatrix from the main pf matrix
convoluted_pf[i,j] = np.sum(np.multiply(w_rotated, f_st)) # First element wise multiplication then summation
convoluted_f = convoluted_pf[a:M+a, b:N+b]
return convoluted_f
import numpy as np
# Creating user defined box kernels
def my_box_kernel(size):
box_kernel = np.ones((size))
box_kernel = box_kernel / np.sum(box_kernel) # normalizing
return box_kernel
Gaussian kernel
Here, $\sigma $ is the standard deviation, controls the spread of the Gaussian function
Box filters are poor approximations to the blurring characteristics of lenses.
Box filters favor blurring along perpendicular directions. In applications involving images with a high level of detail, or with strong geometrical components, the directionality of box filters often produces undesirable results
# Creating user defined lowpass Gaussian Filter Kernels
def my_gaussian_kernel(size, K, sigma):
m,n = size
a = int(( m - 1 ) / 2)
b = int(( n - 1 ) / 2)
row = 0
col = 0
w = np.zeros(size)
for s in range (-a, a+1, 1):
for t in range (-b, b+1, 1):
w[row, col] = K * np.exp(-1 * (s**2 + t**2) / (2 * sigma**2))
col = col + 1
row = row + 1
col = 0
w = w / np.sum(w) # Normalizing
return w
my_gaussian_kernel(size = (5, 5), K = 1, sigma = 1)
array([[0.00296902, 0.01330621, 0.02193823, 0.01330621, 0.00296902], [0.01330621, 0.0596343 , 0.09832033, 0.0596343 , 0.01330621], [0.02193823, 0.09832033, 0.16210282, 0.09832033, 0.02193823], [0.01330621, 0.0596343 , 0.09832033, 0.0596343 , 0.01330621], [0.00296902, 0.01330621, 0.02193823, 0.01330621, 0.00296902]])
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0313(a).tif').convert('L'))
#w = my_box_kernel(size=(11,11))
w = my_gaussian_kernel(size = (5,5), K = 1, sigma = 1)
convoluted_image = my_convolution2d(f, w, padding = 1)
plt.imshow(convoluted_image, cmap='gray')
plt.show()
print(np.shape(f))
print(np.shape(convoluted_image))
(500, 500) (500, 500)
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0313(a).tif').convert('L'))
w1 = my_gaussian_kernel(size = (49,49), K = 1, sigma = 8)
convoluted_image_gaussian1 = my_convolution2d(f, w1, padding = 1)
w2 = my_gaussian_kernel(size = (99,99), K = 1, sigma = 8)
convoluted_image_gaussian2 = my_convolution2d(f, w2, padding = 1)
difference = np.absolute(convoluted_image_gaussian1 - convoluted_image_gaussian2)
ax = plt.figure(figsize = (15, 10))
ax1 = ax.add_subplot(131)
ax2 = ax.add_subplot(132)
ax3 = ax.add_subplot(133)
ax1.imshow(convoluted_image_gaussian1 , cmap='gray', vmin=0, vmax=255)
ax2.imshow(convoluted_image_gaussian2, cmap='gray', vmin=0, vmax=255)
ax3.imshow(difference, cmap='gray', vmin=0, vmax=255)
ax1.set_title(r'Kernel size = $49 \times 49$', fontsize = 20)
ax2.set_title(r'Kernel size = $99 \times 99$', fontsize = 20)
ax3.set_title('Difference Image', fontsize = 20)
plt.show()
Here, for $\sigma = 8$, we tried $49 \times 49$ and $99 \times 99 $. The max difference between two resultant image is -
difference.max()
0.992917601085002
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0313(a).tif').convert('L'))
w_box = my_box_kernel(size=(9,9))
w_gaussian = my_gaussian_kernel(size = (19,19), K = 1, sigma = 3)
convoluted_image_box = my_convolution2d(f, w_box, padding = 1)
convoluted_image_gaussian = my_convolution2d(f, w_gaussian, padding = 1)
difference = np.absolute(convoluted_image_box - convoluted_image_gaussian)
ax = plt.figure(figsize = (15, 10))
ax1 = ax.add_subplot(131)
ax2 = ax.add_subplot(132)
ax3 = ax.add_subplot(133)
ax1.imshow(convoluted_image_box, cmap='gray', vmin=0, vmax=255)
ax2.imshow(convoluted_image_gaussian, cmap='gray', vmin=0, vmax=255)
ax3.imshow(difference , cmap='gray', vmin=0, vmax=255)
ax1.set_title('Box Kernel', fontsize = 20)
ax2.set_title('Gaussian Kernel', fontsize = 20)
ax3.set_title('Difference', fontsize = 20)
plt.show()
Comparatively larger Gaussian kernel is necessary to get a close result of a Box kernel. However, box filters favor blurring along perpendicular directions. In applications involving images with a high level of detail, or with strong geometrical components, the directionality of box filters often produces undesirable results.
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0313(a).tif').convert('L'))
w = my_gaussian_kernel(size = (127,127), K = 1, sigma = 21)
convoluted_image_gaussian_zero = my_convolution2d(f, w, padding = 0)
convoluted_image_gaussian_mirror = my_convolution2d(f, w, padding = 1)
ax = plt.figure(figsize = (15, 10))
ax1 = ax.add_subplot(131)
ax2 = ax.add_subplot(132)
ax1.imshow(convoluted_image_gaussian_zero , cmap='gray', vmin=0, vmax=255)
ax2.imshow(convoluted_image_gaussian_mirror, cmap='gray', vmin=0, vmax=255)
ax1.set_title('Zero Padding', fontsize = 20)
ax2.set_title('Mirror Padding', fontsize = 20)
plt.show()
Zero padding introduces dark borders in the filtered image. So, mirror padding is more applicable when the areas near the border contains image details.
Linear filters cannot cope with the nonlinearities of the image formation model and cannot take into account the nonlinearities of human vision. Furthermore, human vision is very sensitive to high-frequency information. Image edges and image details (e.g. corners and lines) have highfrequency content and carry very important information for visual perception. Filters having good edge and image detail preservation properties are highly suitable for digital image filtering. Most of the classical linear digital image filters have low-pass characteristics. They tend to blur edges and to destroy lines, edges, and other fine image details. These reasons have led researchers to the use of nonlinear filtering techniques.
Order Statistic Filters are filters which rank neighboring pixels in an attempt to remove low frequency effects while retaining edges. There are different types of nonlinear filter such as median filter, max filter, min filter, midpoint filter, etc.
Max filters are Useful for finding the brightest points in an image
Mid-Point filters are very useful for removing randomly distributed noise like Gaussian noise
#
def convolution_nonlinear_filter(f, filter_type, filter_size, padding):
M, N = np.shape(f)
w = np.zeros(filter_size)
m, n = filter_size
a = int((m-1)/2)
b = int((n-1)/2)
if padding == 0:
pf = my_zero_padding(f, w)
elif padding == 1:
pf = my_mirror_padding(f, w)
convoluted_pf = np.zeros(np.shape(pf))
for i in range(a, M+a):
for j in range(b, N+b):
f_st = pf[i-a:i+a+1, j-b:j+b+1]
if filter_type == 'median_filter':
convoluted_pf[i,j] = np.median(f_st)
elif filter_type == 'max_filter':
convoluted_pf[i,j] = f_st.max()
elif filter_type == 'min_filter':
convoluted_pf[i,j] = f_st.min()
elif filter_type == 'midpoint_filter':
convoluted_pf[i,j] = int((f_st.min() + f_st.max()) / 2)
convoluted_f = convoluted_pf[a:M+a, b:N+b]
return convoluted_f
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0404(a).tif').convert('L'))
median_filtered_image = convolution_nonlinear_filter(f, filter_type = 'median_filter', filter_size = (9,9), padding=1)
max_filtered_image = convolution_nonlinear_filter(f, filter_type = 'max_filter', filter_size = (9,9), padding=1)
min_filtered_image = convolution_nonlinear_filter(f, filter_type = 'min_filter', filter_size = (9,9), padding=1)
midpoint_filtered_image = convolution_nonlinear_filter(f, filter_type = 'midpoint_filter', filter_size = (9,9), padding=1)
ax = plt.figure(figsize = (15, 10))
ax1 = ax.add_subplot(231)
ax2 = ax.add_subplot(232)
ax3 = ax.add_subplot(233)
ax4 = ax.add_subplot(234)
ax5 = ax.add_subplot(235)
ax1.imshow(f, cmap='gray', vmin=0, vmax=255)
ax2.imshow(median_filtered_image, cmap='gray', vmin=0, vmax=255)
ax3.imshow(max_filtered_image, cmap='gray', vmin=0, vmax=255)
ax4.imshow(min_filtered_image, cmap='gray', vmin=0, vmax=255)
ax5.imshow(midpoint_filtered_image, cmap='gray', vmin=0, vmax=255)
ax1.set_title('Original Image', fontsize = 20)
ax2.set_title('Median Filtered', fontsize = 20)
ax3.set_title('Max Filtered', fontsize = 20)
ax4.set_title('Min Filtered', fontsize = 20)
ax5.set_title('Midpoint Filtered', fontsize = 20)
plt.show()
# Combining all the filters inside a single low_pass_filter function
def my_lowpass_filter(f, lowpass_filter_type, filter_size, padding):
lft = lowpass_filter_type
if lft == 'box':
w_box = my_box_kernel(filter_size)
g_xy = my_convolution2d(f, w_box, padding)
elif lft == 'gaussian':
sigma = int(filter_size.shape[0] / 6)
w_gaussian = my_gaussian_kernel(filter_size, K = 1, sigma = 3)
g_xy = my_convolution2d(f, w_gaussian, padding)
elif lft == 'median_filter' or lft == 'max_filter' or lft == 'min_filter' or lft == 'midpoint_filter':
g_xy = convolution_nonlinear_filter(f, lowpass_filter_type, filter_size, padding)
return g_xy
def my_sharp(f, type, padding):
sharpen_a = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])
sharpen_b = np.array([[1, 1, 1], [1, -8, 1], [1, 1, 1]])
c = -1
M, N = np.shape(f)
m, n = (3, 3)
if padding == 0:
pf = my_zero_padding(f, sharpen_a)
elif padding == 1:
pf = my_mirror_padding(f, sharpen_a)
M_pad, N_pad = np.shape(pf)
convoluted_pf = np.zeros((M_pad, N_pad))
for i in range(1,M+1):
for j in range(1,N+1):
f_st = pf[i-1:i+2, j-1:j+2] # form the m*n submatrix from the main pf matrix
if type == 'a':
convoluted_pf[i,j] = np.sum(np.multiply(sharpen_a, f_st)) # First element wise multiplication then summation
elif type == 'b':
convoluted_pf[i,j] = np.sum(np.multiply(sharpen_b, f_st)) # First element wise multiplication then summation
convoluted_f = convoluted_pf[1:M+1, 1:N+1]
convoluted_f = f + c * convoluted_f
return convoluted_f
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0404(a).tif').convert('L'))
a_shapened_image = my_sharp(f, type = 'a', padding = 1)
b_shapened_image = my_sharp(f, type = 'b', padding = 1)
ax = plt.figure(figsize = (15, 10))
ax1 = ax.add_subplot(231)
ax2 = ax.add_subplot(232)
ax3 = ax.add_subplot(233)
ax1.imshow(f, cmap='gray', vmin=0, vmax=255)
ax2.imshow(a_shapened_image, cmap='gray', vmin=0, vmax=255)
ax3.imshow(b_shapened_image, cmap='gray', vmin=0, vmax=255)
ax1.set_title('Original Image', fontsize = 20)
ax2.set_title('Sharpen a', fontsize = 20)
ax3.set_title('Sharpen b', fontsize = 20)
plt.show()
Subtracting an unsharp (smoothed) version of an image from the original image to sharpen images is called unsharp masking
Steps -
def my_unsharp_masking(f, lowpass_filter_type, filter_size, padding):
f_bar_xy = my_lowpass_filter(f, lowpass_filter_type, filter_size, padding)
g_mask_xy = f - f_bar_xy
g_xy = f + g_mask_xy
return g_xy
def my_highboost_filtering(f, lowpass_filter_type, filter_size, padding, k):
f_bar_xy = my_lowpass_filter(f, lowpass_filter_type, filter_size, padding)
g_mask_xy = f - f_bar_xy
g_xy = f + k * g_mask_xy
return g_xy
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0404(a).tif').convert('L'))
filtered_image_11 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (5,5), padding = 1, k = 1)
filtered_image_12 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (5,5), padding = 1, k = 5)
filtered_image_13 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (5,5), padding = 1, k = 10)
filtered_image_14 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (9,9), padding = 1, k = 1)
filtered_image_15 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (9,9), padding = 1, k = 5)
filtered_image_16 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (9,9), padding = 1, k = 10)
filtered_image_17 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (13,13), padding = 1, k = 1)
filtered_image_18 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (13,13), padding = 1, k = 5)
filtered_image_19 = my_highboost_filtering(f, lowpass_filter_type = 'median_filter', filter_size = (13,13), padding = 1, k = 10)
ax = plt.figure(figsize = (25, 30))
ax1 = ax.add_subplot(431)
ax2 = ax.add_subplot(432)
ax3 = ax.add_subplot(433)
ax4 = ax.add_subplot(434)
ax5 = ax.add_subplot(435)
ax6 = ax.add_subplot(436)
ax7 = ax.add_subplot(437)
ax8 = ax.add_subplot(438)
ax9 = ax.add_subplot(439)
ax1.imshow(f, cmap='gray', vmin=0, vmax=255)
ax2.imshow(filtered_image_11, cmap='gray', vmin=0, vmax=255)
ax3.imshow(filtered_image_12, cmap='gray', vmin=0, vmax=255)
ax4.imshow(filtered_image_13, cmap='gray', vmin=0, vmax=255)
ax5.imshow(filtered_image_14, cmap='gray', vmin=0, vmax=255)
ax6.imshow(filtered_image_15, cmap='gray', vmin=0, vmax=255)
ax7.imshow(filtered_image_16, cmap='gray', vmin=0, vmax=255)
ax8.imshow(filtered_image_17, cmap='gray', vmin=0, vmax=255)
ax9.imshow(filtered_image_18, cmap='gray', vmin=0, vmax=255)
ax1.set_title('Original Image', fontsize = 20)
ax2.set_title(r'$5\times5$ k = 1', fontsize = 20)
ax3.set_title(r'$5\times5$ k = 5', fontsize = 20)
ax4.set_title(r'$5\times5$ k = 10', fontsize = 20)
ax5.set_title(r'$9\times9$ k = 1', fontsize = 20)
ax6.set_title(r'$9\times9$ k = 5', fontsize = 20)
ax7.set_title(r'$9\times9$ k = 10', fontsize = 20)
ax8.set_title(r'$13\times13$ k = 1', fontsize = 20)
ax9.set_title(r'$13\times13$ k = 5', fontsize = 20)
plt.show()
def my_sharp_sobel(f, padding):
w_x = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
w_y = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
M, N = np.shape(f)
if padding == 0:
pf = my_zero_padding(f, w_x)
elif padding == 1:
pf = my_mirror_padding(f, w_x)
M_pad, N_pad = np.shape(pf)
g_x = np.zeros((M_pad, N_pad))
g_y = g_x
for i in range(1,M+1):
for j in range(1,N+1):
f_st = pf[i-1:i+2, j-1:j+2]
g_x[i,j] = np.sum(np.multiply(w_x, f_st))
g_y[i,j] = np.sum(np.multiply(w_y, f_st))
filtered = np.add(np.absolute(g_x), np.absolute(g_y))
M_xy = filtered[1:M+1, 1:N+1]
# M_xy may be greater than 255, so we need to scale it to [0, 255]
M_xy_scaled = int(M_xy * 255 / M_xy.max)
return M_xy_scaled
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
f = np.array(Image.open('Fig0404(a).tif').convert('L'))
b = my_sharp_sobel(f, padding = 1)
ax = plt.figure(figsize = (25, 10))
ax1 = ax.add_subplot(431)
ax2 = ax.add_subplot(432)
ax1.imshow(f, cmap='gray', vmin=0, vmax=255)
ax2.imshow(b, cmap='gray', vmin=0, vmax=255)