邊緣檢測
一、實驗原理（及部分代碼貼圖）
圖像邊緣信息主要集中在高頻段，通常說圖像銳化或檢測邊緣，實質就是高頻濾波。我們知道微分運算是求信號的變化率，具有加強高頻分量的作用。在空域運算中來說，對圖像的銳化就是計算微分。由於數字圖像的離散信號，微分運算就變成計算差分或梯度。
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()