學習筆記-canny邊緣檢測

時間 2020-02-18

標籤學習筆記 canny 邊緣檢測简体版

原文原文鏈接

Canny邊緣檢測

聲明：閱讀本文須要瞭解線性代數裏面的點乘（圖像卷積的原理），高等數學裏的二元函數的梯度，極大值定義，瞭解機率論裏的二維高斯分佈html

1.canny邊緣檢測原理和簡介算法

2.實現步驟app

3.總結ide

1、 Canny邊緣檢測算法的發展歷史

　　Canny算子是28歲的John Canny在1986年提出的，該文章發表在PAMI頂級期刊(1986. A computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 8, 1986, pp. 679-698)。如今老大爺目前（61歲）在加州伯克利作machine learning,主頁(http://www.cs.berkeley.edu/~jfc/)，大爺就是大爺。函數

　　邊緣檢測是從圖像中提取有用的結構信息的一種技術，若是學過信息論就會知道，一面充滿花紋的牆要比一面白牆的信息量大不少，沒學過也不要緊，直觀上也能理解：充滿花紋的圖像要比單色圖像信息更豐富。爲何要檢測邊緣?由於咱們須要計算機自動的提取圖像的底層（紋理等）或者高層（時間地點人物等）的信息，邊緣能夠說是最直觀、最容易發現的一種信息了。Canny提出了一個對於邊緣檢測算法的評價標準，包括：性能

1) 以低的錯誤率檢測邊緣，也即意味着須要儘量準確的捕獲圖像中儘量多的邊緣。學習

2) 檢測到的邊緣應精肯定位在真實邊緣的中心。spa

3) 圖像中給定的邊緣應只被標記一次，而且在可能的狀況下，圖像的噪聲不該產生假的邊緣。.net

　　簡單來講就是，檢測算法要作到：邊緣要全，位置要準，抵抗噪聲的能力要強。 設計

　　接下來介紹最經典的canny邊緣檢測算法，不少邊緣檢測算法都是在此基礎上進行改進的，學習它有利於一通百通。

2、實現步驟

　　step1：高斯平滑濾波

沒有哪張圖片是沒有噪聲的。————魯迅

　　濾波是爲了去除噪聲，選用高斯濾波也是由於在衆多噪聲濾波器中，高斯表現最好（表現怎麼定義的？最好好到什麼程度？），你也能夠試試其餘濾波器如均值濾波、中值濾波等等。一個大小爲(2k+1)x(2k+1)的高斯濾波器核（核通常都是奇數尺寸的）的生成方程式由下式給出：

‘

　　下面是一個sigma = 1.4，尺寸爲3x3的高斯卷積核的例子，注意矩陣求和值爲1（歸一化）：

　　舉個例子：若圖像中一個3x3的窗口爲A，要濾波的像素點爲e，則通過高斯濾波以後，像素點e的亮度值爲：

　　其中*爲卷積符號，sum表示矩陣中全部元素相加求和，簡單說，就是濾波後的每一個像素值=其原像素中心值及其相鄰像素的加權求和。圖像卷積是圖像處理中很是重要且普遍使用的操做，必定要理解好。

其中高斯卷積核的大小將影響Canny檢測器的性能。尺寸越大，去噪能力越強，所以噪聲越少，但圖片越模糊，canny檢測算法抗噪聲能力越強，但模糊的反作用也會致使定位精度不高，通常狀況下，推薦尺寸5*5，3*3也行。

　　step2：計算梯度強度和方向

邊緣的最重要的特徵是灰度值劇烈變化，若是把灰度值當作二元函數值，那麼灰度值的變化能夠用二元函數的」導數「（或者稱爲梯度）來描述。因爲圖像是離散數據，導數能夠用差分值來表示，差分在實際工程中就是灰度差，說人話就是兩個像素的差值。一個像素點有8鄰域，那麼分上下左右斜對角，所以Canny算法使用四個算子來檢測圖像中的水平、垂直和對角邊緣。算子是以圖像卷積的形式來計算梯度，好比Roberts，Prewitt，Sobel等，這裏選用Sobel算子來計算二維圖像在x軸和y軸的差分值（這些數字的由來？），將下面兩個模板與原圖進行卷積，得出x和y軸的差分值圖，最後計算該點的梯度G和方向θ

