Gaussian Filter Implementation

Gino Prasad

09/20/2022

import cv2
from pylab import imshow
import numpy as np

image_path = "/Users/ginoprasad/image_processing/gaussian blur test image.png"

img = cv2.imread(image_path)

img.shape

(275, 491, 3)

imshow(img)

<matplotlib.image.AxesImage at 0x7f98406ebee0>

fact_cache = [1]
def factorial(n):
    while len(fact_cache) <= n:
        fact_cache.append(fact_cache[-1] * len(fact_cache))
    return fact_cache[n]

sin = 0, 1, 0, -1 0 if even, 1 if odd, 1 if bitshift is even, -1 if bitshift is odd

def sin(x, num_iterations=100):
    sin_nth_deriv = lambda n_: (n_ & 1) * (1 - (((n_ >> 1) & 1) << 1))
    return sum(sin_nth_deriv(n) * (x ** n)  / factorial(n) for n in range(num_iterations))

def cos(x, num_iterations=100):
    sin_nth_deriv = lambda n_: ((n_ + 1) & 1) * (1 - (((n_ >> 1) & 1) << 1))
    return sum(sin_nth_deriv(n) * (x ** n)  / factorial(n) for n in range(num_iterations))

e derivative is 1 always

def e(x, num_iterations=100):
    return sum((x ** n)  / factorial(n) for n in range(num_iterations))

def newton_raphson(function, derivative, x=0.1, num_iterations=100):
    for _ in range(num_iterations):
        x -= function(x) / derivative(x)
    return x

def sqrt(a, num_iterations=100):
    x = 0.1
    function = lambda x: (x * x) - a
    derivative = lambda x: 2 * x
    return newton_raphson(function, derivative, num_iterations=num_iterations)

pi = newton_raphson(sin, cos, x=3.14)

def G(x, y, sigma=1):
    sigma_squared = sigma * sigma
    return (1 / (2 * pi * sigma_squared)) * e(-((x * x) + (y * y)) / (2 * sigma_squared))

def get_gaussian_kernel(kernel_shape, sigma=1):
    kernel_matrix = np.zeros(kernel_shape)
    width, height = kernel_matrix.shape

    for i in range(width):
        for j in range(height):
            x = i - (width // 2)
            y = j - (height // 2)
            kernel_matrix[i, j] = G(x, y, sigma=sigma)

    kernel_matrix /= sum(sum(kernel_matrix))
    return kernel_matrix

imshow(get_gaussian_kernel((5, 5), 1))

<matplotlib.image.AxesImage at 0x7f9810075df0>

def discrete_convolve(f, g):
    width, height = f.shape
    width_max, height_max = g.shape
    ret = np.zeros(g.shape)
    for i in range(width_max):
        for j in range(height_max):
            if i-(width // 2) >= 0 and i+(width//2) < width_max and j-(height // 2) >= 0 and j+(height//2) < height_max:
                spliced_g = g[i-(width // 2):i+(width//2)+1,j-(height // 2):j+(height//2)+1]
                try:
                    assert spliced_g.shape == f.shape
                except AssertionError as e:
                    print(spliced_g.shape)
                    raise e
                ret[i, j] = sum(sum(spliced_g * f))
    return ret

def discrete_convolve_image(kernel, img):
    ret = []
    for channel in range(3):
        channel_img = img[:,:,channel]
        ret.append(discrete_convolve(kernel, channel_img))
    return np.dstack(ret).astype(int)

def gaussian_blur(kernel_shape=(21,21), sigma=3):
    kernel = get_gaussian_kernel(kernel_shape, sigma=sigma)
    return discrete_convolve_image(kernel, img)

imshow(img)

<matplotlib.image.AxesImage at 0x7f98100a4a60>

imshow(gaussian_blur(sigma=3))

<matplotlib.image.AxesImage at 0x7f98416b32b0>

imshow(gaussian_blur(sigma=20))

<matplotlib.image.AxesImage at 0x7f984170bd00>

print(newton_raphson(sin, cos, x=3.14))

3.141592653589794

kernel_matrix = np.zeros((5, 5))
width, height = kernel_matrix.shape

for i in range(width):
    for j in range(height):
        x = i - (width // 2)
        y = j - (height // 2)
        kernel_matrix[i, j] = G(x, y, sigma=1)

kernel_matrix /= sum(sum(kernel_matrix))

sum(sum(kernel_matrix))

0.9999999999999999

kernel_matrix[0:2, 0:1]

array([[0.00296902],
       [0.01330621]])

kernel_matrix

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]])

imshow(kernel_matrix)

<matplotlib.image.AxesImage at 0x7f9831772700>