Canny邊緣檢測算法的實現

本文轉載自查看原文 2017-02-14 15:12 70352

圖像邊緣信息主要集中在高頻段，通常說圖像銳化或檢測邊緣，實質就是高頻濾波。我們知道微分運算是求信號的變化率，具有加強高頻分量的作用。在空域運算中來說，對圖像的銳化就是計算微分。由於數字圖像的離散信號，微分運算就變成計算差分或梯度。圖像處理中有多種邊緣檢測（梯度）算子，常用的包括普通一階差分，Robert算子（交叉差分），Sobel算子等等，是基於尋找梯度強度。拉普拉斯算子（二階差分）是基於過零點檢測。通過計算梯度，設置閥值，得到邊緣圖像。

Canny邊緣檢測算子是一種多級檢測算法。1986年由John F. Canny提出，同時提出了邊緣檢測的三大准則：

低錯誤率的邊緣檢測：檢測算法應該精確地找到圖像中的盡可能多的邊緣，盡可能的減少漏檢和誤檢。
最優定位：檢測的邊緣點應該精確地定位於邊緣的中心。
圖像中的任意邊緣應該只被標記一次，同時圖像噪聲不應產生偽邊緣。

Canny算法出現以后一直是作為一種標准的邊緣檢測算法，此后也出現了各種基於Canny算法的改進算法。時至今日，Canny算法及其各種變種依舊是一種優秀的邊緣檢測算法。而且除非前提條件很適合，你很難找到一種邊緣檢測算子能顯著地比Canny算子做的更好。

關於各種差分算子，還有Canny算子的簡單介紹，這里就不羅嗦了，網上都可以找得到。直接進入Canny算法的實現。Canny算法分為幾步。

1. 高斯模糊。

這一步很簡單，類似於LoG算子（Laplacian of Gaussian）作高斯模糊一樣，主要作用就是去除噪聲。因為噪聲也集中於高頻信號，很容易被識別為偽邊緣。應用高斯模糊去除噪聲，降低偽邊緣的識別。但是由於圖像邊緣信息也是高頻信號，高斯模糊的半徑選擇很重要，過大的半徑很容易讓一些弱邊緣檢測不到。

Lena原圖 Lena高斯模糊，半徑2

2. 計算梯度幅值和方向。

圖像的邊緣可以指向不同方向，因此經典Canny算法用了四個梯度算子來分別計算水平，垂直和對角線方向的梯度。但是通常都不用四個梯度算子來分別計算四個方向。常用的邊緣差分算子（如Rober，Prewitt，Sobel）計算水平和垂直方向的差分Gx和Gy。這樣就可以如下計算梯度模和方向：

梯度角度θ范圍從弧度-π到π，然后把它近似到四個方向，分別代表水平，垂直和兩個對角線方向（0°,45°,90°,135°）。可以以±iπ/8（i=1,3,5,7）分割，落在每個區域的梯度角給一個特定值，代表四個方向之一。

這里我選擇Sobel算子計算梯度。Sobel算法很簡單，到處都可以找到，就不列出代碼來了。相對於其他邊緣算子，Sobel算子得出來的邊緣粗大明亮。

下圖是對上面半徑2的高斯模糊圖像L通道（HSL）應用Sobel算子的梯度模圖，沒有施加任何閥值。

Sobel算子，無閥值

3. 非最大值抑制。

非最大值抑制是一種邊緣細化方法。通常得出來的梯度邊緣不止一個像素寬，而是多個像素寬。就像我們所說Sobel算子得出來的邊緣粗大而明亮，從上面Lena圖的Sobel結果可以看得出來。因此這樣的梯度圖還是很“模糊”。而准則3要求，邊緣只有一個精確的點寬度。非最大值抑制能幫助保留局部最大梯度而抑制所有其他梯度值。這意味着只保留了梯度變化中最銳利的位置。算法如下：

比較當前點的梯度強度和正負梯度方向點的梯度強度。
如果當前點的梯度強度和同方向的其他點的梯度強度相比較是最大，保留其值。否則抑制，即設為0。比如當前點的方向指向正上方90°方向，那它需要和垂直方向，它的正上方和正下方的像素比較。

注意，方向的正負是不起作用的，比如東南方向和西北方向是一樣的，都認為是對角線的一個方向。前面我們把梯度方向近似到水平，垂直和兩個對角線四個方向，所以每個像素根據自身方向在這四個方向之一進行比較，決定是否保留。這一部分的代碼也很簡單，列出如下。pModule，pDirection分別記錄了上一步梯度模值和梯度方向。

pmoddrow = pModule + Width + 1; pdirdrow = pDirection + Width + 1; pstrongdrow = pStrong + Width + 1; for (i = 1; i < Hend - 1; i++) { 　　pstrongd = pstrongdrow; 　　pmodd = pmoddrow; 　　pdird = pdirdrow; 　　for (j = 1; j < Wend - 1; j++) 　　　　{ 　 switch (*pdird) { case 0:        // x direction
            case 4: if (*pmodd > *(pmodd - 1) && *pmodd > *(pmodd + 1)) *pstrongd = 255; break; case 1:        // northeast-southwest direction. Notice the data order on y direction of bmp data
            case 5: if (*pmodd > *(pmodd + Width + 1) && *pmodd > *(pmodd - Width - 1)) *pstrongd = 255; break; case 2:        // y direction
            case 6: if (*pmodd > *(pmodd - Width) && *pmodd > *(pmodd + Width)) *pstrongd = 255; break; case 3:        // northwest-southeast direction. Notice the data order on y direction of bmp data
            case 7: if (*pmodd > *(pmodd + Width - 1) && *pmodd > *(pmodd - Width + 1)) *pstrongd = 255; break; default: ASSERT(0); break; } pstrongd++; pmodd++; pdird++; } pstrongdrow += Width; pmoddrow += Width; pdirdrow += Width; }

