[PYTHON] Digital image processing (spatial filtering)

* Spatial filtering *

It is an area to pixel operation. A filter is applied to each pixel, and a product-sum operation is performed to smooth the image (noise removal) and detect edges. For the image, apply the filter at high speed like this video, and input the calculated value obtained from the filter to the case pixel. (Also called convolution)

Smoothing (reducing the unevenness of pixel values)

* Image to use *

1_av_fil_before.jpg

Library used

python


import numpy as np
from numpy.lib.stride_tricks import as_strided
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import matplotlib.gridspec as gridspec
import seaborn as sns
import math
import cv2

* Average filter *

The average of the pixel values in the region is taken as the value of the pixel of interest.

python


av_filter = np.array([[ 1/25,1/25,1/25,1/25,1/25,],
                      [ 1/25,1/25,1/25,1/25,1/25,],
                      [ 1/25,1/25,1/25,1/25,1/25,],
                      [ 1/25,1/25,1/25,1/25,1/25,],
                      [ 1/25,1/25,1/25,1/25,1/25,]], np.float32)

# ddepth = -1, means destination image has depth same as input image
dst1 = cv2.filter2D(image_1, -1, av_filter)
cv2.imwrite('2_av_fil.jpg', dst1)

↓ It is smoothed as a whole.

2_av_fil.jpg

* Weighted average filter *

Also called * weighted averaging *, a distribution centered on the pixel of interest, with the weight decreasing toward the outside.

python


add_av_fil = np.array([[1, 4, 6, 4, 1],
       [4, 16, 24, 16, 4],
       [6, 24, 36, 24, 6],
       [4, 16, 24, 16, 4],
       [1, 4, 6, 4, 1]], np.float32)/256

dst1 = cv2.filter2D(res, -1, add_av_fil)
cv2.imwrite('3_add_av_fil.jpg', dst1)

↓ It is relatively smoothed while preserving the edges.

3_add_av_fil.jpg

* Gaussian filter *

Weighted according to Gaussian distribution.

h\left(x, y\right) = \frac{1}{2\pi\sigma^2}{\exp{\left(-\frac{x^2+y^2}{2\sigma^2}\right)}}\\

python


def gausian(X, y, siguma):
#     exp(-(x2 + y2)/2σ2)
    h2 = math.e**(-1*((X**2) + (y**2))/(2*(siguma**2)))
#     2πσ2
    h3 = 2*math.pi*(siguma**2)
    h = h2/h3
    return h

sns.set_style('ticks')

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection='3d')

x = y = np.arange(-3,3,0.1)
X,Y = np.meshgrid(x,y)

ax.plot_surface(X,Y, gausian(X,Y,1),cmap='coolwarm')
plt.show()
plt.savefig("gausian_fig.jpg ")

↓ Such an image

gausian_fig.jpg

Handmade

python


def gausian_kernel(size):
    if size%2==0:
        print('kernel size should be odd')
        return
    sigma = (size-1)/2
    
#     [0,size]→[-sigma, sigma]Shift
    x = y = np.arange(0,size) - sigma
    X,Y = np.meshgrid(x,y) 
    print(sigma)
    
    mat = gausian(X,Y,sigma)
    
    #So that the sum is 1
    kernel = mat / np.sum(mat)
    return mat

g_fil_hand = gausian_kernel(3)

#Filter shape
#[[0.05854983 0.09653235 0.05854983]
#[0.09653235 0.15915494 0.09653235]
#[0.05854983 0.09653235 0.05854983]]

dst1 = cv2.filter2D(res, -1, g_fil_hand)
cv2.imwrite('4_g_hand_fil.jpg', dst1)

4_g_hand_fil.jpg

by opencv

python


sigma = 2
blur = cv2.GaussianBlur(res, (0, 0), sigmaX = 3, sigmaY = 3)
cv2.imwrite('5_g_cv_fil.jpg', blur)

5_g_cv_fil.jpg

* Smoothing in a specific direction *

By using the following filters, smoothing is performed only in a specific direction. The picture looks like it's shaken sideways.

python


flat_d_fil = np.zeros([15, 15])

for i in range(15):
    flat_d_fil[i, i] = 1/15

dst1 = cv2.filter2D(res, -1, flat_d_fil)
cv2.imwrite('6_shaky_fil.jpg', dst1)

6_shaky_fil.jpg

* Derivative filter *

For vertical edges

python


diff_col = np.array([[0, 0, 0],
        [1, 0, -1],
        [0, 0, 0]])

