Python實現高斯平滑與Sobel邊緣檢測

平滑濾波與邊緣檢測是圖像處理中非常基礎與重要的部分。平滑濾波器主要有均值濾波，中值濾波，高斯濾波與雙邊濾波等，邊緣檢測主要有Sobel算子，Laplace算子，Canny算子等。本文主要就高斯濾波與Sobel算子進行原理上的介紹，並用Python進行實現。

第一部分，高斯濾波

原理

高斯濾波是一種線性濾波器，能夠較好地平滑與抑制圖像噪聲，與均值中值濾波一樣，高斯濾波也是對圖像像素進行平均的一個過程，但不同的一點是，高斯濾波對圖像像素進行的是加權平均，每一個像素點的值，都有其鄰域內的其他像素點灰度值加權平均得到。在圖像處理領域的高斯濾波其本質上是一種高通濾波器，過濾掉圖像中的細節部分（高頻部分），保留圖像中的平滑部分（低頻部分）。

正如上文說的，高斯濾波處理在圖像上反映出來的是一個加權平均的過程，距離中心點越近的部分，其權重越大，反之越小，而權值的分布，則符合正態分布，對於一個二維圖像來說，則符合二維高斯分布。

高斯函數

一維高斯函數：

\[G(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}} \]

二維高斯函數：

\[G(x,y) = \frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}} \]

二維高斯函數的圖像為：

圖像卷積

原理

為了下一步討論高斯濾波的實現過程，我們首先要介紹一下圖像處理中的卷積是如果工作和實現的。與信號處理中的卷積實現方法不同，在圖像處理中的卷積主要通過模板矩陣與圖像矩陣相乘，得到新的結果矩陣。下圖詳細演示了圖像卷積的過程：

在上圖中，黃色標識的像素標識我們要卷積的目標像素，綠色為經過卷積核得出的結果像素值，其運算過程為：

\[0\times (-1) + 0\times 0 + 0\times 1+0\times (-2)+1\times 0+ 1\times 2 + 1\times(-1) + 1\times 0 + 1\times 1 = 3 \]

邊緣補償

在實際操作中，我們需要對圖像的邊緣進行補償，以便卷積核可以對圖像邊緣部分的像素進行計算。我們常用的邊緣補償方法有以下幾種：

1. 補0：即補償的像素值為0
2. 復制：即復制邊緣一圈的像素值

在本代碼中，我們用補0的方式來進行處理。

代碼

根據上文對圖像處理卷積的理解，我們可以認為高斯卷積也是一個符合高斯分布的卷積核，在本實例中，我們使用的卷積核大小為\(3\times 3\)。

邊緣補償

如上文所述，在卷積處理開始之前，要對邊緣進行補償，我們這里使用全0補償法，代碼如下：

def padding(image,ksize):# 輸入目標圖像，卷積核大小

    h = image.shape[0] # 獲取圖像尺寸
    w = image.shape[1]
    c = image.shape[2]

    pad = ksize // 2 # 需要補償的邊緣大小

    out_p = np.zeros((h+2*pad,w+2*pad,c)) # 創造一個補償后大小的全0矩陣

    out_copy = image.copy()

    out_p[pad:pad+h,pad:pad+w,0:c] = out_copy.astype(np.uint8) # 將原圖像復制入目標圖像中

    return out_p

高斯矩陣

下來是構建高斯卷積核，根據上圖中的二維高斯函數公式，我們的代碼為：

    # 高斯卷積核

    kernel = np.zeros((ksize,ksize)) # 創建一個卷積核大小的全0二維矩陣
    for x in range(-pad,-pad+ksize):
        for y in range(-pad,-pad+ksize):
            kernel[y+pad,x+pad] = np.exp(-(x**2+y**2)/(2*(sigma**2))) # 給卷積核內賦值
    kernel /= (sigma*np.sqrt(2*np.pi)) # 計算平均值
    kernel /=  kernel.sum() # 加總

要注意的是，高斯方程中有一個很重要的參數就是\(sigma\)值，在本示例中，如果未指定確定值，則\(sigma\)由下面的公式自動確定：

\[sigma = 0.3 \times ((ksize-1)\times 0.5-1) + 0.8 \]

其中ksize為高斯卷積核大小（本例中為\(3\)，即高斯卷積核為\(3\times 3\)）。

完整代碼

最后，再根據圖像中卷積的定義，對目標進行高斯卷積。

整個過程的完整代碼如下：

import numpy as np 
from PIL import Image
from cv2 import cv2
import matplotlib.pyplot as plt
import math




def padding(image,ksize):

    h = image.shape[0]
    w = image.shape[1]
    c = image.shape[2]

    pad = ksize // 2 

    out_p = np.zeros((h+2*pad,w+2*pad,c))

    out_copy = image.copy()

    out_p[pad:pad+h,pad:pad+w,0:c] = out_copy.astype(np.uint8)

    return out_p