　　計算梯度的模和方向屬於高等數學部分的內容，若是不理解應該補習一下數學基本功，圖像處理常常會用到這個概念。

　　這部分我實現了下，首先了解opencv的二維濾波函數：dst=cv.filter2D(src, ddepth, kernel[, dst[, anchor[, delta[, borderType]]]])

　　dst：　　輸出圖片

　　src：　　輸入圖片

　　ddepth：輸出圖片的深度, 詳見 combinations,若是填-1，那麼就跟與輸入圖片的格式相同。

　　kernel：單通道、浮點精度的卷積核。

　　如下是默認參數：

　　anchor：內核的基準點(anchor)，其默認值爲(-1,-1)表示位於kernel的中心位置。基準點即kernel中與進行處理的像素點重合的點。舉個例子就是在上面的step1中，e=H*A獲得的e是放在原像素的3*3的哪個位置，通常來講都是放在中間位置，設置成默認值就好。

　　delta ：在儲存目標圖像前可選的添加到像素的值，默認值爲0。（沒用過）

　　borderType：像素向外逼近的方法，默認值是BORDER_DEFAULT,即對所有邊界進行計算。（沒用過）

　　上代碼

 1 import cv2
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 
 6 img=cv2.imread("images/luxun.png",cv2.IMREAD_GRAYSCALE)  # 讀入圖片
 7 sobel_x = np.array([[-1, 0, 1],[-2,0,+2],[-1, 0, 1]])    # sobel的x方向算子
 8 sobel_y = np.array([[1, 2, 1],[0,0,0],[-1, -2, -1]])     # sobel的x方向算子
 9 sobel_x=cv2.flip(sobel_x,-1)          # cv2.filter2D()計算的是相關，真正的卷積須要翻轉，在進行相關計算。
10 sobel_x=cv2.flip(sobel_y,-1)
11 # cv2.flip()第二個參數：等於0：沿着x軸反轉。大於0：沿着y軸反轉。小於零：沿着x軸，y軸同時反轉
12 
13 # 卷積 opencv是用濾波器函數實現的
14 img_x=cv2.filter2D(img,-1, sobel_x)
15 img_y=cv2.filter2D(img,-1, sobel_y)
16 # 畫圖 plt不支持中文，可是能夠經過如下方法設置修復
17 plt.rcParams['font.sans-serif']=['SimHei']
18 plt.rcParams['axes.unicode_minus'] = False
19 
20 plt.subplot(221), plt.imshow(img_x, 'gray'),plt.title('sobel_x')
21 plt.subplot(222), plt.imshow(img_y, 'gray'),plt.title('sobel_y')
22 plt.subplot(223), plt.imshow(img_y+img_x, 'gray'),plt.title('sobel')
23 plt.subplot(224), plt.imshow(img, 'gray'),plt.title('原圖')
24 plt.show()

View Code

　　運行效果：

　　　　須要注意一點：在圖像處理領域，卷積運算的定義是先將核關於x軸和y軸反轉，然在作相關運算。然而工程實踐中每每跳過反轉，用相關運算代替卷積（好比opencv）。若是你須要嚴格的卷積運算，應該注意原函數的具體實現方式。sobel算子天生關於中心對稱，因此反轉與否並不影響結果（我在代碼裏用cv2.flip（）進行了反轉操做）。

　　　　在以後的實現中，我發現用opencv自帶的濾波函數算出來的梯度是歸一化到（0-255）的，引入其餘的庫也很麻煩，所以本身寫了個簡單的二位卷積函數來實現梯度計算。因此上面的圖適合看效果，並不適合在程序中使用，卷積函數的代碼以下：

 1 def conv2d(src,kernel):  # 輸入必須爲方形卷積核
 2     # 本函數仍然是相關運算，沒有反轉。若是非要嚴格的卷積運算，把下面一行代碼的註釋取消。
 3     #kernel=cv2.flip(kernel,-1)
 4     [rows,cols] = kernel.shape
 5     border=rows//2                                  # 向下取整 得到卷積核邊長
 6     [rows,cols]=src.shape
 7     dst = np.zeros(src.shape)                       # 採用零填充再卷積，卷積結果不會變小。
 8     # print("圖像長：",rows,"寬：",cols,"核邊界",border)
 9     # print(border,rows-border,border,cols-border)