dst1 = cv2.filter2D(res, -1, 8*diff_col)
cv2.imwrite('7_diff_col.jpg', dst1)

7_diff_col.jpg

For horizontal edges

python


diff_row = np.array([[0, -1, 0],
        [0, 0, 0],
        [0, 1, 0]])

dst1 = cv2.filter2D(res, -1, 8*diff_row)
cv2.imwrite('8_diff_row.jpg', dst1)

8_diff_row.jpg

* Pruwit filter *

Combination of differential filter and average filter Pick up the gradient of the pixel value in the horizontal direction + Smooth the gradient of the pixel value in the vertical direction

python


pulu_fil = np.array([[1, 0, -1],
        [1, 0, -1],
        [1, 0, -1]])

dst1 = cv2.filter2D(res, -1, 2*pulu_fil)
cv2.imwrite('9_pulu_fil.jpg', dst1)

9_pulu_fil.jpg

* Sobel filter *

Combination of differential filter and weighted average filter Achieves a smoother degree of blurring than Pruwit. Pick up the gradient of the pixel value in the horizontal direction + Smooth the gradient of the pixel value in the vertical direction

python


sobel_fil = np.array([[1, 0, -1],
        [2, 0, -2],
        [1, 0, -1]])

dst1 = cv2.filter2D(res, -1, sobel_fil)
cv2.imwrite('10_sobel_fil.jpg', dst1)

10_sobel_fil.jpg

* Second derivative filter *

Feeling that the differential filter is doubled

python


todiff_fil = np.array([[0, 0, 0],
                      [1, -2, 1],
                      [0, 0, 0]])

dst1 = cv2.filter2D(res, -1, 2*todiff_fil)
cv2.imwrite('11_to_diff_fil.jpg', dst1)

* Laplacian filter *

It feels like the vertical + horizontal of the second derivative.

python


lap_fil = np.array([[0, 1, 0],
        [1, -4, 1],
        [0, 1, 0]])

dst1 = cv2.filter2D(res,ddepth=cv2.CV_16S,kernel=lap_fil)
cv2.imwrite('12_lap_hand_fil.jpg', dst1+100)

plt.plot(ds[400])
plt.savefig("lap_fig_1.jpg ")

12_lap_fil.jpg

lap_fig_1.jpg

with opencv

python


ddepth = cv2.CV_16S
kernel_size = 3
dst = cv2.Laplacian(res, ddepth, ksize=kernel_size)
cv2.imwrite('13_lap_cv_fil.jpg', dst)

plt.plot(dst[400]+100)
plt.savefig("lap_fig_2.jpg ")

13_lap_cv_fil.jpg

* Log filter (Laplacian Gaussian) *

h\left(x, y\right) = \frac{x^2+y^2-2\sigma^2}{2\pi\sigma^6}{\exp{\left(-\frac{x^2+y^2}{2\sigma^2}\right)}}\\

Handmade

python


def Log_fil(X, y, siguma):
    
#     X2+y2-2σ2
    h1 = ((X**2) + (y**2) - (2*(siguma**2))) 
    
#     exp(-(x2 + y2)/2σ2)
    h2 = math.e**(-1*((X**2) + (y**2))/(2*(siguma**2)))
    
#     2πσ6
    h3 = 2*math.pi*(siguma**6)

    h = h1*h2/h3
    return h

sns.set_style('ticks')

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection='3d')

x = y = np.arange(-3,3,0.1)
X,Y = np.meshgrid(x,y)

ax.plot_surface(X,Y,Log_fil(X,Y,1),cmap='coolwarm')
plt.savefig("log_fig.jpg ")
plt.show()

log_fig.jpg

python


def L_of_g_kernel(size, sigma):
    if size%2==0:
        print('kernel size should be odd')
        return
    ad_size = (size-1)/2
    
    # [0,size]→[-sigma, sigma]Shift
    x = y = np.arange(0,size) - ad_size
    X,Y = np.meshgrid(x,y) 
    print(sigma)
    
    mat = Log_fil(X,Y,sigma)
    print(mat)
    
    #So that the sum is 0
    kernel = mat - np.sum(mat)/size**2
    
    return kernel

log_fil =  L_of_g_kernel(5, 1)

#Filter shape
#[[ 0.01749015  0.0391927   0.04307856  0.0391927   0.01749015]
# [ 0.0391927   0.         -0.09653235  0.          0.0391927 ]
# [ 0.04307856 -0.09653235 -0.31830989 -0.09653235  0.04307856]
# [ 0.0391927   0.         -0.09653235  0.          0.0391927 ]
# [ 0.01749015  0.0391927   0.04307856  0.0391927   0.01749015]]

