边缘检测
一、实验原理（及部分代码贴图）
图像边缘信息主要集中在高频段，通常说图像锐化或检测边缘，实质就是高频滤波。我们知道微分运算是求信号的变化率，具有加强高频分量的作用。在空域运算中来说，对图像的锐化就是计算微分。由于数字图像的离散信号，微分运算就变成计算差分或梯度。
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()