10     temp=[]
11     for i in range(border,rows-border):
12         for j in range(border,cols-border):
13             temp=src[i-border:i+border+1,j-border:j+border+1]  # 從圖像獲取與核匹配的圖像
14             # 切片語法：索引位置包括開頭但不包括結尾 [start: end: step]
15             dst[i][j]=(kernel*temp).sum()  # 計算卷積
16     return dst

View Code

　　　　小技巧：用plt顯示二維矩陣，鼠標移到某個像素就會顯示座標（x，y）和灰度值，浮點數也能夠顯示。這能夠很方便的看某個數據（像素點）是否有問題。

　　　　梯度和幅值的計算效果以下：

　　　　　　能看出來sobel算子計算的邊緣很粗很亮，比較明顯，可是不夠精確，咱們的目標是精確到一個像素寬，至於梯度相位就很難看出什麼特徵，而且梯度相位其實是爲了下一步打基礎的。下面附上代碼：

1 img_x=conv2d(img,sobel_x)  # 使用我本身的寫的卷積計算梯度
2 img_y=conv2d(img,sobel_y)
3 G=np.sqrt(img_x*img_x+img_y*img_y)  # 梯度幅值
4 theta=np.arctan(img_y,(img_x+0.0000000000001))*180/np.pi  # 化爲角度，分母+極小值是爲了不除以0
5 # plt.imshow(theta, 'gray'),plt.title('梯度相位')
6 plt.imshow(G, 'gray'),plt.title('梯度幅值')
7 plt.show()

View Code

　　step3:非極大值抑制

