Topic 1 — Part I) Bilinear Interpolation

import numpy as np
import matplotlib.pyplot as plt
import utils_ as ut

ut.set_params()

# Loading the image
image = plt.imread('tiger.jpg')

# Linear Interpolation
def linear_interpolation(y, x, xi):
    i = np.clip(np.searchsorted(x, xi), 1, len(x)-1)
    dx1 = (x[i]-xi) / (x[i]-x[i-1])
    dx2 = (xi-x[i-1]) / (x[i]-x[i-1])
    yi = dx1*y[i-1] + dx2*y[i]
    return yi

# Bilinear Interpolation
def bilinear_interpolation(image, new_shape):
    x_new = np.linspace(0, image.shape[0]-1, new_shape[0])
    y_new = np.linspace(0, image.shape[1]-1, new_shape[1])
    img_x = np.apply_along_axis(linear_interpolation, 0, image, np.arange(image.shape[0]), x_new)
    img_y = np.apply_along_axis(linear_interpolation, 1, img_x, np.arange(image.shape[1]), y_new)
    return img_y

# Applying Bilinear Interpolation
new_image = bilinear_interpolation(image, (1500, 1500))

ut.plt_img(image, 'Original Image')

ut.plt_img(new_image, 'Bilinear Interpolated Image')

Topic 1 — Part II) Histogram Operations

a) Histogram Equalization

# Loading the image and plotting it along its histogram
image = plt.imread('valley.jpg')
ut.present(image, 'Original Image')

# Histogram Equalization
# from https://docs.opencv.org/4.x/d5/daf/tutorial_py_histogram_equalization.html
def histEq(img):
    pdf = np.histogram(img.flatten(), bins=256, range=(0, 256))[0]
    cdf = pdf.cumsum()
    cdf_m = np.ma.masked_equal(cdf, 0)
    cdf_m = (cdf_m-cdf_m.min())*255/(cdf_m.max()-cdf_m.min())
    T = np.ma.filled(cdf_m, 0).astype('uint8')
    pdf_new = np.histogram(T[img].flatten(), bins=256, range=(0, 256))[0]
    return T[img], pdf, pdf_new

# Apply Histogram Equalization
img_histeq = histEq(image)
ut.plt_modh(img_histeq, 'Equalized')

b) Histogram Stretching (Min-Max Normalization)

This is a modified hist-stretch function that implements a general purpose stretching application that clips narrow tails on histograms.

# Histogram Stretching
def histStretch(img, clip=0.0, lamb=255):
    pdf = np.histogram(img.flatten(), bins=256, range=(0, 255))[0]

    ### Implementation of the `autocontrast` in Photoshop ###
    gam = (pdf.sum()*clip/1000)
    # Cusum Control For the Histogram Boundaries
    cdf = pdf.cumsum()
    cdf_r = pdf[::-1].cumsum()[::-1]
    low = np.argmin(np.abs(cdf-gam))-1
    high = np.argmin(np.abs(cdf_r-gam))
    img_ = np.clip(img, low, high)
    # This is necessary because a single 0 or 255 value would
    # prevent the image from being stretched to that limit.
    ###___________________________________________________###

    img_stretch = np.round(((img_-img_.min())/(img_.max()-img_.min())) * lamb)\
        .astype(np.uint8)
    new_pdf = np.histogram(img_stretch.flatten(), bins=256, range=(0, 255))[0]
    return img_stretch, pdf, new_pdf

# Applying Histogram Stretching
img_stretch = histStretch(image)
ut.plt_modh(img_stretch, 'Hist-Stretched')

Histogram stretching normalizes the image to a range of 0 to 255 allowing a better contrast along the spectrum while histogram equalization increases the contrast between the darkest and brightest regions of the image losing the details in the middle of the spectrum.

Some extra stuff: Most common applications of the min-max normalization include a small clipping to better present the contrast (similar to histogram equalization), but without losing details or general lighting effects:

# Histogram Stretching with clipping
img_autocontrast = histStretch(image, clip=4)
ut.plt_modh(img_autocontrast, 'Hist-Stretched with Clipping')

