高通濾波法、微分運算元法、神經網路方法實現邊緣檢測

邊緣檢測（Edge detection）是影像處理和計算機視覺中的基本問題，邊緣檢測的目的是標識數字影像中亮度變化明顯的點。本文使用多種不同的方法，實現對 Lena 肖像的邊緣檢測，研究分析各演算法的效果和優缺點。所涉及的方法如下：

• 高通濾波法

  • 理想高通濾波器
  • Butterworth 高通濾波器
  • 指數高通濾波器

• 微分運算元法

  • Roberts 運算元
  • Sobel 運算元
  • Laplacian 運算元
  • Canny 運算元

• 神經網路方法

  • HED 演算法

高通濾波法

影像中的邊緣或線條等細節部分與影像頻譜的高頻分量相對應，因此採用高通濾波讓高頻分量順利通過，使影像的邊緣或線條細節變得清楚，實現邊緣提取和影像銳化。

常見的高通濾波器包括：理想高通濾波器、Butterworth 高通濾波器、指數高通濾波器等。

理想高通濾波器

理想高通濾波器的傳遞函式 \(H(u, v)\) 滿足下式：

\[H(u,v) = \begin{cases} 1, &D(u,v)>D_0 \cr 0, &D(u,v)>D_0 \end{cases} \]

理想高通濾波器只是一種理想狀況下的濾波器，不能用實際的電子器件實現。

Butterworth 高通濾波器

Butterworth 高通濾波器的傳遞函式 \(H(u,v)\) 如下：

\[H(u,v)=\frac{1}{1+[D_0/D(u,v)]^{2n}} \]

式中，\(n\) 為階數，\(D_0\) 為截止頻率。

Butterworth 高通濾波器在高低頻率間的過渡比較平滑，所以由其得到的輸出影像的振鈴現象不明顯。

指數高通濾波器

指數高通濾波器的傳遞函式 \(H(u,v)\) 如下：

\[H(u,v)=\exp\{-[\frac{D_0}{D(u,v)}]^n\} \]

式中，變數 \(n\) 控制從原點算起的傳遞函式 \(H(u,v)\) 的增長率。

指數高通濾波器的另一種常用的傳遞函式如下式所示：

\[H(u,v)=\exp\{[\ln(\frac{1}{\sqrt 2})][\frac{D_0}{D(u,v)}]^n\} \]

程式碼實現

為了在頻率域中實現高通濾波，先通過傅立葉變換得到影像的頻譜，根據不同濾波器的不同傳遞函式，對頻率進行相應的過濾，最後再對其進行傅立葉反變換，得到濾波後的影像。

傅立葉變換

img = plt.imread('images/lena.bmp')
fft_shift = np.fft.fftshift(np.fft.fft2(img))	# 變換後將零頻分量移到頻譜中心
fft_img = np.log(np.abs(fft_shift))				# 視覺化

實現三種濾波器

