This notebook consists of calculation of integral image, mean grey value of the image, 2D filtering, generating a noisy image, denoising, separability of filters, Fourier transform on the images, and template matching.
Published on December 08, 2021
import cv2 as cv import numpy as np import random import os import sys from numpy.random import randint import matplotlib.pyplot as plt
img_path = "data"
def show_img(img, title, is_gray=False): _ = plt.title(title) if is_gray: _ = plt.imshow(img, cmap="gray") else: _ = plt.imshow(img)
img = cv.imread(os.path.join(img_path, "bonn.png")).astype(np.float32) / 255.0 img = cv.cvtColor(img, cv.COLOR_BGR2RGB) if img is None: print("Could not read the image.") show_img(img, "Original Image")
img_gray = cv.cvtColor(img, cv.COLOR_RGB2GRAY) show_img(img_gray, "Intensity Image", is_gray=True)
Multiplying the intensity image I by 0.5 and subtracting it from each color channel. Make sure that the values do not become negative, i.e. the new (R, G, B) values are (max(R - 0.5I, 0), max(G - 0.5I, 0), max(B - 0.5I, 0)).
Using pixel-wise operations in a nested loop:
def reduce_intensity_with_loops(img, img_gray): new_img = img.copy() for i in range(img.shape[0]): for j in range(img.shape[1]): new_img[i,j] = np.clip(img[i,j] - 0.5 * img_gray[i,j], 0, None) return new_img %timeit reduce_intensity_with_loops(img, img_gray) new_img = reduce_intensity_with_loops(img, img_gray) show_img(new_img, "Reduced Intensity Image", is_gray=False)
3.01 s ± 24 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Using one-line python statement:
def reduce_intensity_one_line(img, img_gray): return np.clip(img - 0.5 * np.expand_dims(img_gray, 2), 0, None) %timeit reduce_intensity_one_line(img, img_gray) new_img = reduce_intensity_one_line(img, img_gray) show_img(new_img, "Reduced Intensity Image", is_gray=False)
1.86 ms ± 48.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Extract a 16 × 16 image patch out of the original image centered at the middle of the image, display it, and copy the content to a random location of the image.
PATCH_SIZE = 16 MAX_ELLIPSE_SIZE = 50 MIN_ELLIPSE_SIZE = 10 MAX_RECT_SIZE = 100 MIN_RECT_SIZE = 20
def extract_patch(img): patch_top, patch_bottom = img.shape[0]//2 - PATCH_SIZE//2 , img.shape[0]//2 + PATCH_SIZE//2 patch_left, patch_right = img.shape[1]//2 - PATCH_SIZE//2 , img.shape[1]//2 + PATCH_SIZE//2 return img[patch_top:patch_bottom, patch_left:patch_right, :] img_patch = extract_patch(img) show_img(img_patch, "Center Patch", is_gray=False)
def place_rand_location(img): rand_coord = randint([img.shape[0]-PATCH_SIZE, img.shape[1]-PATCH_SIZE], size=2) new_img = img.copy() new_img[rand_coord[0]:rand_coord[0]+PATCH_SIZE, rand_coord[1]:rand_coord[1]+PATCH_SIZE] = img_patch return new_img, rand_coord new_img, rand_coord = place_rand_location(img) show_img(new_img, 'Center Patch Placed Random {}, {}' .format(rand_coord[0], rand_coord[1]), is_gray=False)
Draw 10 random rectangles and 10 random ellipses on the image using rectangle and ellipse and display it. Fill the shapes with colors of your choice.
def sample_color(saturation=255, value=255): rand_hue = randint(255,size=1) rand_color = cv.cvtColor(np.uint8([[[rand_hue, saturation, value]]]), cv.COLOR_HSV2BGR )/255 rand_color = tuple(rand_color[0,0].tolist()) return rand_color
def draw_random_rectangles_ellipses(img): new_img = img.copy() for i in range(10): rand_axes = randint(MIN_ELLIPSE_SIZE, MAX_ELLIPSE_SIZE, size=2) rand_center = randint((img.shape[1]-rand_axes[0]*2, img.shape[0]-rand_axes[1]*2), size=2) + rand_axes rand_angle = randint(0,360,size=1) rand_color = sample_color(saturation=127, value=255) new_img = cv.ellipse( new_img, rand_center, rand_axes, float(rand_angle), startAngle=0, endAngle=360, color=rand_color, thickness=-1 ) size = randint(MIN_RECT_SIZE, MAX_RECT_SIZE,size=2) rand_corner1 = randint((img.shape[1]-size[0]*2, img.shape[0]-size[1]*2),size=2) + size rand_corner2 = (rand_corner1[0] + size[0], rand_corner1[1] + size[1]) rand_hue = randint(255,size=1) rand_color = sample_color(saturation=255, value=255) new_img = cv.rectangle(new_img, rand_corner1, rand_corner2, color=rand_color, thickness=-1) return new_img new_img = draw_random_rectangles_ellipses(img) show_img(new_img, "Images with Rectangles and Ellipses", is_gray=False)
/tmp/ipykernel_18922/4251451081.py:3: DeprecationWarning: setting an array element with a sequence. This was supported in some cases where the elements are arrays with a single element. For example `np.array([1, np.array([2])], dtype=int)`. In the future this will raise the same ValueError as `np.array([1, [2]], dtype=int)`. rand_color = cv.cvtColor(np.uint8([[[rand_hue, saturation, value]]]), cv.COLOR_HSV2BGR )/255
def integral_image(img): integral_img = np.zeros((img.shape[0]+1, img.shape[1]+1), dtype=np.float32) for i in range(1, integral_img.shape[0]): for j in range(1, integral_img.shape[1]): integral_img[i,j] = integral_img[i,j-1] + integral_img[i-1,j] - integral_img[i-1,j-1]+img[i-1,j-1] return integral_img integral_img = integral_image(img_gray) show_img(integral_img, "integral_img", is_gray=True)
Computing an integral image using the function integral of opencv library
integral_img = cv.integral(img_gray) show_img(integral_img, "integral_img", is_gray=True)
Average gray value of all pixels in the selection.
In the following task, it can be seen that computing the mean grey value of the image using the integral image is much faster than by summing up each pixel value in the image.
def sum_image(img): return np.sum(img) sum_img = sum_image(img_gray) print(sum_img)
83389.51
def random_crop(img_shape, crop_size): i_start = random.randrange( img_shape[0] - crop_size[0] ) j_start = random.randrange( img_shape[1] - crop_size[1] ) i_stop = i_start + crop_size[0] j_stop = j_start + crop_size[1] return i_start, i_stop, j_start, j_stop
Select 10 random squares of size 100x100 within the image and compute the mean grey value using the two versions below.
1. Using integral image:
def mean_grey_value_by_integ_img(img, integral_img): for i in range(10): i_start, i_stop, j_start, j_stop = random_crop(img.shape, crop_size=(100,100)) pix_sum = integral_img[i_stop, j_stop] - integral_img[i_stop, j_start] - integral_img[i_start, j_stop] + integral_img[i_start, j_start] pix_mean = pix_sum / (100*100) return pix_mean integral_img = integral_image(img_gray) print("Mean grey value: ", mean_grey_value_by_integ_img(img_gray, integral_img))
Mean grey value: 0.3255134521484375
2. By summing up each pixel value in the image
def mean_gray_value_by_summing_pixels(img): for i in range(10): i_start, i_stop, j_start, j_stop = random_crop(img.shape, crop_size=(100,100)) img_cropped = img[i_start:i_stop, j_start:j_stop] pix_mean = np.mean(img_cropped) return pix_mean %timeit mean_gray_value_by_summing_pixels(img_gray) print("Mean grey value: ", mean_gray_value_by_summing_pixels(img_gray))
225 µs ± 6.13 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) Mean grey value: 0.44674042
print("Elapsed time of my implementation:") integral_img = integral_image(img_gray) %timeit mean_grey_value_by_integ_img(img_gray, integral_img) print("\nElapsed time of OpenCV function:") integral_img = cv.integral(img_gray) %timeit mean_grey_value_by_integ_img(img_gray, integral_img)
Elapsed time of my implementation: 47 µs ± 1.09 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) Elapsed time of OpenCV function: 25.3 µs ± 701 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
def gaussian_kernel(sigma, size=(5,5)): kernel = np.zeros(size) x = 0 - ((kernel.shape[0])//2) for i in range(kernel.shape[0]): y = 0 - ((kernel.shape[1])//2) for j in range(kernel.shape[1]): kernel[i,j] = np.exp(-(x**2+y**2)/(2*sigma**2)) y += 1 x += 1 return kernel/np.sum(kernel)
kernel1 = gaussian_kernel(sigma=2, size=(17,17)) kernel2 = gaussian_kernel(sigma=2*np.sqrt(2), size=(17,17)) fig, (axs1, axs2) = plt.subplots(1, 2) _ = axs1.imshow(kernel1, cmap="gray") _ = axs1.set_title("kernel with sigma 2") _ = axs2.imshow(kernel1, cmap="gray") _ = axs2.set_title("kernel with sigma 2*sqrt(2)")
img_padded = np.pad(img_gray, pad_width=(1), mode='constant', constant_values=0) show_img(img_padded, "img with padding", is_gray=True)
Filtering the image twice with a Gaussian kernel with sigma=2
filtered_img = cv.filter2D(img_padded, ddepth=-1, kernel=kernel1) filtered_img = cv.filter2D(filtered_img, ddepth=-1, kernel=kernel1) show_img(filtered_img, "filtered_img with kernel 1", is_gray=True)
Filtering the image once with a Gaussian kernel with sigma=2*sqrt(2)
filtered_img2 = cv.filter2D(img_padded, ddepth=-1, kernel=kernel2) show_img(filtered_img2, "filtered_img with kernel 1", is_gray=True)
Absolute pixel-wise difference between the results
def pixelwise_dif(img1, img2): #np.max(np.abs(np.subtract(img1,img2))) return np.max(np.absolute((img1*255).astype(np.int32) - (img2*255).astype(np.int32))) print("maximum pixel error: ", pixelwise_dif(filtered_img, filtered_img2))
maximum pixel error: 1
def add_salt_n_pepper_noise(img): img_flat = img.copy().flatten() mask = np.random.choice(len(img_flat), int(len(img_flat)*0.3), replace=False) noise = np.random.choice([0.0,1.0], len(mask)) img_flat[mask] = noise return img_flat.reshape(img.shape)
img_noisy = add_salt_n_pepper_noise(img_gray) show_img(img_noisy, 'Image with Salt and Pepper Noise', is_gray=True)
kernel_size = 9 g_kernel = gaussian_kernel(sigma=2, size=(kernel_size,kernel_size)) gaussian_filtered = cv.filter2D(img_noisy, ddepth=-1, kernel=g_kernel) print("maximum pixel error between Gaussian filtered image and original image:\n", pixelwise_dif(gaussian_filtered, img_gray)) median_filtered = cv.medianBlur(np.uint8(img_noisy*255), kernel_size) #, cv.BORDER_CONSTANT print("maximum pixel error between Median filtered image and original image:\n", pixelwise_dif(median_filtered, img_gray)) bilateral_filtered = cv.bilateralFilter(img_noisy, kernel_size, kernel_size * 2, kernel_size / 2) print("maximum pixel error between Bilateral filtered image and original image:\n", pixelwise_dif(bilateral_filtered, img_gray))
maximum pixel error between Gaussian filtered image and original image: 157 maximum pixel error between Median filtered image and original image: 249 maximum pixel error between Bilateral filtered image and original image: 166
fig, axs = plt.subplots(1, 3, figsize=(25,25)) _ = axs[0].imshow(gaussian_filtered, cmap="gray") _ = axs[0].set_title('Gaussian Filtering on Noisy Image') _ = axs[1].imshow(median_filtered, cmap="gray") _ = axs[1].set_title('Median Filtering on Noisy Image') _ = axs[2].imshow(bilateral_filtered, cmap="gray") _ = axs[2].set_title('Bilateral Filtering on Noisy Image')
K1 = np.array([[0.0113, 0.0838, 0.0113], [0.0838, 0.6193, 0.0838], [0.0113, 0.0838, 0.0113]]) K2 = np.array([[-0.8984, 0.1472, 1.1410], [-1.9075, 0.1566, 2.1359], [-0.8659, 0.0573, 1.0337]])
fig, axs = plt.subplots(1, 2) _ = axs[0].imshow(K1); _ = axs[0].set_title("K1") _ = axs[1].imshow(K2); _ = axs[1].set_title("K2")
Filtering the image using the two 2D kernels given above
filtered_img1 = cv.filter2D(img_gray, ddepth=-1, kernel=K1) show_img(filtered_img1, 'Filtered Image with K1', is_gray=True) _ = plt.title('2D Filtered Image with K1') _ = plt.imshow(filtered_img1, cmap="gray", vmin=0, vmax=1)
filtered_img2 = cv.filter2D(img_gray, ddepth=-1, kernel=K2) _ = plt.title('2D Filtered Image with K2') _ = plt.imshow(filtered_img2, cmap="gray", vmin=0, vmax=1)
Using the class SVD of OpenCV to seperate each kernel. If a kernel is not separable, used an approximation by taking only the highest singular value. Filtered the images with the obtained 1D kernels.
# w: calculated singular values # u: left singular vectors # vt: right singular vectors w, u, vt = cv.SVDecomp(K1) vertical_kernel = np.sqrt(w[0]) * u[:, 0] horizontal_kernel = np.sqrt(w[0]) * vt[0] filtered_img3 = cv.sepFilter2D(img_gray, ddepth=-1, kernelX=horizontal_kernel, kernelY=vertical_kernel) _ = plt.title('Filtered Image with SVD using kernel 1') _ = plt.imshow(filtered_img3, cmap="gray", vmin=0, vmax=1)
# w: calculated singular values # u: left singular vectors # vt: right singular vectors w, u, vt = cv.SVDecomp(K2) vertical_kernel = np.sqrt(w[0]) * u[:, 0] horizontal_kernel = np.sqrt(w[0]) * vt[0] filtered_img4 = cv.sepFilter2D(img_gray, ddepth=-1, kernelX=horizontal_kernel, kernelY=vertical_kernel) _ = plt.title('Filtered Image with SVD using kernel 2') _ = plt.imshow(filtered_img4, cmap="gray", vmin=0, vmax=1)
print("maximum pixel error between 2D and SVD versions of kernel 1:\n", pixelwise_dif(filtered_img1, filtered_img3)) print("maximum pixel error between 2D and SVD versions of kernel 2:\n", pixelwise_dif(filtered_img2, filtered_img4))
maximum pixel error between 2D and SVD versions of kernel 1: 1 maximum pixel error between 2D and SVD versions of kernel 2: 30
The convolution in spatial domain can be computed by multiplication in the frequency domain.
image = cv.imread("data/einstein.jpeg", 0) kernel = cv.getGaussianKernel(ksize=7, sigma=1) kernel = kernel * kernel.T show_img(image, "Einstein Image", is_gray=True)
conv_result = cv.filter2D(image, ddepth=-1, kernel=kernel) show_img(conv_result, "conv_result", is_gray=True)
def get_convolution_using_fourier_transform(image, kernel): img_fft2 = np.fft.fft2(image/255) kernel_fft2 = np.fft.fft2(kernel, s=img_fft2.shape) result = np.fft.ifft2(img_fft2 * kernel_fft2) return np.uint8(result.real * 255)
fft_result = get_convolution_using_fourier_transform(image, kernel) show_img(fft_result, "fft_result", is_gray=True)
def imgs_absolute_mean_error(img1, img2): return np.mean(np.absolute(np.int32(img1) - np.int32(img2)))
print("mean absolute difference of the two blurred images:\n", imgs_absolute_mean_error(conv_result, fft_result))
mean absolute difference of the two blurred images: 11.698823529411765
Template matching using normalized cross-correlation as similarity measures.
image = cv.imread("data/lena.png", 0) template = cv.imread("data/eye.png", 0) fig, axs = plt.subplots(1, 2) _ = axs[0].imshow(image, cmap="gray") _ = axs[1].imshow(template, cmap="gray")
def normalized_cross_correlation(image, template): ncc_img = np.zeros(image.shape) normalized_template = template - np.mean(template) for i in range(image.shape[0] - template.shape[0]): for j in range(image.shape[1] - template.shape[1]): patch = image[i:i+template.shape[0], j:j+template.shape[1]] normalized_patch = patch - np.mean(patch) ncc_img[i+template.shape[0]//2 , j+template.shape[1]//2] = np.sum(normalized_template * normalized_patch) / np.sqrt(np.sum(normalized_template**2) * np.sum(normalized_patch**2)) return ncc_img result_ncc = normalized_cross_correlation(image, template) _ = plt.imshow(result_ncc, cmap="gray", vmin=0, vmax=1) #show_img(result_ncc, "result_ncc", is_gray=True)
threshold = 0.7 similarity_img = np.float32(result_ncc>=threshold) _ = plt.imshow(similarity_img, cmap="gray", vmin=0, vmax=1); _ = plt.title("Similarity Img")
indexes = np.where(result_ncc>=threshold) img_colored = cv.cvtColor(image, cv.COLOR_GRAY2BGR) for i,j in zip(indexes[0], indexes[1]): corner1 = (j - template.shape[1]//2, i - template.shape[0]//2) corner2 = (j + template.shape[1]//2, i + template.shape[0]//2) img_colored = cv.rectangle(img_colored, corner1, corner2, color=(255,0,0), thickness=1) _ = plt.imshow(img_colored); _ = plt.title("Eye Detection")