The result matches what Photoshop produces with their 'sophisticated' autocontrast tool. (although their tool has contingency and conditioning for overall dark or bright images that I didn't attempt to implement here.) A better implementation of histogram stretching would also include a hist-mean shift option that would adjust the overall brightness of the image.

Topic 2 — Hybrid Images

import numpy as np
from scipy.ndimage import gaussian_filter
import cv2 as cv
import matplotlib.pyplot as plt
import utils_ as ut

ut.set_params()

# Loading the Images
img_eye = plt.imread('eye.jpg')
img_cam = plt.imread('camera.jpg')

# Converting the Images to Grayscale
img_eye_gray = np.mean(img_eye, axis=2)
img_cam_gray = np.mean(img_cam, axis=2)

def hybridizer(img_eye, img_cam, sigma_eye=1, sigma_cam=1, ratio=0.5):
    # Processing the Eye Image
    img_eye_lowpass = gaussian_filter(img_eye, sigma=sigma_eye)
    # Processing the Camera Image
    img_cam_highpass = img_cam - gaussian_filter(img_cam, sigma=sigma_cam)
    # Creating the Hybrid Image
    img_hybrid = img_eye_lowpass*ratio + img_cam_highpass*(1-ratio)
    # Standardizing the Image
    img_hybrid_norm = cv.normalize(img_hybrid, None, 0, 255, cv.NORM_MINMAX).astype(np.uint8)

    return img_hybrid_norm

# Creating the Hybrid Image
img_hybrid = hybridizer(img_eye_gray, img_cam_gray,
                        sigma_eye=1, sigma_cam=1)

_= plt.imshow(img_hybrid, cmap='gray', vmin=0, vmax=255),\
   plt.axis('off'), plt.title('Hybrid Image'), plt.show()

This hybrid is not as good as one would expect, it appears sigma=1 is not enough blurring to produce a good low or high pass filter, let's explore the effect of varying these parameters.

# Creating 4 more hybrid images with varying sigma values
img_hybrid_1 = hybridizer(img_eye_gray, img_cam_gray,
                          sigma_eye=2, sigma_cam=2)
img_hybrid_2 = hybridizer(img_eye_gray, img_cam_gray,
                          sigma_eye=3, sigma_cam=3)
img_hybrid_3 = hybridizer(img_eye_gray, img_cam_gray,
                          sigma_eye=4, sigma_cam=4)
img_hybrid_4 = hybridizer(img_eye_gray, img_cam_gray,
                          sigma_eye=5, sigma_cam=5)

# plot of the 4 hybrid images in a 2 by 2 grid
fig, ax = plt.subplots(2, 2, figsize=(13, 10))
ax[0, 0].imshow(img_hybrid_1, cmap='gray', vmin=0, vmax=255),\
ax[0, 0].set_title('sigma_eye=2, sigma_cam=2'), ax[0, 0].axis('off')
ax[0, 1].imshow(img_hybrid_2, cmap='gray', vmin=0, vmax=255),\
ax[0, 1].set_title('sigma_eye=3, sigma_cam=3'), ax[0, 1].axis('off')
ax[1, 0].imshow(img_hybrid_3, cmap='gray', vmin=0, vmax=255),\
ax[1, 0].set_title('sigma_eye=4, sigma_cam=4'), ax[1, 0].axis('off')
ax[1, 1].imshow(img_hybrid_4, cmap='gray', vmin=0, vmax=255),\
ax[1, 1].set_title('sigma_eye=5, sigma_cam=5'), ax[1, 1].axis('off')
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()

It appears cam_sigma=4 and eye_sigma=3 are the best parameters to filter the images for a good hybrid:

img_hybrid_best = hybridizer(img_eye_gray, img_cam_gray,
                              sigma_eye=3, sigma_cam=4, ratio=0.6)
_= plt.imshow(img_hybrid_best, cmap='gray', vmin=0, vmax=255), \
plt.axis('off'), plt.title('Hybrid Image'), plt.show()

Up close, the hybrid image clearly shows a camera and the eye is blurred and if you lose your eye focus or look from a distance, only the camera is visible.

Topic 3 — Transformations, Thresholding, and Denoising

import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt
import utils_ as ut

ut.set_params()

a. Reading and displaying the image b. Showing the size of the original image

img = plt.imread('horse.jpg')
print("Image Size: ", img.shape[0],
      "px by ", img.shape[1],"px", sep='')
_= plt.imshow(img), plt.axis('off')

Image Size: 563px by 1000px

...b. Resizing and showing the size of the resized image

# resize image to 1/3 of original size
img_resized = cv.resize(img, (0,0), fx=1/3, fy=1/3)
print("Resized Image Size: ", img_resized.shape[0],
      "px by ", img_resized.shape[1],"px", sep='')

Resized Image Size: 188px by 333px

c. Converting the image to grayscale and also a black and white image with 0.5 level.

img_gray = ut.toGray(img)  #np.dot(img[...,:3],[0.299,0.587,0.114])
img_bw = cv.threshold(img_gray, 256/2, 255, cv.THRESH_BINARY)[1]

_= plt.imshow(img_bw, cmap='gray'), plt.axis('off')

d. Histogram of the grayscale image.

ut.present(img_gray, 'Grayscale Image')

This histogram presents 3 distinct peaks, first peak is likely the horse, the second peak is likely the mountains, and the third peak contains the the sky and highlights in the grass and snow on the mountains. This also indicates that the image can be segmented easily. This histogram is also not spread across the spectrum which means the image could benefit from normalization:

img_stretched = ut.histStretch(img_gray, 3) # reused from Q2
ut.plt_modh(img_stretched, 'Autocontrast')

img_linear = 255 - img_gray
ut.present(img_linear, 'Linear Transformation')

This generally only serves in the presentation of the image, it allows certain light details to be more visible when inverted to dark shades and vise versa.

1. Log-transformation:

img_20_log = np.log(img_gray + 1) * 20
ut.present(img_20_log, r'20$\times$Log Transformation')

img_40_log = np.log(img_gray + 1) * 40
ut.present(img_40_log, r'40$\times$Log Transformation')

Log-transforms are useful for image blending and dynamic range image enhancement, as well as in some high-pass filtering applications.

1. Thresholding:

img_100_thresh = np.where(img_gray > 100, 255, 0)
ut.present(img_100_thresh, 'Thresholding', 15)

img_150_thresh = np.where(img_gray > 150, 255, 0)
ut.present(img_150_thresh, 'Thresholding', 15)

The thresholding is useful for segmenting the image into two parts, separating the foreground and background (given they don't have overlapping shades). It also allows the image to be presented as a binary array useful for high throughput models and ANN applications.

1. Histogram shifting:

img_shift_ = np.clip(img_gray, 25-50, 225-50) #clip between L-S, U-S
img_shift = (img_shift_.astype(int)+50).astype(np.uint8) #add S
ut.present(img_shift, 'Histogram Shifting', 3)

Histogram shifting is useful for adjusting overall brightness of the image, which is one of the steps in gamma correction.

1. Histogram stretching:

img_stretch = ut.histStretch(img_gray, lamb=200)[0]
ut.present(img_stretch, 'Histogram Stretching')

Histogram stretching and shrinking is used in increasing and reducing the overall contrast of the image which is another step in gamma correction.

f. Noise and Noise Reduction:

# Making a noisy image
img_noisy_ = img_gray + np.random.normal(0, scale=10, size=img_gray.shape)
img_noisy = img_noisy_.clip(0, 255).astype(np.uint8)
ut.present(img_noisy, 'Noisy Image')

# Denoising the noisy image using convolution mask
kernel = np.ones(shape=(3,3))*1/9
img_denoised = cv.filter2D(img_noisy, -1, kernel)
ut.present(img_denoised, 'Denoised Image')

# Side by side comparison of the noisy and denoised images
ut.plt_sbs(img_noisy, img_denoised)

g. Sharpening the image:

# Sharpening the original color image using convolution mask
kernel = np.array([[-1, -1, -1],
                   [-1, 9, -1],
                   [-1, -1, -1]])
img_sharpened = cv.filter2D(img, -1, kernel)
_= plt.imshow(img_sharpened), plt.axis('off')

Topic 4 — Edge Detection

import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv
import utils_ as ut

ut.set_params()

a. Loading and convertening the image to grayscale

# Loading the Image
image = plt.imread('lotus-flower.jpg')
img_gray = ut.toGray(image)

b. Using the Otsu's method to find the threshold value

# Otsu's thresholding method
def OtsuThresh(img):
    hist = np.histogram(img, bins=256, range=(0, 255))[0]
    w1, w2 = np.cumsum(hist)+1, np.cumsum(hist[::-1])[::-1]+1
    mu1 = np.cumsum(hist * np.arange(0, 256)) / w1
    mu2 = (np.cumsum((hist * np.arange(0, 256))[::-1]) / w2[::-1])[::-1]
    sigma = w1 * w2 * (mu1 - mu2) ** 2
    filtered = np.where(img<=sigma.argmax(), 0, 255)
    return sigma.argmax(), filtered

# Or simply use the opencv function for identical results:
# thresh,_ = cv.threshold(img_gray, 0, 255, cv.THRESH_BINARY + cv.THRESH_OTSU)

print("Otsu's threshold value:\nmy `OtsuThresh`:",
      OtsuThresh(img_gray)[0], "\n         OpenCV:",
      cv.threshold(img_gray, 0, 255, cv.THRESH_BINARY + cv.THRESH_OTSU)[0])

Otsu's threshold value:
my `OtsuThresh`: 95
         OpenCV: 95.0

c. Image segmentation using k-means clustering

# K-means clustering segmentation of colored image
# Received inspiration from: https://www.programcreek.com/python/example/77061/cv2.KMEANS_RANDOM_CENTERS

def kmSeg(img, k):
    # img_rgb = cv.cvtColor(img, cv.COLOR_BGR2RGB) #used if cv.imread
    img_rgb = img
    img_rgb = img_rgb.reshape((img_rgb.shape[0]*img_rgb.shape[1], 3))
    img_rgb = np.float32(img_rgb)
    criteria = (cv.TERM_CRITERIA_EPS + cv.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    _, label, center = cv.kmeans(img_rgb, k, None, criteria, 10, cv.KMEANS_RANDOM_CENTERS)
    center = np.uint8(center)
    res = center[label.flatten()]
    res = res.reshape((img.shape))
    return res

img_seg_2 = kmSeg(image, 2)
img_seg_3 = kmSeg(image, 3)
img_seg_4 = kmSeg(image, 4)
img_seg_5 = kmSeg(image, 5)

# plotting the 4 segmented images in a 2x2 grid
fig, ax = plt.subplots(2, 2, figsize=(11, 8.5))
fig.suptitle('K-means Segmentation of Image', fontsize=16)
ax[0, 0].imshow(img_seg_2), ax[0, 0].axis('off')
ax[0, 0].set_title('2 clusters')
ax[0, 1].imshow(img_seg_3), ax[0, 1].axis('off')
ax[0, 1].set_title('3 clusters')
ax[1, 0].imshow(img_seg_4), ax[1, 0].axis('off')
ax[1, 0].set_title('4 clusters')
ax[1, 1].imshow(img_seg_5), ax[1, 1].axis('off')
ax[1, 1].set_title('5 clusters'), fig.tight_layout()
plt.show()

d. Edge detection using the Sobel operator

# Sobel operator with some pre/post-processing
def edgeSobel(img, cutoff, k=3, blur=0):
    if blur > 0:
        img = cv.GaussianBlur(img, (blur, blur), 0)
    sobel_x = cv.Sobel(img, cv.CV_64F, 1, 0, ksize=k)
    sobel_y = cv.Sobel(img, cv.CV_64F, 0, 1, ksize=k)
    sobel_ = np.sqrt(sobel_x**2 + sobel_y**2)
    edges_ = np.where(sobel_ > cutoff, 255, 0).astype(np.uint8)
    edges_ = cv.GaussianBlur(edges_, (5, 5), 0)
    edges_ = np.where(edges_>200, 255, 0)
    return edges_.astype(np.uint8), sobel_

img_sob_25 = edgeSobel(img_gray, 25) # Sobel with cutoff=25
_= plt.imshow(img_sob_25[0], cmap='gray'), plt.axis('off')

img_sob_55 = edgeSobel(img_gray, 55) # Sobel with cutoff=25
_= plt.imshow(img_sob_55[0], cmap='gray'), plt.axis('off')

Increasing clipping value removes smaller noise artifacts and smaller edges. Increasing too much will cause the important edges to be lost.

Adding prior bluring will remove noise and edges and reducing the required clipping for a better result:

img_sob_better = edgeSobel(img_gray, 45, blur=5) # cutoff 45, blur 5
_= plt.imshow(img_sob_better[0], cmap='gray'), plt.axis('off')

# Convolves each compass direction to the image
# Returns 8 filtered images.
def compass_filter(img, compass):
    img = cv.GaussianBlur(img, (5,5), 0)
    edges = []
    for i in compass:
        edges.append(cv.filter2D(img, -1, i))

    return np.array(edges)

# Creating the Kernel Compasses
kern_Robins = ut.kernel_compass(np.array([[-1,1,1],[-1,-2,1],[-1,1,1]]))
kern_Krisch = ut.kernel_compass(np.array([[-3,-3,5],[-3,0,5],[-3,-3,5]]))

img_Robins = compass_filter(image, kern_Robins)
ut.plot_compass(img_Robins, "Edge Detection Using Robinson Compass")

img_Kirsch = compass_filter(image, kern_Krisch)
ut.plot_compass(img_Kirsch, "Edge Detection Using Kirsch Compass")

kern_Sobel = ut.kernel_compass(np.array([[-1,0,1],[-2,0,2],[-1,0,1]]))

img_Sobel = compass_filter(image, kern_Sobel)
ut.plot_compass(img_Sobel, "Edge Detection Using Sobel Compass")

So far all results are similar. Other filtering is required to improve the results. Considering Sobel is the most commonly used kernel in other edge detection combination algorithm, its use should be further investigated.

utils_.py

1    import matplotlib.pyplot as plt
2    import seaborn as sns
3    import numpy as np
4    
5    
6    def set_params():
7        sns.set_theme(style="whitegrid", rc={"axes.spines.right": False, "axes.spines.top": False})
8        plt.rcParams['figure.dpi'] = 250
9        plt.rcParams['savefig.dpi'] = 600
10   
11   def plt_img(img, title):
12       figsz = (img.shape[1] / 250, img.shape[0] / 250)
13       plt.figure(figsize=figsz)
14       _ = plt.imshow(img, cmap='gray', vmin=0, vmax=255), \
15           plt.title(title), plt.axis('off'), \
16           plt.figtext(0.5, 0.01, 'Image Size: {}'.format(img.shape),
17                       ha='center'), plt.tight_layout(), plt.show()
18   
19   def toGray(img):
20       img = np.dot(img[..., :3], [0.299, 0.587, 0.114])
21       return np.round(img).astype('uint8')
22   
23   def present(img, title, cols=1.0001):
24       plt.figure(figsize=(10, 11))
25   
26       plt.subplot(211)
27       plt.imshow(img, cmap='gray', vmin=0, vmax=255)
28       plt.axis('off')
29       plt.title(title)
30   
31       plt.subplot(212)
32       sns.histplot(img.ravel(), bins=256, lw=0, color="#05006c", binwidth=cols)
33       plt.xlim([0, 256])
34       plt.title('Histogram')
35       plt.subplots_adjust(wspace=0, hspace=0.1)
36   
37       plt.tight_layout()
38       plt.show()
39   
40   def plt_sbs(img_noisy, img_denoised, title1='Noisy Image', title2='Denoised Image'):
41       plt.figure(figsize=(10, 8))
42   
43       plt.subplot(221)
44       plt.imshow(img_noisy, cmap='gray', vmin=0, vmax=255)
45       plt.axis('off')
46       plt.title(title1)
47   
48       plt.subplot(222)
49       plt.imshow(img_denoised, cmap='gray', vmin=0, vmax=255)
50       plt.axis('off'), plt.title(title2), plt.subplot(212)
51       plt.imshow(np.hstack([img_noisy[300:400, 100:200], img_denoised[300:400, 100:200]]),
52                  cmap='gray', vmin=0, vmax=255)
53       plt.axis('off')
54       plt.title('zoomed')
55   
56       plt.subplots_adjust(wspace=0, hspace=0)
57   
58       plt.tight_layout()
59       plt.show()
60   
61   
62   def plt_modh(img_heq, modification='Equalized'):
63       image, pdf, pdf_new = img_heq
64   
65       plt.subplots(figsize=(8, 10))
66   
67       ax1 = plt.subplot(221)
68       ax1.bar(range(256), pdf, width=1.5, color='#05006c',
69               edgecolor='#05006c', linewidth=0.7)
70       ax1.set_xlim([0, 256])
71       ax1.set_title('Histogram of Original Image')
72   
73       ax2 = plt.subplot(222, sharey=ax1)
74       ax2.bar(range(256), pdf_new, width=1.5, color='#710019',
75               edgecolor='#710019', linewidth=0.7)
76       ax2.set_xlim([0, 256])
77       ax2.set_title('Histogram of {} Image'.format(modification))
78       ax2.tick_params('y', labelleft=False)
79   
80       ax3 = plt.subplot(212)
81       ax3.imshow(image, cmap='gray')
82       ax3.axis('off')
83       ax3.set_title('{} Image'.format(modification))
84   
85       plt.tight_layout()
86       plt.show()
87   
88   # Misc: functions from different Qs to be used across the project.
89   def histStretch(img, clip=0.0, lamb=255):
90       pdf = np.histogram(img.flatten(), bins=256, range=(0, 255))[0]
91   
92       ### Implementation of the `autocontrast` in Photoshop ###
93       gam = (pdf.sum()*clip/1000)
94       # Cusum Control For the Histogram Boundaries
95       cdf = pdf.cumsum()
96       cdf_r = pdf[::-1].cumsum()[::-1]
97       low = np.argmin(np.abs(cdf-gam))-1
98       high = np.argmin(np.abs(cdf_r-gam))
99       img_ = np.clip(img, low, high)
100      # This is necessary because a single 0 or 255 value would
101      # prevent the image from being stretched to that limit.
102      ###___________________________________________________###
103  
104      img_stretch = np.round(((img_-img_.min())/(img_.max()-img_.min())) * lamb)\
105          .astype(np.uint8)
106      new_pdf = np.histogram(img_stretch.flatten(), bins=256, range=(0, 255))[0]
107      return img_stretch, pdf, new_pdf
108  
109  # To create compass kernels
110  # Some parts from https://stackoverflow.com/questions/41502529/
111  def outer_slice(x):
112      return np.r_[x[0],x[1:-1,-1],x[-1,:0:-1],x[-1:0:-1,0]]
113  
114  def rotate_steps(x, shift):
115      out = np.empty_like(x)
116      N = x.shape[0]
117      idx = np.arange(x.size).reshape(x.shape)
118      for n in range((N+1)//2):
119          sliced_idx = outer_slice(idx[n:N-n,n:N-n])
120          out.ravel()[sliced_idx] = np.roll(np.take(x,sliced_idx), shift)
121      return out
122  
123  def kernel_compass(kern_N):
124      compass = [kern_N]
125      for i in range(7):
126          compass.append(rotate_steps(kern_N, i+1))
127      return compass
128  
129  def plot_compass(processed, title):
130      img_plot = np.vstack([np.hstack([processed[1], processed[2], processed[3]]),
131                            np.hstack([processed[0], np.sum(processed / 6, axis=0), processed[4]]),
132                            np.hstack([processed[7], processed[6], processed[5]])]).astype(np.uint8)
133  
134      _= plt.imshow(img_plot), plt.title(title), plt.axis('off')
135