def distance(shape):		# 計算每個畫素到中心原點的距離
    n, m = shape
    u = np.arange(n)
    v = np.arange(m)
    u, v = np.meshgrid(u, v)
    return np.sqrt((u - n//2)**2 + (v - m//2)**2)
    
def ideal_filter(shape, d0):
    d = distance(shape)
    mask = d > d0
    return mask.astype(int)

def butterworth_filter(shape, d0, order=1):
    d = distance(shape)
    mask = 1 / (1 + (d0 / d)**(2 * order))
    return mask

def exponential_filter(shape, d0, order=1):
    d = distance(shape)
    mask = np.exp(-(d0 / d)**order)
    return mask

濾波後進行傅立葉反變換

ifft_shift = np.fft.ifftshift(fft_shift * mask)
ifft_img = np.abs(np.fft.ifft2(ifft_shift))

執行結果及分析

下圖中，從上到下依次展示了使用理想高通濾波器、Butterworth 高通濾波器和指數高通濾波器對 Lena 進行邊緣檢測，在不同的截止頻率 \(D_0\) 下（10、20、40、80）所得到的頻譜圖以及濾波後的影像輸出。

通過對比可以發現，使用理想高通濾波器得到的結果有明顯的振鈴現象，而 Butterworth 高通濾波器和指數高通濾波器的結果相近，均具有較好的效果。從頻譜圖中也可以看出，理想高通濾波器對頻率的截斷非常陡峭，在臨界點發生了突變，而後兩者的濾波比較平滑，是一種在高低頻率間逐漸過渡的過程。

微分運算元法

針對由於平均或積分運算而引起的影像模糊，可用微分運算來實現影像的銳化。微分運算是求訊號的變化率，有加強高頻分量的作用，從而使影像輪廓清晰。

常見的邊緣檢測運算元包括：Roberts 運算元、Sobel 運算元、Laplacian 運算元、Canny 運算元等。各種運算元的存在就是對這種導數分割原理進行的例項化計算，是為了在計算過程中直接使用的一種計算單位。實際使用時，通常用各種運算元對應的模板對原圖進行卷積運算，從而提取出影像的邊緣資訊。

上述各運算元的具體定義就不展開贅述了，具體可以參考相關書籍或文章。這裡通過一張表對這些運算元進行簡要的介紹和比較。

運算元	介紹及優缺點比較
Roberts	一種最簡單的運算元，採用對角線方向相鄰兩畫素之差近似梯度幅值檢測邊緣。檢測垂直邊緣的效果好於斜向邊緣，定位精度高，但是對噪聲敏感，對具有陡峭邊緣且含噪聲少的影像效果較好。
Sobel	根據畫素點上下左右四鄰域灰度加權差檢測邊緣，類似區域性平均運算，因此對噪聲具有平滑作用，對灰度漸變和噪聲較多的影像處理效果比較好，對邊緣定位比較準確。
Laplacian	屬於二階微分運算元，在只考慮邊緣點的位置而不考慮周圍的灰度差時適合用該運算元進行檢測。對噪聲非常敏感，只適用於無噪聲影像。存在噪聲的情況下，使用該運算元檢測邊緣之前需要先進行低通濾波，因此通常把 Laplacian 運算元和平滑運算元結合起來生成一個新的模板。
Canny	該運算元功能比前面幾種都要好，不容易受噪聲的干擾，能夠檢測到真正的弱邊緣，但是實現起來較為麻煩，是一個具有濾波、增強、檢測的多階段的優化運算元。在進行處理前，Canny 運算元先利用高斯平滑濾波器來平滑影像以除去噪聲。Canny 分割演算法採用一階偏導的有限差分來計算梯度幅值和方向，在處理過程中，該運算元還將經過一個非極大值抑制的過程，最後採用兩個閾值來連線邊緣。

程式碼實現

Roberts、Sobel 以及 Laplacian 運算元方法均為手工實現：先根據不同運算元對應的模板，分別定義對單個畫素塊的卷積運算函式，然後在影像上滑動模板，呼叫該函式計算每一個卷積塊，最終得到經過各運算元微分後的輸出影像。對於 Canny 運算元，由於其實現過程比較複雜，這裡選擇直接呼叫 OpenCV 中的 Canny() 函式。

實現前三種運算元

def roberts_operator(block):
    kernel1 = np.array([[1,0], [0,-1]])
    kernel2 = np.array([[0,-1], [1,0]])
    return np.abs(np.sum(block[1:,1:] * kernel1)) + np.abs(np.sum(block[1:,1:] * kernel2))

def sobel_operator(block, orientation):		# 水平和垂直兩個方向
    if orientation == 'horizontal':
        kernel = np.array([[-1,-2,-1], [0,0,0], [1,2,1]])
    elif orientation == 'vertical':
        kernel = np.array([[-1,0,1], [-2,0,2], [-1,0,1]])
    else:
        raise('Orientation Error')
    return np.abs(np.sum(block * kernel))

def laplacian_operator(block):
    kernel = np.array([[0,-1,0], [-1,4,-1], [0,-1,0]])
    return np.abs(np.sum(block * kernel))

滑動計算卷積塊

def operator_process(img, operator_type, orientation=None):
    n, m = img.shape
    res = np.zeros((n, m))
    for i in range(1, n-1):
        for j in range(1, m-1):
            if operator_type == 'roberts':
                res[i][j] = roberts_operator(img[i-1:i+2, j-1:j+2])
            elif operator_type == 'sobel':
                res[i][j] = sobel_operator(img[i-1:i+2, j-1:j+2], orientation)
            elif operator_type == 'laplacian':
                res[i][j] = laplacian_operator(img[i-1:i+2, j-1:j+2])
            else:
                raise('Operator Type Error')
    return res

呼叫 OpenCV 中的 Canny 運算元

canny_res = cv2.Canny(img, threshold1, threshold2)		# 設定高低閾值引數

執行結果及分析

使用 Roberts 運算元、Laplacian 運算元以及水平和垂直兩個方向的 Sobel 運算元進行邊緣檢測所得到的結果如下圖所示：

可以看到，Roberts 運算元簡單但有效，已經能實現比較好的邊緣檢測效果；而 Laplacian 運算元的效果相對較差，邊緣不是很清晰，還出現了很多噪點；Sobel 運算元整體表現也較好，可以明顯地看出水平和垂直兩個方向上結果的差別，前者對水平邊緣響應最大（如眼眶、嘴脣、下巴處），後者對垂直邊緣響應最大（如鼻樑、兩側臉頰處）。

下圖顯示了使用 OpenCV 自帶的 Canny 運算元在不同閾值下的輸出結果，其中低閾值 threshold1 依次取 30、50、80、120、150，而高閾值 threshold2 根據 Canny 演算法的推薦，均取為 threshold1 的 3 倍。

低於 threshold1 的畫素被認為不是邊緣，高於 threshold2 的畫素被認為是邊緣，介於二者之間的則會根據相鄰的畫素點進一步確定。從圖中也可以看出，閾值設定得越高，對邊緣的過濾就越嚴格，輸出結果中的邊緣線條也越發稀疏。

噪聲測試

為了進一步研究在有噪聲的情況下各運算元的邊緣檢測效果，首先使用 skimage 庫中的函式，對 Lena 加入高斯噪聲和椒鹽噪聲。

img_noise1 = skimage.util.random_noise(img, mode='gaussian')
img_noise2 = skimage.util.random_noise(img, mode='s&p')

然後分別在這兩張有噪聲的影像上應用 Roberts 運算元、Sobel 運算元（水平方向）、Laplacian 運算元和 Canny 運算元，得到的結果如下：

其中，第一行是加入高斯噪聲後各運算元的輸出，第一行是加入椒鹽噪聲後各運算元的輸出。

可以發現，Roberts 運算元和 Sobel 運算元都有一定的抗噪聲能力，從圖中依然可以看出部分邊緣資訊，而在 Laplacian 運算元則對噪聲非常敏感，其輸出結果完全看不出任何邊緣。對於 Canny 運算元，需要將閾值設定得很高，才能得到比較好的效果，否則大量噪點也會被認為是邊緣。經過多次嘗試，最終將 threshold1 設為 180，將 threshold2 設為 3 * threshold1，得到了上圖中最後一列的結果。

實驗過程中還發現，在相同的閾值條件下，Canny 運算元對高斯噪聲的抵抗能力比對椒鹽噪聲的抵抗能力強。下圖展示了當 threshold1 取值為 100、120、140、160、180 時，Canny 運算元對加入了高斯噪聲和椒鹽噪聲的影像的邊緣檢測效果：

HED 演算法

2015年，Saining Xie 等人提出了一種基於卷積神經網路的邊緣檢測演算法——Holistically-Nested Edge Detection（HED）演算法。模型使用 VGG-16 作為骨幹網路進行多尺度多層級的特徵學習。其中 Holistically 的意思是”整體地“，表示該演算法試圖訓練一個 image-to-image 的網路；Nested 則強調在生成的輸出過程中，通過不斷的整合和學習，得到更精確的邊緣預測圖。

HED 演算法具有很多優點。單就預測過程來說，對一張圖片進行邊緣檢測的速度是很快的。使用時，根據圖片的實際內容，可以通過調整合適的超引數從而得到最優尺度的邊緣資訊。另外，由於 HED 本身就是基於神經網路的，因而可以很方便地嵌入其他網路模型中，直接參與各種學習任務的訓練過程。

在效果上，HED 演算法也具有優越性。論文中，作者把 HED 和傳統 Canny 演算法進行對比。如下圖所示，可以看到 HED 的效果明顯優於 Canny 演算法。關於 HED 演算法的更多內容請參考論文原文。

程式碼實現

這裡採用的是 HED 演算法的另一個 PyTorch 實現版本，直接使用了作者提供的預訓練模型進行預測。

執行結果及分析

模型在彩色 Lena 影像上的執行結果如下圖所示：

在新增了高斯噪聲和椒鹽噪聲後的影像上，執行結果分別如下所示：

新增高斯噪聲後的 HED 輸出影像

新增椒鹽噪聲後的 HED 輸出影像

可見，雖然 HED 演算法的邊緣線條比較粗，但整體表現還是相當優秀的，尤其是在存在噪聲的情況下，該演算法的效果比前述幾種基於微分運算元的方法都要好。

總結

本文通過使用高通濾波法、微分運算元法、神經網路方法三大類，共計 8 種不同的方法對 Lena 影像進行了邊緣檢測，將各種方法得到的結果進行橫向比較，並對它們的優缺點和適用場景進行了一定的討論。對於具有引數的演算法，進一步根據引數取值的不同進行了縱向比較，觀察引數對於輸出結果的影響。此外，還通過對原影像新增高斯噪聲和椒鹽噪聲進行噪聲測試，研究各演算法對兩種噪聲的敏感性。

實驗結果表明，不同演算法由於原理或核心函式的不同，均具有各自的優缺點和適用場景。使用時應根據影像內容和實際需求進行選擇取捨，並通過調整相關引數從而到達最佳的效果。

完整原始碼請見 GitHub 倉庫

參考資料

1. 王一丁，李琛，王蘊紅.《數字影像處理》.西安電子科技大學出版社
2. Xie, Saining, and Zhuowen Tu. "Holistically-nested edge detection." Proceedings of the IEEE international conference on computer vision. 2015.
3. PyTorch-HED：https://github.com/sniklaus/pytorch-hed
4. 網際網路上的一些部落格、文章、資料