python


dst1 = cv2.filter2D(res,ddepth=cv2.CV_16S,kernel=log_fil)
cv2.imwrite('14_log_fil_hand.jpg', dst1*20+200)
plt.savefig("log_fill_1.jpg ")
plt.plot(dst1[400]*20+200)

14_log_fil_hand.jpg

log_fill_1.jpg

by opencv

python


ddepth = cv2.CV_16S
kernel_size = 5

# Apply Gaussian Blur
blur = cv2.GaussianBlur(res,(5, 5), sigmaX = 1, sigmaY = 1)

# Apply Laplacian operator in some higher datatype
dst = cv2.Laplacian(blur, ddepth, ksize=kernel_size)
cv2.imwrite('15_log_fil_cv.jpg', dst+200)

plt.plot(dst[400])
plt.savefig("log_fill_2.jpg ")

15_log_fil_cv.jpg

log_fill_2.jpg

Extraction of 0 intersections

python


def Zero_crossing(image, thresh):
    z_c_image = np.zeros(image.shape)
    
    # For each pixel, count the number of positive
    # and negative pixels in the neighborhood
    
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            negative_count = 0
            positive_count = 0
            neighbour = [image[i+1, j-1],image[i+1, j],image[i+1, j+1],image[i, j-1],image[i, j+1],image[i-1, j-1],image[i-1, j],image[i-1, j+1]]
            d = max(neighbour)
            e = min(neighbour)
            for h in neighbour:
                if h>0:
                    positive_count += 1
                elif h<0:
                    negative_count += 1
 
 
            # If both negative and positive values exist in 
            # the pixel neighborhood, then that pixel is a 
            # potential zero crossing
            
            z_c = ((negative_count > 0) and (positive_count > 0) and (d-e > thresh))
            
            # Change the pixel value with the maximum neighborhood
            # difference with the pixel
 
            if z_c:
                z_c_image[i, j] = 200

 
    return z_c_image

dst2 = Zero_crossing(dst1, 10)
cv2.imwrite('16_log_fil_zero.jpg', dst2)

16_log_fil_zero.jpg

* Detection by Canny *

1, smoothing

Smoothing with Gaussian filter

2, Processing with vertical and horizontal differential filters

You are using the sobel filter.

3, Calculation of gradient amount and gradient direction

M(x,y) = \sqrt{G_x^2+G_y^2}
\Theta = atan2 \left(G_y/G_x \right)

4, Threshold processing

Only a certain range is detected in order to remove false positives.

python