View Code

下圖是非最大值抑制的結果。可見邊緣寬度已經大大減小。但是這個圖像中因為沒有應用任何閥值，還含有大量小梯度模值的點，也就是圖中很暗的地方。下面，閥值要上場了。

非最大值抑制結果

4. 雙閥值。

一般的邊緣檢測算法用一個閥值來濾除噪聲或顏色變化引起的小的梯度值，而保留大的梯度值。Canny算法應用雙閥值，即一個高閥值和一個低閥值來區分邊緣像素。如果邊緣像素點梯度值大於高閥值，則被認為是強邊緣點。如果邊緣梯度值小於高閥值，大於低閥值，則標記為弱邊緣點。小於低閥值的點則被抑制掉。這一步算法很簡單。

5. 滯后邊界跟蹤。

至此，強邊緣點可以認為是真的邊緣。弱邊緣點則可能是真的邊緣，也可能是噪聲或顏色變化引起的。為得到精確的結果，后者引起的弱邊緣點應該去掉。通常認為真實邊緣引起的弱邊緣點和強邊緣點是連通的，而又噪聲引起的弱邊緣點則不會。所謂的滯后邊界跟蹤算法檢查一個弱邊緣點的8連通領域像素，只要有強邊緣點存在，那么這個弱邊緣點被認為是真是邊緣保留下來。

這個算法搜索所有連通的弱邊緣，如果一條連通的弱邊緣的任何一個點和強邊緣點連通，則保留這條弱邊緣，否則抑制這條弱邊緣。搜索時可以用廣度優先或者深度優先算法，我在這里實現了應該是最容易的深度優先算法。一次連通一條邊緣的深度優先算法如下：

准備一個棧s，一個隊列q，設聯通指示變量connected為假。從圖像的第一個點開始，進入2。
如果這個點是弱邊界點並且沒有被標記，把它標記，並把它作為第一個元素放入棧s中，同時把它放入記錄連通曲線的隊列q，進入3。如果這個點不是弱邊界或者已經被標記過，到圖像的下一個點，重復2。
從棧s中取出一個元素，查找它的8像素領域。如果一個領域像素是弱邊界並且沒有被標記過，把這個領域像素標記，並加入棧s中，同時加入隊列q。同時查找領域對應的強邊界圖，如果有一個像素是強邊界，表示這條弱邊界曲線和強邊界聯通，設置connected為真。重復3直到棧中沒有元素了。如果connected為假，則依次從隊列q中取出每個元素，清空標記。如果connected為真，保留標記。
清空隊列q，設置connected為假，移動到圖像的下一個點，到2。

// 5. Edge tracking by hysteresis
    stack<CPoint> s; queue<CPoint> q; BOOL connected = FALSE; long row_idx = Width; for (i = 1; i < Height - 1; i++, row_idx += Width) { for (j = 1; j < Width - 1; j++) { pweakd = pWeak + row_idx + j; if (*pweakd == 255) { s.push(CPoint(j, i)); q.push(CPoint(j, i)); *pweakd = 1;        // Label it

                while (!s.empty()) { CPoint p = s.top(); s.pop(); // Search weak edge 8-point neighborhood
                    pweakd = pWeak + p.y*Width + p.x; if (*(pweakd - Width - 1) == 255) { CPoint np = CPoint(p.x - 1, p.y - 1); s.push(np); q.push(np); *(pweakd - Width - 1) = 1;        // Label it
 } if (*(pweakd - Width) == 255) { CPoint np = CPoint(p.x, p.y - 1); s.push(np); q.push(np); *(pweakd - Width) = 1;        // Label it
 } if (*(pweakd - Width + 1) == 255) { CPoint np = CPoint(p.x + 1, p.y - 1); s.push(np); q.push(np); *(pweakd - Width + 1) = 1;        // Label it
 } if (*(pweakd - 1) == 255) { CPoint np = CPoint(p.x - 1, p.y); s.push(np); q.push(np); *(pweakd - 1) = 1;        // Label it
 } if (*(pweakd + 1) == 255) { CPoint np = CPoint(p.x + 1, p.y); s.push(np); q.push(np); *(pweakd + 1) = 1;        // Label it
 } if (*(pweakd + Width - 1) == 255) { CPoint np = CPoint(p.x - 1, p.y + 1); s.push(np); q.push(np); *(pweakd + Width - 1) = 1;        // Label it
 } if (*(pweakd + Width) == 255) { CPoint np = CPoint(p.x, p.y + 1); s.push(np); q.push(np); *(pweakd + Width) = 1;        // Label it
 } if (*(pweakd + Width + 1) == 255) { CPoint np = CPoint(p.x + 1, p.y + 1); s.push(np); q.push(np); *(pweakd + Width + 1) = 1;        // Label it
 } // Search strong edge 8-point neighborhood
                    if (connected == FALSE) { pstrongd = pStrong + p.y*Width + p.x; for (int m = -1; m <= 1; m++) { for (int n = -1; n <= 1; n++) { if (*(pstrongd + m*Width + n) == 255) { connected = TRUE; goto next; } } } } next: continue; } // No more element in the stack
                if (connected == FALSE) { // The weak edge is not connected to any strong edge. Suppress it.
                    while (!q.empty()) { CPoint p = q.front(); q.pop(); pWeak[p.y*Width + p.x] = 0; } } else { // Clean the queue
                    while (!q.empty()) q.pop(); connected = FALSE; } } } } // Add the connected weak edges (labeled) into strong edge image. // All strong edge pixels are labeled 255, otherwise 0.
    for (i = 0; i < len; i++) { if (pWeak[i] == 1) pStrong[i] = 255; }