Experments with OpenCV and Sk-Image¶
The document covers some of the most commonly known image processing and analysis functionalities.
In [1]:
import gdal
import numpy as np
from osgeo import gdal
import pandas as pd
# For visualization
% matplotlib inline
from matplotlib import pyplot as plt
import cv2
import skimage
from skimage import io
from skimage import external
import rasterio
Image¶
In [2]:
img_gray = cv2.imread('high_res_img.tif',0)
#subset_area
img_a= img_gray[3750:4000, 1750:2000]
Histogram Equalization¶
Histogram equalization is a technique of enhancing image contrast by spreading out the most frequent pixel values. There are two kinds:
- Global Histogram Equalization (GHE):- Increases the over all contrast of the entire image
- Local Histogram equalization(LHE):- emphasizes local graylevel variations
In [5]:
# Load modules
from skimage.exposure import equalize_hist
from skimage.morphology import disk
from skimage.filters import rank
#Set plot dimension
plt.subplots(2, 3, figsize=(20,16))
#GHE
global_hist_eq = equalize_hist(img_a)
#LHE
selem = disk(30)
local_hist_eq = rank.equalize(img_a, selem=selem)
#plot
plt.subplot(231),plt.imshow(img_a, cmap='gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(234),plt.hist(img_a), plt.ylabel('Frequency'), plt.xlabel('Pixel Values')
plt.subplot(232), plt.imshow(global_hist_eq, cmap='gray')
plt.title('Global Equalization'), plt.xticks([]), plt.yticks([])
plt.subplot(235), plt.hist(global_hist_eq), plt.ylabel('Frequency'), plt.xlabel('Pixel Values')
plt.subplot(233), plt.imshow(local_hist_eq, cmap='gray')
plt.title('Local Equalization'), plt.xticks([]), plt.yticks([])
plt.subplot(236), plt.hist(local_hist_eq), plt.ylabel('Frequency'), plt.xlabel('Pixel Values')
plt.show()
Noise Removal¶
In [6]:
plt.subplots(1, 2, figsize=(20,12))
img_g = skimage.filters.gaussian(img_a)
plt.subplot(121),plt.imshow(img_a, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.title('Original')
plt.subplot(122),plt.imshow(img_g, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.title('Gaussian')
plt.show()
Canny algorithm¶
In [7]:
plt.subplots(1, 2, figsize=(20,12))
edges =cv2.Canny(img_a, 0,0.4,0.5)
plt.subplot(121), plt.imshow(img_a, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.title('Original')
plt.subplot(122), plt.imshow(edges), plt.xticks([]), plt.yticks([])
plt.title('Canny')
plt.show()
Sobel X&Y, Laplacian Filters¶
In [8]:
plt.subplots(2, 2, figsize=(18,18))
laplacian = cv2.Laplacian(img_g, cv2.CV_64F)
sobelx = cv2.Sobel(img_g, cv2.CV_64F, 1,0, ksize=5)
sobely = cv2.Sobel(img_g, cv2.CV_64F, 0,1, ksize=5)
plt.subplot(2,2,1),plt.imshow(img_a,cmap = 'gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(2,2,2),plt.imshow(laplacian,cmap = 'Blues')
plt.title('Laplacian'), plt.xticks([]), plt.yticks([])
plt.subplot(2,2,3),plt.imshow(sobelx,cmap = 'Blues')
plt.title('Sobel X'), plt.xticks([]), plt.yticks([])
plt.subplot(2,2,4),plt.imshow(sobely,cmap = 'Blues')
plt.title('Sobel Y'), plt.xticks([]), plt.yticks([])
plt.show()
Fourier Transform¶
In [27]:
plt.subplots(1, 2, figsize=(20,12))
rows, cols = img_a.shape
crow,ccol = rows/2 , cols/2
#fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
f = np.fft.fft2(img_a)
fshift = np.fft.fftshift(f)
fshift[125-30:125+30, 125-30:125+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
plt.subplot(121),plt.imshow(img_a, cmap = 'gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap="Reds")
plt.title('Fourier Transform'), plt.xticks([]), plt.yticks([])
plt.show()
Histogram of Oriented Gradients (HOG)¶
In [34]:
plt.subplots(1, 2, figsize=(20,12))
from skimage.feature import hog
d, hog_img = hog(img_a, orientations=9, pixels_per_cell=(2,2),
cells_per_block=(2, 2),block_norm='L1', visualize=True, multichannel=False)
plt.subplot(121),plt.imshow(img_a, cmap='gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(hog_img, cmap='Greys')
plt.title('HOG'), plt.xticks([]), plt.yticks([])
plt.show()
Local Binary Patterns(LBP)¶
In [48]:
plt.subplots(1, 2, figsize=(20,12))
from skimage.feature import local_binary_pattern
radius = 1
n_points = 16 * radius
METHOD = 'uniform'
plt.rcParams['font.size'] = 9
lbp_img = local_binary_pattern(img_a, n_points, radius, METHOD)
plt.subplot(121),plt.imshow(img_a, cmap='gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(lbp_img, cmap='gray')
plt.title('LBP'), plt.xticks([]), plt.yticks([])
plt.show()
Local Maxima¶
In [51]:
from scipy import ndimage as ndi
from skimage.feature import peak_local_max
from skimage import data, img_as_float
image_max = ndi.maximum_filter(img_a, size=20, mode='constant')
coordinates = peak_local_max(img_a, min_distance=20)
fig, axes = plt.subplots(1, 3, figsize=(16, 12), sharex=True, sharey=True)
ax = axes.ravel()
ax[0].imshow(img_a, cmap=plt.cm.gray)
ax[0].axis('off')
ax[0].set_title('Original')
ax[1].imshow(image_max, cmap=plt.cm.Greys)
ax[1].axis('off')
ax[1].set_title('Maximum filter')
ax[2].imshow(img_a, cmap=plt.cm.gray)
ax[2].autoscale(False)
ax[2].plot(coordinates[:, 1], coordinates[:, 0], 'r.')
ax[2].axis('off')
ax[2].set_title('Peak local max')
fig.tight_layout()
plt.show()
Contours with thresholds¶
In [52]:
plt.subplots(1, 3, figsize=(20,12))
ret,thresh = cv2.threshold(img_a, 127,255,0)
image, contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
cntr_img = cv2.drawContours(img_a, contours, -1, (0,255,0), 3)
ret_lhe,thresh_lhe = cv2.threshold(local_hist_eq, 125,255,0)
image, contours, hierarchy = cv2.findContours(thresh_lhe,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
cntr_img_lhe =cv2.drawContours(local_hist_eq, contours, -1, (125,155,0),3)
plt.subplot(131),plt.imshow(img_a, cmap='gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(cntr_img, cmap='Blues')
plt.title('Contour of Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(cntr_img_lhe, cmap='Blues')
plt.title('Contour for LHE Image'), plt.xticks([]), plt.yticks([])
plt.show()
Denoising the image¶
Deonising an image is a technique of smoothening interior of grains, while keeping sharp edges between grains.
In [53]:
plt.subplots(1,2, figsize=(20,16))
from skimage import restoration
from skimage import img_as_float
img_float = img_as_float(img_a)
img_denoised= restoration.denoise_nl_means(img_float, h=0.001, multichannel=False)
plt.subplot(121),plt.imshow(img_a,cmap='Blues'),plt.xticks([]), plt.yticks([])
plt.title('Orignal')
plt.subplot(122),plt.imshow(img_denoised,cmap='Blues'), plt.xticks([]),plt.yticks([])
plt.title('Denoised')
plt.show()
Matplotlib Contours¶
In [65]:
plt.subplots(1,2, figsize=(20,16))
# Try to threshold the image
plt.subplot(121),plt.imshow(img_a, cmap='gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_a, cmap='gray')
plt.contour(img_a, levels=[0], colors='yellow')
plt.title('Contours'), plt.xticks([]), plt.yticks([])
plt.show()
Canny Algorithm For Edge Detection¶
In [78]:
plt.subplots(1,2, figsize=(20,16))
# Orignal
from skimage import feature
plt.subplot(121), plt.imshow(img_a, cmap='gray')
plt.title('Original'), plt.xticks([]), plt.yticks([])
# Denoised
edges_d = feature.canny(img_a, sigma=0.1,
low_threshold=0.001,
high_threshold=0.009)
plt.subplot(122), plt.imshow(img_denoised, cmap='gray')
plt.contour(edges_d), plt.title('Canny Detected Edge')
plt.xticks([]), plt.yticks([])
plt.show()
Top-Hat Transform¶
Top-hat transform is an operation that extracts small elements and details from given images.
In [82]:
plt.subplots(1,2, figsize=(20,16))
from scipy import ndimage
from skimage import morphology
wat_img_a= ndimage.black_tophat(img_float, 10)
wat_img_a -= 0.3* img_float
wat_img_a= morphology.dilation(wat_img_a)
plt.subplot(121),plt.imshow(wat_img_a, cmap='gray')
plt.title('Original')
wat_img= ndimage.black_tophat(img_denoised, 10)
wat_img -= 0.3* img_denoised
wat_img= morphology.dilation(wat_img)
plt.subplot(122),plt.imshow(wat_img, cmap='gray')
plt.title('Top Hat Transformed')
plt.show()