def gaussian(image,ksize,sigma):

    """
    1. padding
    2. 定義高斯濾波公式與卷積核
    3. 卷積過程

    高斯卷積卷積核是按照二維高斯分布規律產生的，公式為：

    G(x,y) = (1/(2*pi*(sigma)^2))*e^(-(x^2+y^2)/2*sigma^2)

    唯一的未知量是sigma，在未指定sigma的前提下，可以通過下列參考公式讓程序自動選擇合適的
    sigma值：

    sigma =  0.3 *((ksize-1)*0.5-1) + 0.8

    @ 如果mode為default，則返回abs值，否則返回unit8值
    """
        

    pad = ksize//2

    out_p = padding(image,ksize) # padding之后的圖像
    # print(out_p)

    h = image.shape[0]
    w = image.shape[1]
    c = image.shape[2]

    # 高斯卷積核

    kernel = np.zeros((ksize,ksize))
    for x in range(-pad,-pad+ksize):
        for y in range(-pad,-pad+ksize):
            kernel[y+pad,x+pad] = np.exp(-(x**2+y**2)/(2*(sigma**2)))
    kernel /= (sigma*np.sqrt(2*np.pi))
    kernel /=  kernel.sum()

    # print(kernel)

    tmp = out_p.copy()

    # print(tmp)

    for y in range(h):
        for x in range(w):
            for z in range(c):

                out_p[pad+y,pad+x,z] = np.sum(kernel*tmp[y:y+ksize,x:x+ksize,z])


    out = out_p[pad:pad+h,pad:pad+w].astype(np.uint8)
    # print(out)

    return out
	
if __name__ == "__main__":

    path = "../bilder/lena.png"

    img = cv2.imread(path)
    gaussian_img = gaussian(img,3,0.8)
	cv2.imshow('Original Image',img)
	cv2.imshow('Gaussian Image',gaussian_img)
	cv2.waitkey()

結果

我們用圖像處理中著名的"lena"圖來進行結果展示。

通過結果我們可以看出，高斯圖像相比原圖更加平滑，比如頭發上的細節，經過高斯處理后，平滑了許多，如果想要更多的效果，可以通過增大sigma值來實現。

第二部分 Sobel算子

原理

Sobel算子是非常常用的邊緣檢測算子，那么，既然是邊緣檢測，首先就需要了解什么是邊緣。

何為邊緣？

我們人眼可以非常容易的分辨哪里是邊緣哪里不是，但對於計算機來說，則沒有那么容易。因為對於程序來說，一個灰度圖像的本質是一個二維矩陣，矩陣中的每一個值對應的是圖像每一個像素點的灰度值。那么，對於邊緣來說，就是灰度變化最激烈的地方，換句話說，就是灰度梯度值最大的點，因為梯度就代表了變化程度。

Sobel算子也是一種加權算子，主要思想在於，各個像素對中心像素點的影響不是等價的，會隨着距離的不同而不同。Sobel分別在x和y兩個方向上對像素進行求導。

\[G(x) = \begin{bmatrix} -1&0&1\\ -2&0&2\\ -1&0&1 \end{bmatrix} * I \]

\[G(y) = \begin{bmatrix} -1&-2&-1\\ 0&0&0\\ 1&2&1 \end{bmatrix} * I \]

代碼

在本例實現的代碼中，我們定義了一個參數operator_type來確定是進行x還是y方向上的求導。如果operator_type = 'horizontal'則是x方向，operator_type = 'Vertical'則是y方向。不同求導方向對於不同方向上邊緣的敏感性不同，具體的區別會在后面的結果清楚的呈現。

def sobel(image,operator_type,mode):


    if operator_type == "horizontal": # 定義求導方向

        sobel = np.array([[-1,-2,-1],[0,0,0],[1,2,1]])
    elif operator_type == "vertical":
        sobel = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
    else:
        print("type error") # 如果未定義，則輸出錯誤提醒

    h = image.shape[0]
    w = image.shape[1]
    c = image.shape[2]
    out_p = padding(image,ksize=3)# padding，邊緣補償

    tmp = out_p.copy()

    for y in range(h):
        for x in range(w):
            for z in range(c):

                if mode == 'default': # 為后面的canny算子進行准備，默認對結果像素取絕對值。

                    out_p[1+y,1+x,z] = abs(np.sum(sobel*tmp[y:y+3,x:x+3,z]))
                else:
                    out_p[1+y,1+x,z] = np.sum(sobel*tmp[y:y+3,x:x+3,z])

    out = out_p[1:1+h,1:1+w].astype(np.uint8)

    out = out.clip(0,255) # 裁剪像素，將結果像素值限定在0-255范圍內。

    if mode == "default": # 設定閾值，當前閾值設定為70
        threshold = 70
        for i in range(h):
            for j in range(w):
                for z in range(c):
                    if out[i,j,z] < threshold:
                        out[i,j,z] = 0
    # print(out)

    return out

結果

下面是Sobel在不同方向上的運行結果。

原圖

x方向上的結果

y方向上的結果

我們可以清楚的看出來求導方向對檢測結果的影響，在Sobel算子上，不同方向的運算對方向上邊緣的檢測效果有很大的差別。但很多時候我們需要檢測的是整個圖像的所有邊緣而非特定方向的邊緣，這時候我們就可以用Canny算子，在下一篇博客我會進行Canny算子的相關介紹和實現。

結論

以上就是Python對傳統方法圖像處理中高斯平滑與Sobel算子的原理介紹與代碼實現，代碼其實還有很大的優化空間。