　　　　　　sobel算子檢測出來的邊緣太粗了，咱們須要抑制那些梯度不夠大的像素點，只保留最大的梯度，從而達到瘦邊的目的。這些梯度不夠大的像素點極可能是某一條邊緣的過渡點。按照高數上二位函數的極大值的定義，即對點（x0，y0）的某個鄰域內全部（x，y）都有f（x，y）≤(f（x0，y0），則稱f在（x0，y0）具備一個極大值，極大值爲f（x0，y0）。簡單方案是判斷一個像素點的8鄰域與中心像素誰更大，但這很容易篩選出噪聲，所以咱們須要用梯度和梯度方向來輔助肯定。

　　　　　　以下圖所示，中心像素C的梯度方向是藍色直線，那麼只需比較中心點C與dTmp1和dTmp2的大小便可。因爲這兩個點的像素不知道，假設像素變化是連續的，就能夠用g一、g2和g三、g4進行線性插值估計。設g1的幅值M(g1)，g2的幅值M(g2),則M(dtmp1)=w*M(g2)+(1-w)*M(g1) ，其w=distance(dtmp1,g2)/distance(g1,g2) 。也就是利用g1和g2到dTmp1的距離做爲權重，來估計dTmp1的值。w在程序中能夠表示爲tan（θ）來表示，具體又分爲四種狀況（下面右圖）討論。

　　　　　　以下圖，通過非極大值抑制能夠很明顯的看出去除了不少點，邊緣也變得很細。在程序實現中，要注意opencv的默認座標系是從左到右爲x軸，從上到下是y軸，原點位於左上方，計算g一、g二、g三、g4的位置的時候，必定要當心（坑了我好久）。通過非極大值抑制能夠看出來圖片的邊緣明顯變細，不少看起來黑色的部分其實有值的，只是由於值過小了看不清楚，而這些黑色的部分多是噪聲或者其餘緣由形成的局部極大值，下一步咱們就要用雙閾值來限定出強邊緣和弱邊緣，儘量的減小噪聲的檢出。代碼附上：

 1 # step3:非極大值抑制
 2 anchor=np.where(G!=0)  # 獲取非零梯度的位置
 3 Mask=np.zeros(img.shape)
 4 
 5 for i in range(len(anchor[0])):
 6     x=anchor[0][i]
 7     y=anchor[1][i]
 8     center_point=G[x,y]
 9     current_theta=theta[x,y]
10     dTmp1=0
11     dTmp2=0
12     W=0
13     if current_theta>=0 and current_theta<45:
14         #        g1         第一種狀況
15         # g4  C  g2
16         # g3
17         g1 = G[x + 1, y - 1]
18         g2 = G[x + 1, y]
19         g3 = G[x - 1, y + 1]
20         g4 = G[x - 1, y]
21         W=abs(np.tan(current_theta*np.pi/180))  # tan0-45範圍爲0-1
22         dTmp1= W*g1+(1-W)*g2
23         dTmp2= W*g3+(1-W)*g4
24 
25     elif current_theta>=45 and current_theta<90:
26         #     g2 g1         第二種狀況
27         #     C
28         # g3  g4
29 
30         g1 = G[x + 1, y - 1]
31         g2 = G[x, y - 1]
32         g3 = G[x - 1, y + 1]
33         g4 = G[x, y + 1]
34         W = abs(np.tan((current_theta-90) * np.pi / 180))
35         dTmp1= W*g1+(1-W)*g2
36         dTmp2= W*g3+(1-W)*g4
37 
38     elif current_theta>=-90 and current_theta<-45:
39         # g1  g2            第三種狀況
40         #     C
41         #     g4 g3
42         g1 = G[x - 1, y - 1]
43         g2 = G[x, y - 1]
44         g3 = G[x + 1, y + 1]
45         g4 = G[x, y + 1]
46         W = abs(np.tan((current_theta-90) * np.pi / 180))
47         dTmp1= W*g1+(1-W)*g2
48         dTmp2= W*g3+(1-W)*g4
49 
50     elif current_theta>=-45 and current_theta<0:
51         # g3                第四種狀況
52         # g4  C  g2
53         #        g1
54         g1 = G[x + 1, y + 1]
55         g2 = G[x + 1, y]
56         g3 = G[x - 1, y - 1]
57         g4 = G[x - 1, y]
58         W = abs(np.tan(current_theta * np.pi / 180))
59         dTmp1= W*g1+(1-W)*g2
60         dTmp2= W*g3+(1-W)*g4
61 
62     if dTmp1<center_point and dTmp2<center_point:   # 記錄極大值結果
63             Mask[x,y]=center_point
64 #Mask=(Mask-Mask.min())/(Mask.max()-Mask.min())*256 #歸一化
65 plt.imshow(Mask,'gray'),plt.title('Mask')
66 plt.show()

View Code

　　step4：用雙閾值算法檢測和鏈接邊緣

　　　　　　雙閾值法很是簡單，咱們假設兩類邊緣：通過非極大值抑制以後的邊緣點中，梯度值超過T1的稱爲強邊緣，梯度值小於T1大於T2的稱爲弱邊緣，梯度小於T2的不是邊緣。能夠確定的是，強邊緣必然是邊緣點，所以必須將T1設置的足夠高，以要求像素點的梯度值足夠大（變化足夠劇烈），而弱邊緣多是邊緣，也多是噪聲，如何判斷呢？當弱邊緣的周圍8鄰域有強邊緣點存在時，就將該弱邊緣點變成強邊緣點，以此來實現對強邊緣的補充。實際中人們發現T1:T2=2:1的比例效果比較好，其中T1能夠人爲指定，也能夠設計算法來自適應的指定，好比定義梯度直方圖的前30%的分界線爲T1，我實現的是人爲指定閾值。檢查8鄰域的方法叫邊緣滯後跟蹤，鏈接邊緣的辦法還有區域生長法等等。強邊緣、弱邊緣、綜合效果、和opencv的canny函數對好比下：

3、總結

　　　　實現結果仍是很打擊的，我檢測到的邊緣過於斷續，沒有opencv實現的效果好。查了一下opencv的源碼，這裏猜想兩個可能的緣由：源碼裏梯度的方向被近似到四個角度之一（0，45，90，135），但我用線性插值的的結果是梯度方向更精確，而過於精確-->過於嚴格-->容易受到噪聲干擾，因此在非極大值抑制這以後，我比opencv少了更多的點，最終致使了邊緣不夠連續；第二個緣由多是邊緣鏈接算法效果不夠好，把圖象放大來看，我產生的邊緣傾向於對角線上鏈接，而opencv的邊緣傾向於折線鏈接，所以opencv的邊緣更完整連續，而個人邊緣更細，更容易斷續。

　　　　限於時間，暫時研究到這裏，但願各位多多指正！感謝全部我參考過的博客和文檔！

  1 import cv2
  2 import numpy as np
  3 import matplotlib.pyplot as plt
  4 
  5 # 畫圖 plt不支持中文，可是能夠經過如下方法設置修復
  6 plt.rcParams['font.sans-serif']=['SimHei']
  7 plt.rcParams['axes.unicode_minus'] = False
  8 
  9 def conv2d(src,kernel):  # 輸入必須爲方形卷積核
 10     # 本函數仍然是相關運算，沒有反轉。若是非要嚴格的卷積運算，把下面一行代碼的註釋取消。
 11     #kernel=cv2.flip(kernel,-1)
 12     [rows,cols] = kernel.shape
 13     border=rows//2                                  # 向下取整 得到卷積核邊長
 14     [rows,cols]=src.shape
 15     dst = np.zeros(src.shape)                       # 採用零填充再卷積，卷積結果不會變小。
 16     # print("圖像長：",rows,"寬：",cols,"核邊界",border)
 17     # print(border,rows-border,border,cols-border)
 18     temp=[]
 19     for i in range(border,rows-border):
 20         for j in range(border,cols-border):
 21             temp=src[i-border:i+border+1,j-border:j+border+1]  # 從圖像獲取與核匹配的圖像
 22             # 切片語法：索引位置包括開頭但不包括結尾 [start: end: step]
 23             dst[i][j]=(kernel*temp).sum()  # 計算卷積
 24     return dst
 25 
 26 
 27 # step0:讀入圖片
 28 img=cv2.imread("images/luxun.png",cv2.IMREAD_GRAYSCALE)  # 讀入圖片
 29 
 30 # step1:高斯濾波
 31 img=cv2.GaussianBlur(img,(5,5),0)
 32 
 33 # step2:計算梯度強度和方向
 34 sobel_x = np.array([[-1, 0, 1],[-2,0,+2],[-1, 0, 1]])   # sobel的x方向算子
 35 sobel_y = np.array([[1, 2, 1],[0,0,0],[-1, -2, -1]])    # sobel的y方向算子
 36 
 37 # img_x=cv2.filter2D(img,-1, sobel_x)  # 這個濾波器會將卷積結果歸一化到0-255，沒法計算梯度方向。
 38 # img_y=cv2.filter2D(img,-1, sobel_y)  # 而真正的圖像卷積可能會出現負數，所以只能本身寫個卷積。
 39 img_x=conv2d(img,sobel_x)  # 使用我本身的寫的卷積計算梯度
 40 img_y=conv2d(img,sobel_y)
 41 G=np.sqrt(img_x*img_x+img_y*img_y)  # 梯度幅值
 42 theta=np.arctan(img_y,(img_x+0.0000000000001))*180/np.pi  # 化爲角度，分母+極小值是爲了不除以0
 43 
 44 # plt.imshow(theta, 'gray'),plt.title('梯度相位')
 45 # plt.imshow(G, 'gray'),plt.title('梯度幅值')
 46 # plt.show()
 47 # exit()
 48 
 49 # step3:非極大值抑制
 50 anchor=np.where(G!=0)  # 獲取非零梯度的位置
 51 Mask=np.zeros(img.shape)
 52 for i in range(len(anchor[0])):
 53     x=anchor[0][i]  # 取出第i個非零梯度的x座標
 54     y=anchor[1][i]
 55     center_point=G[x,y]
 56     current_theta=theta[x,y]
 57     dTmp1=0
 58     dTmp2=0
 59     W=0
 60     if current_theta>=0 and current_theta<45:
 61         #        g1         第一種狀況
 62         # g4  C  g2
 63         # g3
 64         g1 = G[x + 1, y - 1]
 65         g2 = G[x + 1, y]
 66         g3 = G[x - 1, y + 1]
 67         g4 = G[x - 1, y]
 68         W=abs(np.tan(current_theta*np.pi/180))  # tan0-45範圍爲0-1
 69         dTmp1= W*g1+(1-W)*g2
 70         dTmp2= W*g3+(1-W)*g4
 71 
 72     elif current_theta>=45 and current_theta<90:
 73         #     g2 g1         第二種狀況
 74         #     C
 75         # g3  g4
 76         g1 = G[x + 1, y - 1]
 77         g2 = G[x, y - 1]
 78         g3 = G[x - 1, y + 1]
 79         g4 = G[x, y + 1]
 80         W = abs(np.tan((current_theta-90) * np.pi / 180))
 81         dTmp1= W*g1+(1-W)*g2
 82         dTmp2= W*g3+(1-W)*g4
 83     elif current_theta>=-90 and current_theta<-45:
 84         # g1  g2            第三種狀況
 85         #     C
 86         #     g4 g3
 87         g1 = G[x - 1, y - 1]
 88         g2 = G[x, y - 1]
 89         g3 = G[x + 1, y + 1]
 90         g4 = G[x, y + 1]
 91         W = abs(np.tan((current_theta-90) * np.pi / 180))
 92         dTmp1= W*g1+(1-W)*g2
 93         dTmp2= W*g3+(1-W)*g4
 94     elif current_theta>=-45 and current_theta<0:
 95         # g3                第四種狀況
 96         # g4  C  g2
 97         #        g1
 98         g1 = G[x + 1, y + 1]
 99         g2 = G[x + 1, y]
100         g3 = G[x - 1, y - 1]
101         g4 = G[x - 1, y]
102         W = abs(np.tan(current_theta * np.pi / 180))
103         dTmp1= W*g1+(1-W)*g2
104         dTmp2= W*g3+(1-W)*g4
105 
106     if dTmp1<center_point and dTmp2<center_point:   # 記錄極大值結果
107             Mask[x,y]=center_point
108 # plt.imshow(Mask,'gray'),plt.title('Mask')
109 # plt.show()
110 # exit()
111 
112 # step4:雙閾值選取
113 high_threshold=100
114 low_threshold=high_threshold/2
115 strong_edge=np.zeros(G.shape)   # 強邊緣
116 weak_edge=np.zeros(G.shape)     # 弱邊緣
117 
118 xNum = [1, 1, 0, -1, -1, -1, 0, 1]  # 8鄰域偏移座標
119 yNum = [0, 1, 1, 1, 0, -1, -1, -1]
120 [rows, cols] = G.shape
121 for i in range(rows):
122     for j in range(cols):
123         current_point=Mask[i,j]    
124         if current_point>0:
125             if current_point>high_threshold:   # 強邊緣提取
126                 strong_edge[i,j]=255
127             elif current_point<high_threshold and current_point>low_threshold:  # 弱邊緣提取
128                 # step6:順便進行邊緣鏈接
129                 change = True
130                 while change:
131                     change = False
132                     for k in range(8):
133                         xx=i+xNum[k]
134                         yy=j+yNum[k]
135                         if Mask[xx,yy]>high_threshold:
136                             weak_edge[i, j] = 255
137                             break  # 跳出八鄰域循環
138 output=strong_edge+weak_edge    # 強弱邊緣綜合效果
139 
140 img_edge = cv2.Canny(img, 50, 100)  # opencv實現效果
141 
142 # 顯示效果
143 plt.subplot(141), plt.imshow(strong_edge, 'gray'),plt.title('strong_edge')
144 plt.subplot(142), plt.imshow(weak_edge, 'gray'),plt.title('weak_edge')
145 plt.subplot(143), plt.imshow(output, 'gray'),plt.title('result')
146 plt.subplot(144), plt.imshow(img_edge, 'gray'),plt.title('opencv')
147 plt.show()