class Canny_hand:
    def __init__(self, image, sigma, h_thresh, l_thresh):
        self.image = image.astype(np.float32)#Data type specification
        self.h, self.w = image.shape
        self.sigma = sigma
        self.h_thresh = h_thresh
        self.l_thresh = l_thresh
        self.sobelx = np.array([[-1, 0, 1],
                          [-2, 0, 2],
                          [-1, 0, 1]])

        self.sobely = np.array([[-1, -2, -1],
                          [0,  0, 0],
                          [1, 2, 1]])

    def do(self,):
        blur = cv2.GaussianBlur(self.image,(0, 0),  sigmaX = self.sigma, sigmaY = self.sigma)
        Gx = cv2.filter2D(blur, -1, self.sobelx)
        Gy = cv2.filter2D(blur, -1, self.sobely)
        G_value = np.sqrt(np.square(Gy) + np.square(Gx)) 
        G_angle = np.arctan2(Gy, Gx)
        G_angle_set = self.set_angle(G_angle)
        G_max_group = self.extract_max(G_value, G_angle_set)
        result = self.thresh_hold(G_max_group, min=self.l_thresh, max=self.h_thresh)

        return result

    def set_angle(self, G_angle):
        pai = math.pi
        
        G_angle[np.where((G_angle >= -pai/8) & (G_angle < pai/8))] = 0
        G_angle[np.where((G_angle >= pai/8) & (G_angle < 3*pai/8))] = 45
        G_angle[np.where((G_angle >= 3*pai/8) & (G_angle < 5*pai/8))] = 90
        G_angle[np.where((G_angle >= 5*pai/8) & (G_angle < 7*pai/8))] = 135
        G_angle[np.where((G_angle >= 7*pai/8) & (G_angle < -7*pai/8))] = 180
        G_angle[np.where((G_angle >= -7*pai/8) & (G_angle < -5*pai/8))] = 225
        G_angle[np.where((G_angle >= -5*pai/8) & (G_angle < -3*pai/8))] = 270
        G_angle[np.where((G_angle >= -3*pai/8) & (G_angle < -pai/8))] = 315

        return G_angle

    def extract_max(self, G_value, G_angle_set):
        result = G_value.copy()

        for y in range(1, self.h - 1):
            for x in range(1, self.w - 1):

                if G_angle_set[y][x] == 0:
                    if (G_value[y][x] < G_value[y][x+1]):
                        result[y][x] = 0
                elif G_angle_set[y][x] == 45:
                    if (G_value[y][x] < G_value[y+1][x+1]):
                        result[y][x] = 0
                elif G_angle_set[y][x] == 90:
                    if (G_value[y][x] < G_value[y+1][x]):
                        result[y][x] = 0
                elif G_angle_set[y][x] == 135:
                    if (G_value[y][x] < G_value[y+1][x-1]):
                        result[y][x] = 0
                elif G_angle_set[y][x] == 180:
                    if (G_value[y][x] < G_value[y][x-1]):
                        result[y][x] = 0
                elif G_angle_set[y][x] == 225:
                    if (G_value[y][x] < G_value[y-1][x-1]):
                        result[y][x] = 0
                elif G_angle_set[y][x] == 270:
                    if (G_value[y][x] < G_value[y-1][x]):
                        result[y][x] = 0
                elif G_angle_set[y][x] == 315:
                    if (G_value[y][x] < G_value[y-1][x+1]):
                        result[y][x] = 0

        return result

    def thresh_hold(self, in_put, min=75, max=150, d=1):

        for y in range(0, self.h):
            for x in range(0, self.w):

                if in_put[y][x] >= max:
                    in_put[y][x] = 255

                elif in_put[y][x] < min:
                    in_put[y][x] = 0

                else:
                    if np.max(in_put[y-d:y+d+1, x-d:x+d+1]) >= max:
                        in_put[y][x] = 255
                    else:
                        in_put[y][x] = 0

        return in_put

canny = Canny_hand(res, 0.3,100, 200 )
can_do = canny.do()
cv2.imwrite('17_canny_hand.jpg', can_do)

17_canny_hand.jpg

by opencv

python


edges = cv2.Canny(res,0.3,  100, 200)
cv2.imwrite('17_canny.jpg', dst2)

17_canny.jpg

* Sharpening *

Original image + (original image-image after smoothing) x constant

python


# Apply Gaussian Blur
blur = cv2.GaussianBlur(res,(0, 0), sigmaX = 3, sigmaY = 3)
cv2.imwrite('shap_blur.jpg', blur)

diff = res - blur
cv2.imwrite('18_shape_diff.jpg', diff)

shapen = res + diff*3
cv2.imwrite('19_shape_shapen.jpg', shapen)

shap_blur.jpg 18_shape_diff.jpg 19_shape_shapen.jpg

* Local averaging *

It can be averaged while leaving a relatively edge. As an example, the following 9 patterns are calculated. Select the one with the least dispersion.

aria_ave.jpg.png

python


def vari(data):
    count = len(data)
    av = sum(data)/count
    void = []
    total = 0
    for d in data:
        diff = (d - av)**2
        total += diff
    pvariance = total/count

    return pvariance


