import numpy as np
import cv2
import matplotlib.pyplot as plt
def sobel(img, threshold):
# sobel边缘检测器：
# 计算两个方向的图像梯度
# 计算梯度的幅值
# 给定梯度幅值的阈值
G_x = np.array([[1, 0, -1],[2, 0, -2],[1, 0, -1]]) #垂直方向的检测
G_y = np.array([[1, 2, 1],[0, 0, 0],[-1, -2, -1]]) #水平方向的检测
# axis的值没有设定，返回矩阵的元素个数
# axis = 0，返回该二维矩阵的行数
# axis = 1，返回该二维矩阵的列数
rows = np.size(img, 0)
columns = np.size(img, 1)
mag = np.zeros(img.shape)
for i in range(0, rows - 2):
for j in range(0, columns - 2):
v = sum(sum(G_x * img[i:i+3, j:j+3])) # 垂直方向的检测
h = sum(sum(G_y * img[i:i+3, j:j+3])) # 水平方向的检测
mag[i+1, j+1] = np.sqrt((v ** 2) + (h ** 2)) #在计算每个像素的时候，将水平和垂直的值作平方和开根的处理
for p in range(0, rows):
for q in range(0, columns):
if mag[p, q] < threshold: # 给定梯度幅值的阈值进行检测
mag[p, q] = 0
return mag
def canny(img):
# 算法实现步骤
# 1. 应用高斯滤波来平滑(模糊)图像，目的是去除噪声
# 2. 计算滤波后图像的差分值/一阶梯度值
# 3. 计算各个像素点的梯度幅值以及梯度方向值
# 4. 应用非最大抑制技术NMS来消除边误检
# 5. 应用双阈值滞后阈值处理的方法来决定可能的（潜在的）边界
###### step1 首先进行高斯平滑 G[i, j] = (1/(2*pi*sigma**2))*exp(-((i-k-1)**2 + (j-k-1)**2)/2*sigma**2)
sigma = 1.4
length = 3 #卷积核的阶数设置为3阶
k = length // 2
gaussian = np.zeros([length, length])
for i in range(length):
for j in range(length):
gaussian[i, j] = np.exp(-((i-k) ** 2 + (j-k) ** 2) / (2 * sigma ** 2))
gaussian /= (2 * np.pi * sigma ** 2)
#gaussian = gaussian / np.sum(gaussian)
# 使用三阶的高斯滤波器进行平滑滤波得到新的图像
W, H = img.shape
new_image = np.zeros([W - k * 2, H - k * 2]) #行数、列数改变
for i in range(W - 2 * k):
for j in range(H - 2 * k):
# 卷积运算
new_image[i, j] = np.sum(img[i:i+length, j:j+length] * gaussian)
new_image = np.uint8(new_image)
###### step2. 计算滤波后图像的差分值/一阶梯度值 && step3. 计算各个像素点的梯度幅值以及梯度方向值
Gx = np.array([[1, 0, -1],[2, 0, -2],[1, 0, -1]]) #垂直方向的检测
Gy = np.array([[1, 2, 1],[0, 0, 0],[-1, -2, -1]]) #水平方向的检测
W, H = new_image.shape
gradients = np.zeros([W - 2, H - 2])
direction = np.zeros([W - 2, H - 2])
for i in range(W - 2):
for j in range(H - 2):
dx = np.sum(new_image[i:i+3, j:j+3] * Gx)
dy = np.sum(new_image[i:i+3, j:j+3] * Gy)
gradients[i, j] = np.sqrt(dx ** 2 + dy ** 2)
if dx == 0:
direction[i, j] = np.pi / 2
else:
direction[i, j] = np.arctan(dy / dx)
gradients = np.uint8(gradients)
##### step4. 应用非最大抑制技术NMS来消除边误检
#method and target
# 目的是将模糊（blurred）的边界变得清晰（sharp）。通俗的讲，就是保留了每个像素点上梯度强度的极大值，而删掉其他的值。对于每个像素点，进行如下操作：
# a) 将其梯度方向近似为以下值中的一个（0,45,90,135,180,225,270,315）（即上下左右和45度方向）
# b) 比较该像素点，和其梯度方向正负方向的像素点的梯度强度
# c) 如果该像素点梯度强度最大则保留，否则抑制（删除，即置为0）
W, H = gradients.shape
nms = np.copy(gradients[1:-1, 1:-1])
for i in range(1, W - 1):
for j in range(1, H - 1):
theta = direction[i, j]
weight = np.tan(theta)
if theta > np.pi / 4:
d1 = [0, 1]
d2 = [1, 1]
weight = 1 / weight
elif theta >= 0:
d1 = [1, 0]
d2 = [1, 1]
elif theta >= - np.pi / 4:
d1 = [1, 0]
d2 = [1, -1]
weight *= -1
else:
d1 = [0, -1]
d2 = [1, -1]
weight = -1 / weight
g1 = gradients[i + d1[0], j + d1[1]]
g2 = gradients[i + d2[0], j + d2[1]]
g3 = gradients[i - d1[0], j - d1[1]]
g4 = gradients[i - d2[0], j - d2[1]]
grade_count1 = g1 * weight + g2 * (1 - weight)
grade_count2 = g3 * weight + g4 * (1 - weight)
if grade_count1 > gradients[i, j] or grade_count2 > gradients[i, j]:
nms[i - 1, j - 1] = 0
##### step5. 应用双阈值滞后阈值处理的方法来决定可能的（潜在的）边界
#双阈值方法：设定一个阈值上界maxVal和阈值下界minVal，图像中的像素点如果大于阈值上界则认为必然是边界（称为强边界，strong edge），小于阈值下界则认为必然不是边界，
# 两者之间的则认为是候选项（称为弱边界，weak edge），需进行进一步处理——如果与确定为边缘的像素点邻接，则判定为边缘；否则为非边缘。
#滞后技术：和强边界相连的弱边界认为是边界，其他的弱边界则被抑制。
# 由真实边缘引起的弱边缘像素将连接到强边缘像素，而噪声响应未连接。
# 为了跟踪边缘连接，通过查看弱边缘像素及其8个邻域像素，只要其中一个为强边缘像素，则该弱边缘点就可以保留为真实的边缘。
###概括：要设置两个阈值，minVal 和 maxVal。梯度大于 maxVal 的任何边缘肯定是真边缘，而 minVal 以下的边缘肯定是非边缘，因此被丢弃。位于这两个阈值之间的边缘会基于其连通性而分类为边缘或非边缘，如果它们连接到“可靠边缘”像素，则它们被视为边缘的一部分。否则，也会被丢弃。
W, H = nms.shape
output = np.zeros([W, H])
# 定义高低阈值
TL = 0.1 * np.max(nms)
TH = 0.2 * np.max(nms)
for i in range(1, W-1):
for j in range(1, H-1):
# 双阈值选取
if (nms[i, j] < TL):
output[i, j] = 0
elif (nms[i, j] > TH):
output[i, j] = 1
# 滞后检测技术
elif (nms[i-1, j-1:j+1] < TH).any() or (nms[i+1, j-1:j+1] < TH).any()or (nms[i, [j-1, j+1]] < TH).any():
output[i, j] = 1
return output
def Log(img):
#算法步骤：
# step1 对图像做高斯滤波
# step2 求平滑后的图像的拉普拉斯（Laplacian）图
# step3 通过检测滤波结果的零交叉（Zero crossings）获得图像或物体的边缘
###代码思路：1.将高斯滤波器直接进行拉普拉斯运算，2.与初始图象卷积之后再进行零交叉验证
#part1
#method1 采用表达式进行构建
# sigma = 2
# length = 5 #卷积核的阶数设置为3阶
# k = length // 2
# LogMatrix = np.zeros([length, length])
# for i in range(length):
# for j in range(length):
# LogMatrix[i, j] = -(np.exp(-((i-k) ** 2 + (j-k) ** 2) / (2 * sigma ** 2)))*(2-((i-k) ** 2 + (j-k) ** 2)/(sigma ** 2))/(((2 * np.pi) ** 0.5) * sigma**3)
# LogMatrix = LogMatrix / np.sum(LogMatrix)
#method2 采用常规模板进行构建
LogMatrix = np.array([[0,0,-1,0,0],[0,-1,-2,-1,0],[-1,-2,16,-2,-1],[0,-1,-2,-1,0],[0,0,-1,0,0]]) #垂直方向的检测
W, H = img.shape
length = 5
k = 5 // 2
new_image = np.zeros([W - k * 2, H - k * 2])
for i in range(W - 2 * k):
for j in range(H - 2 * k):
# 卷积运算
new_image[i, j] = np.sum(img[i:i+length, j:j+length] * LogMatrix)
#part2 零交叉验证
# 以p为中心的一个3*3领域，p点处的零交叉意味着至少有两个相对的领域像素的符号不同。
# 在这个基础上还要判断中心点p的值介于两个符号不同像素的值之间，就是p的值要小于两个符号不同的像素的绝对值。
# 如果是，那才可以判断为零交叉点。
output = np.ones_like(new_image)
rows,columns = new_image.shape
for y in range(1,rows-1):
for x in range(1, columns - 1) :
if 200 > abs(new_image[y - 1, x]-new_image[y + 1, x]) and (new_image[y - 1, x]) * (new_image[y + 1, x]) < 0:
output[y, x] = 0
if 200 > abs(new_image[y, x - 1]-new_image[y, x + 1]) and (new_image[y, x - 1]) * (new_image[y, x + 1]) < 0 :
output[y, x] = 0
if 200 > abs(new_image[y + 1, x - 1]-new_image[y - 1, x + 1]) and (new_image[y + 1, x - 1]) * (new_image[y - 1, x + 1]) < 0 :
output[y, x] = 0
if 200 > abs(new_image[y - 1, x - 1]-new_image[y + 1, x + 1]) and (new_image[y - 1, x - 1]) * (new_image[y + 1, x + 1]) < 0 :
output[y, x] = 0
return output
if __name__ == '__main__':
img = cv2.imread('C:\\Users\\tracerX\\Desktop\\photo course\\lena.jpg', 0) # read an image
sobelImg = sobel(img, 100)
cannyImg = canny(img)
LogImage = Log(img)
plt.subplot(321) #在同一窗口2行2列的第一个窗口显示
plt.title("origin image")
plt.imshow(img,cmap=plt.cm.gray)
plt.subplot(322) #在同一窗口2行2列的第二个窗口显示
plt.title("sobel")
plt.imshow(sobelImg,cmap=plt.cm.gray)
plt.subplot(325) #在同一窗口2行2列的第三个窗口显示
plt.title("canny")
plt.imshow(cannyImg,cmap=plt.cm.gray)
plt.subplot(326) #在同一窗口2行2列的第四个窗口显示
plt.title("Log")
plt.imshow(LogImage,cmap=plt.cm.gray)
plt.show()