影象的邊緣檢測（canny、sobel、prewitt的比較）
邊緣檢測
一、實驗原理（及部分程式碼貼圖）
影象邊緣資訊主要集中在高頻段，通常說影象銳化或檢測邊緣，實質就是高頻濾波。我們知道微分運算是求訊號的變化率，具有加強高頻分量的作用。在空域運算中來說，對影象的銳化就是計算微分。由於數字影象的離散訊號，微分運算就變成計算差分或梯度。
Canny實現運算元流程
1.高斯平滑
類似於LoG運算元作高斯模糊一樣，主要作用就是去除噪聲。因為噪聲也集中於高頻訊號，很容易被識別為偽邊緣。應用高斯模糊去除噪聲，降低偽邊緣的識別。但是由於影象邊緣資訊也是高頻訊號，高斯模糊的半徑選擇很重要，過大的半徑很容易讓一些弱邊緣檢測不到。
2.計算梯度幅值和方向：
影象的邊緣可以指向不同方向，因此經典Canny演算法用了4個梯度運算元來分別計算水平，垂直和對角線方向的梯度。但是通常都不用四個梯度運算元來分別計算四個方向。常用的邊緣差分運算元（如Rober，Prewitt，Sobel）計算水平和垂直方向的差分Gx和Gy。這樣就可以如下計算梯度模和方向。本次我選擇選擇Sobel運算元計算梯度。
3：非極大值抑制
沿著梯度方向進行非極大值抑制，第一步先算出梯度的朝向，然後比較當前畫素點的當前方向的梯度值在3X3區域的該方向上是否為最大值。為了更好的比較梯度的大小，首先，假設沿垂直方向和水平方向的相鄰畫素之間梯度變化是連續的，然後，比較該畫素點沿X和Y方向的梯度大小的方向，如果|dy|>|dx|且dy與dx方向相同，設權值w=|dx|/|dy|，並對相關鄰域該方向上的兩對點進行加權。最後將該點的梯度值與另外兩對點的加權梯度值進行非最大值抑制操作。
4：雙閥值檢測和邊緣連線
設定兩個閥值，分別為TL和TH。其中大於TH的都被檢測為邊緣，而低於TL的都被檢測為非邊緣。對於中間的畫素點，如果與確定為邊緣的畫素點鄰接，則判定為邊緣；否則為非邊緣。
Sobel運算元實現流程
Sobel運算元利用一階差分的思想來創造運算元，並對中心點的梯度重點考慮。
其中中心點 f(x, y) 是重點考慮的，它的權重應該多一些，所以改進成下面這樣的
然後分別計算偏 x 方向的 dx，偏 y 方向的 dy，利用公式
求出最終邊緣影象。（或者絕對值相加）
Prewitt運算元實現流程
通常用 f '(x) = f(x + 1) - f(x - 1) 近似計算一階差分。可以提出係數：[-1, 0, 1]，這個就是模板。在二維情況下是：
最後，利用X方向和Y方向梯度的最大值或者均值和來計算最終邊緣矩陣圖。
二、實驗結果
從網路上下載影象
實驗1
對原圖分別用canny，sobel，prewitt進行邊緣提取，並將自己的程式碼實現與opencv庫的實現作比較。實驗結果：
實驗2
對原圖分別用canny，sobel，prewitt進行邊緣提取，並將自己的程式碼實現與opencv庫的實現作比較。實驗結果：
實驗3
對原圖分別用canny，sobel，prewitt進行邊緣提取，並將自己的程式碼實現與opencv庫的實現作比較。實驗結果：
三、實驗分析
Prewitt運算元是一種中規中矩的影象邊緣檢測的微分運算元，利用一階差分進行檢測。由於運算元採用 33 模板對區域內的畫素值進行計算，所以相比於2X2的模板，它在水平和垂直的梯度方向更加明顯、但是對於有噪聲的影象和細邊緣影象的檢測魯棒性較差，沒有考慮相鄰點的距離遠近對當前畫素點的影響，速度相對較快。
Sobel運算元在Prewitt運算元的基礎上增加了權重的概念，認為相鄰點的距離遠近對當前畫素點的影響是不同的，距離越近的畫素點對應當前畫素的影響越大，從而實現影象銳化並突出邊緣輪廓。Sobel運算元對邊緣定位不是很準確，影象的邊緣不止一個畫素，邊緣容易出現多畫素寬度；當對精度要求不是很高時，是一種較為常用的邊緣檢測方法。
Canny邊緣檢測運算元是一種多級檢測演算法。原理在第一節已經介紹。Canny運算元基於以下三個基本目標：
低錯誤率：所有邊緣都應被找到，並且不應有虛假響應。
邊緣點應被很好地定位：已定位的邊緣必須儘可能接近真實邊緣。也就是說，由檢測子標記為邊緣的一點和真實邊緣的中心距離最小。
單個邊緣點響應：對於每個真實的邊緣點，檢測子應只返回一個點。也就是說真實邊緣周圍的區域性最大數應是最小的。
相比於sobel和prewitt運算元檢測出的邊緣較為細緻，檢測精度較高，但是速度較慢。
四、實驗程式碼
import matplotlib.pyplot as plt
import numpy as np
import math
import cv2
def canny_my(gray):
sigma1 = sigma2 = 0.88
sum = 0
gaussian = np.zeros([5, 5])
# 生成二維高斯濾波矩陣
for i in range(5):
for j in range(5):
gaussian[i, j] = math.exp((-1 / (2 * sigma1 * sigma2)) * (np.square(i - 3)
+ np.square(j - 3))) / (
2 * math.pi * sigma1 * sigma2)
sum = sum + gaussian[i, j]
# 歸一化處理
gaussian = gaussian / sum
# print(gaussian)
# print("step1")
# step1.高斯濾波
W, H = gray.shape
new_gray = np.zeros([W - 5, H - 5])
for i in range(W - 5):
for j in range(H - 5):
new_gray[i, j] = np.sum(gray[i:i + 5, j:j + 5] * gaussian) # 與高斯矩陣卷積實現濾波
# plt.imshow(new_gray, cmap="gray")
# print("step2")
# step2.增強 通過求梯度幅值
W1, H1 = new_gray.shape
dx = np.zeros([W1 - 1, H1 - 1])
dy = np.zeros([W1 - 1, H1 - 1])
d = np.zeros([W1 - 1, H1 - 1])
for i in range(W1 - 1):
for j in range(H1 - 1):
dx[i, j] = new_gray[i, j + 1] - new_gray[i, j]
dy[i, j] = new_gray[i + 1, j] - new_gray[i, j]
d[i, j] = np.sqrt(np.square(dx[i, j]) + np.square(dy[i, j])) # 影象梯度幅值作為影象強度值
# plt.imshow(d, cmap="gray")
# setp3.非極大值抑制 NMS
W2, H2 = d.shape
NMS = np.copy(d)
# print("step3")
NMS[0, :] = NMS[W2 - 1, :] = NMS[:, 0] = NMS[:, H2 - 1] = 0
for i in range(1, W2 - 1):
for j in range(1, H2 - 1):
if d[i, j] == 0:
NMS[i, j] = 0
else:
gradX = dx[i, j]
gradY = dy[i, j]
gradTemp = d[i, j]
# 如果Y方向幅度值較大
if np.abs(gradY) > np.abs(gradX):
weight = np.abs(gradX) / np.abs(gradY)
grad2 = d[i - 1, j]
grad4 = d[i + 1, j]
# 如果x,y方向梯度符號相同
# C 為當前畫素，與g1-g4 的位置關係為：
# g1 g2
# C
# g4 g3
if gradX * gradY > 0:
grad1 = d[i - 1, j - 1]
grad3 = d[i + 1, j + 1]
# 如果x,y方向梯度符號相反
# g2 g1
# C
# g3 g4
else:
grad1 = d[i - 1, j + 1]
grad3 = d[i + 1, j - 1]
# 如果X方向幅度值較大
else:
weight = np.abs(gradY) / np.abs(gradX)
grad2 = d[i, j - 1]
grad4 = d[i, j + 1]
# 如果x,y方向梯度符號相同
# g3
# g2 C g4
# g1
if gradX * gradY > 0:
grad1 = d[i + 1, j - 1]
grad3 = d[i - 1, j + 1]
# 如果x,y方向梯度符號相反
# g1
# g2 C g4
# g3
else:
grad1 = d[i - 1, j - 1]
grad3 = d[i + 1, j + 1]
gradTemp1 = weight * grad1 + (1 - weight) * grad2
gradTemp2 = weight * grad3 + (1 - weight) * grad4
if gradTemp >= gradTemp1 and gradTemp >= gradTemp2:
NMS[i, j] = gradTemp
else:
NMS[i, j] = 0
# plt.imshow(NMS, cmap = "gray")
# print("step4")
# step4. 雙閾值演算法檢測、連線邊緣
W3, H3 = NMS.shape
DT = np.zeros([W3, H3])
# 定義高低閾值
TL = 0.1 * np.max(NMS)
TH = 0.3 * np.max(NMS)
for i in range(1, W3 - 1):
for j in range(1, H3 - 1):
if (NMS[i, j] < TL):
DT[i, j] = 0
elif (NMS[i, j] > TH):
DT[i, j] = 1
# 檢測相鄰元素是否有強畫素點
elif ((NMS[i - 1, j - 1:j + 1] < TH).any() or NMS[i + 1, j - 1:j + 1].any()
or (NMS[i, [j - 1, j + 1]] < TH).any()):
DT[i, j] = 1
return DT
def sobel_my(img):
W, H = img.shape
new_image = np.zeros((W - 3, H - 3))
new_imageX = np.zeros((W - 3, H - 3))
new_imageY = np.zeros((W - 3, H - 3))
s_suanziX = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]) # X方向
s_suanziY = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
for i in range(W - 3):
for j in range(H - 3):
new_imageX[i, j] = abs(np.sum(img[i:i + 3, j:j + 3] * s_suanziX))
new_imageY[i, j] = abs(np.sum(img[i:i + 3, j:j + 3] * s_suanziY))
new_image[i, j] = np.sqrt(np.square(new_imageX[i, j]) + np.square(new_imageY[i, j]))
return new_imageX, new_imageY, new_image # 無方向運算元處理的影象
def prewitt(img):
W, H = img.shape
new_image = np.zeros((W - 3, H - 3))
new_imageX = np.zeros((W - 3, H - 3))
new_imageY = np.zeros((W - 3, H - 3))
s_suanziX = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]]) # X方向
s_suanziY = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]])
for i in range(W - 3):
for j in range(H - 3):
new_imageX[i, j] = abs(np.sum(img[i:i + 3, j:j + 3] * s_suanziX))
new_imageY[i, j] = abs(np.sum(img[i:i + 3, j:j + 3] * s_suanziY))
new_image[i, j] = max(new_imageX[i, j], new_imageY[i, j])
return new_image
if __name__ == '__main__':
img = cv2.imread('22.png')
# cv2.imread讀取的圖片是以bgr的形式儲存，需要轉換成rgb格式
lena_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
# 灰度化處理
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
blur = cv2.GaussianBlur(lena_img, (5, 5), 0) # 用高斯濾波處理原影象降噪
canny = cv2.Canny(blur, 50, 150)
DT = canny_my(gray)
x = cv2.Sobel(gray, cv2.CV_16S, 1, 0) # 對x求一階導
y = cv2.Sobel(gray, cv2.CV_16S, 0, 1) # 對y求一階導
# Sobel函式求完導數後會有負值，還有會大於255的值。而原影象是uint8，即8位無符號數，所以Sobel建立的影象位數不夠，會有截斷。因此要使用16位有符號的資料型別，即cv2.CV_16S。處理完影象後，再使用cv2.convertScaleAbs()函式將其轉回原來的uint8格式，否則影象無法顯示。
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Sobel = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)
SX, SY, SM = sobel_my(gray)
kernelx = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]], dtype=int)
kernely = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], dtype=int)
Prewitt_x = cv2.filter2D(gray, cv2.CV_16S, kernelx)
Prewitt_y = cv2.filter2D(gray, cv2.CV_16S, kernely)
absPX = cv2.convertScaleAbs(Prewitt_x)
absPY = cv2.convertScaleAbs(Prewitt_y)
Prewitt = cv2.addWeighted(absPX, 0.5, absPY, 0.5, 0)
PW = prewitt(gray)
images = [gray, canny, DT, absX, absY, Sobel, SX, SY, SM, Prewitt, PW]
titles = ["Lena", "opencv Canny", "my Canny", "soblex", "sobley", "soble", "SobleX_my", "SobleY_my", "Soble_my",
"Prewitt", "Prewitt_my"]
plt.figure(figsize=(20, 20))
for i in range(len(images)):
plt.subplot(4, 3, i + 1)
plt.imshow(images[i], "gray")
plt.title(titles[i])
plt.show()