def selective_ave(image):
    av_image = np.zeros(image.shape)
    
    for i in range(2, image.shape[0] - 2):
        for j in range(2, image.shape[1] - 2):

            center_neighbour = [image[i+1, j-1],image[i+1, j],image[i+1, j+1],image[i, j-1],
                image[i, j+1],image[i-1, j-1],image[i-1, j],image[i-1, j+1]]

            l_u_neighbour = [image[i-2, j-2],image[i-2, j-1],image[i-1, j-2],image[i-1, j-1],
                image[i-1, j],image[i, j-1],image[i, j]]
            
            up_neighbour = [image[i-2, j-1],image[i-2, j],image[i-2, j+1],
                image[i-1, j-1],image[i-1, j],image[i-1, j+1],image[i, j]]

            r_u_neighbour = [image[i+2, j+2],image[i+2, j+1],image[i+1, j+2],image[i+1, j+1],
                image[i+1, j],image[i, j+1],image[i, j]]

            right_neighbour = [image[i-1, j+2],image[i, j+2],image[i+1, j+2],
                image[i-1, j+1],image[i, j+1],image[i+1, j+1],image[i, j]]
            
            left_neighbour = [image[i-1, j-2],image[i, j-2],image[i+1, j-2],
                image[i-1, j-1],image[i, j-1],image[i+1, j-1],image[i, j]]

            l_d_neighbour = [image[i+2, j-2],image[i+2, j-1],image[i+1, j-2],
                image[i+1, j-1],image[i+1, j],image[i, j-1],image[i, j]]

            down_neighbour = [image[i+2, j-1],image[i+2, j],image[i+2, j+1],
                image[i+1, j-1],image[i+1, j],image[i+1, j+1],image[i, j]]

            r_d_neighbour = [image[i+2, j+2],image[i+2, j+1],image[i+1, j+2],
                image[i+1, j+1],image[i+1, j],image[i, j+1],image[i, j]]
            
    
            nb_list = [center_neighbour, l_u_neighbour, up_neighbour,
                r_u_neighbour, right_neighbour, left_neighbour,
                l_d_neighbour, down_neighbour, r_d_neighbour]

            number = 0
            vari_max = 0

            for nb in nb_list:
                pre_vari = vari(nb)
                
                if vari_max < pre_vari:
                    vari_max = pre_vari
                    max_index = number
                
                number += 1
            
            result = statistics.mean(nb_list[max_index])

            av_image[i, j] = result
 
    return av_image

select_av =  selective_ave(res)
cv2.imwrite('20_select_av.jpg', select_av)

20_select_av.jpg

* Bilateral filter (filter using two Gaussian distributions) *

Implementation of smoothing with a Gaussian filter. Weights are distributed according to the distance to the point of interest (the farther the distance, the lighter) Up to this point, it is the same as a normal Gaussian filter. Here, further, the weighting of the Gaussian distribution is adopted for the difference in the pixel value between the pixels. That is, (the larger the difference (the farther the pixel value), the lighter the weight to be smoothed, so the edge portion (the larger the difference in pixel value) is less likely to be smoothed. As a result, the edges are preserved.

python


bi = cv2.bilateralFilter(res,15, 20, 20)
cv2.imwrite('21_bilateral_fil.jpg', bi)

21_bilateral_fil.jpg

* Non-local mean filter *

Whereas the bilateral filter weights the difference in pixel values, Here, the difference in similarity with other small areas (surrounding image areas) is weighted. That is, if there are many similar images in the periphery, the smoothing becomes heavy, and if they are not similar, the smoothing becomes light.

python


dst = cv2.fastNlMeansDenoising(res, None, 10, 7, 21)
cv2.imwrite('22_non_l_m_fil.jpg', dst)

22_non_l_m_fil.jpg

* Median filter *

At the time of filtering, the median value of the area of interest is output. The median is (M × N + 1) / 2, counting from the pixel with the smallest pixel value for an image of M × N pixels. The value of the second pixel. This filter is effective for spike-like noise such as sprinkled with sesame seeds.

python


median = cv2.medianBlur(noicy_image, 3)
cv2.imwrite('24_median_fil.jpg', median)

Before 23_spike_noise.jpg

After 24_median_fil.jpg

Referenced site http://optie.hatenablog.com/entry/2018/03/21/185647 https://theailearner.com/2019/05/25/laplacian-of-gaussian-log/ https://algorithm.joho.info/programming/python/opencv-canny-edge-detecte-py/

that's all

Recommended Posts

Digital image processing (spatial filtering)
Read digital image processing
[Image processing] Posterization
Image processing 100 knocks ①
Image processing with MyHDL
Flat Field image processing
First Python image processing
Image processing with Python
Image Processing with PIL
Image processing with Python (Part 2)
opencv-python Introduction to image processing
Image processing with PIL (Pillow)
100 image processing knocks !! (011 --020) Early game
100 image processing knocks !! (001 --010) Carefully and carefully
Image processing with Python (Part 1)
Image processing with Python (Part 3)
Image processing by python (Pillow)
Image Processing Collection in Python
Image expansion and contraction processing
[Python] Image processing with scikit-image
Real-time image processing basics with opencv
[Python] Using OpenCV with Python (Image Filtering)
Personal notes for python image processing
Image processing with Python 100 knocks # 3 Binarization
Image processing 100 knocks Q9, Q10 (filter) speedup
Environmentally friendly scraping using image processing
Image processing with Python 100 knocks # 2 Grayscale
Image processing | predicting species from images