点击查看代码
from PIL import Image
from matplotlib import pyplot as plt
import math
from PIL import ImageFilter
#Image转为list类型
im = Image.open("pic.jpg")
im_gray = im.convert('L') #rgb转灰度图
#im_gray.show() #显示灰度图
#im_gray =im_gray.filter(ImageFilter.SMOOTH)
#im_gray.show()
im_array = im_gray.load()#从 PIL 灰度图像`im_gray`获取像素访问对象`im_array`，可以用`[x,y]`直接读取每个位置灰度值。注意：`x`是水平坐标（宽），`y`是垂直坐标（高）。
im_list =[[0 for i in range(im.size[0])] for j in range(im.size[1])]
for i in range(im.size[1]):
for j in range(im.size[0]):
im_list[i][j] = im_array[j,i]
plt.figure(1)#创建编号为 1 的画布窗口，用来绘图。
plt.imshow(im_list,'gray')
#plt.show()
#求梯度（剃度模值和方向）
grad = [[[0,0]for i in range(im.size[0])] for j in range(im.size[1])]
#Sobel 算子是 **3×3 窗口**，计算当前像素梯度时，**需要用到当前像素上下左右总共 9 个邻域像素**。
#这就是 **Sobel-X 3×3 卷积模板**，`右边一列像素求和 − 左边一列像素求和`，中间列不参与。
#Sobel 算子中间行（SobelX）/ 中间列（SobelY）×2，是**增加中心附近像素的权重，做平滑降噪**。
for j in range(1,im.size[1] -1):
for i in range(1,im.size[0] -1):
#x方向梯度
grad_x = im_list[j + 1][i + 1] + im_list[j - 1][i + 1] + 2 * im_list[j][i + 1] - \
im_list[j - 1][i - 1] - im_list[j + 1][i - 1] - 2 * im_list[j][i - 1]
#y方向梯度
grad_y = im_list[j + 1][i - 1] + im_list[j + 1][i + 1] + 2 * im_list[j + 1][i] - \
im_list[j - 1][i - 1] - im_list[j - 1][i + 1] - 2 * im_list[j - 1][i]
grad_x = math.floor(grad_x / 4)
grad_y = math.floor(grad_y / 4)
#合梯度
grad[j][i][0] =math.floor(math.sqrt(grad_x*grad_x + grad_y*grad_y))
if(grad[j][i][0] > 255):
grad[j][i][0] = 255
# if(grad[j][i][0] < 50):
# grad[j][i][0] = 0
#梯度方向
# ========== 梯度方向量化 ==========
# 梯度方向 θ = atan2(Gy, Gx)，返回弧度，范围 [-π ~ π]
# atan2(y,x)：由y、x求角度，图像坐标系y轴向下
if(grad_x == 0):
grad[j][i][1] = 2 # Gx=0 → 只有y方向梯度，90度
else:
theta = math.atan2(grad_y, grad_x)
# π/8 = 22.5°；3π/8=67.5°
if(math.fabs(theta) < math.pi / 8):
grad[j][i][1] = 0 # |θ|<22.5°，接近水平0度
elif(theta > 0):
if(math.fabs(theta) < math.pi * 3 / 8):
grad[j][i][1] = 1 # 22.5°<θ<67.5° →45度
else:
grad[j][i][1] = 2 # θ>67.5° →90度
else:
if (math.fabs(theta) < math.pi * 3 / 8):
grad[j][i][1] = 3 # -67.5°<θ<-22.5° →135度
else:
grad[j][i][1] = 2 # θ<-67.5° →90度
#显示梯度图
img = [[0 for i in range(im.size[0])] for j in range(im.size[1])]
for i in range(im.size[1]):
for j in range(im.size[0]):
img[i][j] = grad[i][j][0]
plt.figure(2)
plt.imshow(img, cmap='gray')
plt.show()
#非极大值抑制,使边缘更清晰和细
#不从 0 开始，不取最后一行 / 列，**跳过图像边界**，防止索引越界，访问不到邻居像素。
img2 = img
for j in range(1, im.size[1] - 1):
for i in range(1, im.size[0] - 1):
dir = grad[j][i][1]
grad_now = grad[j][i][0]
if (dir == 0): #梯度方向为x
if(grad_now <grad[j][i + 1][0] or grad_now < grad[j][i - 1][0]):
# grad[j][i][0] == 0
img2[j][i] = 0
print('0')
elif(dir == 1):#45度方向
if(grad_now < grad[j + 1][i + 1][0] or grad_now < grad[j - 1][i - 1][0]):
# grad[j][i][0] == 0
img2[j][i] = 0
print('1')
elif(dir == 2): #y方向
if(grad_now < grad[j + 1][i][0] or grad_now < grad[j - 1][i][0]):
# grad[j][i][0] == 0
img2[j][i] = 0
print('2')
else: #145度
if(grad_now< grad[j + 1][i - 1][0] or grad_now < grad[j - 1][i + 1][0]):
# grad[j][i][0] == 0
img2[j][i] = 0
print('3')
#显示非极大值抑制后的图像
plt.figure(3)
plt.imshow(img2, cmap='gray')
plt.show()
#双阈值检测
low = 30
high = 60
for i in range(im.size[1]):
for j in range(im.size[0]):
if(img2[i][j] < low):
img2[i][j] = 0
elif(img2[i][j] > high):
img2[i][j] = 255
#双阈值检测后的图像
plt.figure(4)
plt.imshow(img2, cmap='gray')
plt.show()
#抑制孤立的弱边缘
#含自身 (0,0)，不影响逻辑，也可以去掉。
listx = [-1, 0, 1]
listy = [-1, 0, 1]
#x、y 偏移量，组合起来就是当前像素周围**8 个相邻像素**（8 邻域）：
#(-1,-1),(-1,0),(-1,1),
#(0,-1), (0,0), (0,1),
#(1,-1), (1,0), (1,1)
img3 = img2
for j in range(1, im.size[1] - 1):
for i in range(1, im.size[0] - 1):
flag = 1
# flag 标记：
#flag=1：8 邻域**没有强边缘 255**，这个弱边缘是孤立的，需要删掉
# flag=0：8 邻域**存在强边缘 255**，这个弱边缘和主边缘相连，保留
for dx in range(len(listx)):
for dy in range(len(listy)):
#下标索引
#通过下标去取出偏移
#`listx = [-1, 0, 1]`，`len(listx)=3`
#循环组合：dx=0,1,2；dy=0,1,2，一共 3×3=9 次循环，遍历当前像素和周围 8 邻域（包含自己本身）。
nj = j + listy[dy]
ni = i + listx[dx]
if(img2[j][i] == 255):
flag = 0
break
if(not flag):
break
if(flag):
img3[j][i] = 0
#抑制孤立的弱边缘后的图像
# print(&#39;ssssssss&#39;)
plt.figure(5)
plt.imshow(img3, cmap='gray')
plt.show()
浙公网安备 33010602011771号