Canny邊緣檢測演算法的實現

阿新 • • 發佈：2019-01-11

影象邊緣資訊主要集中在高頻段，通常說影象銳化或檢測邊緣，實質就是高頻濾波。我們知道微分運算是求訊號的變化率，具有加強高頻分量的作用。在空域運算中來說，對影象的銳化就是計算微分。由於數字影象的離散訊號，微分運算就變成計算差分或梯度。影象處理中有多種邊緣檢測（梯度）運算元，常用的包括普通一階差分，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;
}

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

非最大值抑制結果

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;
    }

下面是對Lena圖計算Canny邊緣檢測的梯度模圖和二值化圖，高斯半徑2，高閥值100，低閥值50。