在这里插入代码片
```import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import scipy.signal as signal
import cv2 as cv
import random
import easygui as g
import imutils
import time
import math
def Canny(img):
# Gray scale
def BGR2GRAY(img):
b = img[:, :, 0].copy()
g = img[:, :, 1].copy()
r = img[:, :, 2].copy()
# Gray scale
out = 0.2126 * r + 0.7152 * g + 0.0722 * b
out = out.astype(np.uint8)
return out
# Gaussian filter for grayscale
def gaussian_filter(img, K_size=3, sigma=1.3):
if len(img.shape) == 3:
H, W, C = img.shape
gray = False
else:
img = np.expand_dims(img, axis=-1)
H, W, C = img.shape
gray = True
## Zero padding
pad = K_size // 2
out = np.zeros([H + pad * 2, W + pad * 2, C], dtype=np.float)
out[pad: pad + H, pad: pad + W] = img.copy().astype(np.float)
## prepare Kernel
K = np.zeros((K_size, K_size), dtype=np.float)
for x in range(-pad, -pad + K_size):
for y in range(-pad, -pad + K_size):
K[y + pad, x + pad] = np.exp(- (x ** 2 + y ** 2) / (2 * sigma * sigma))
# K /= (sigma * np.sqrt(2 * np.pi))
K /= (2 * np.pi * sigma * sigma)
K /= K.sum() #得到了我们自己的高斯模板 k
tmp = out.copy()
# filtering
for y in range(H):
for x in range(W):
for c in range(C):
out[pad + y, pad + x, c] = np.sum(K * tmp[y: y + K_size, x: x + K_size, c])
out = np.clip(out, 0, 255)
out = out[pad: pad + H, pad: pad + W] #原图经过pading操作后比原图大pading,我们这里裁减掉padding,得到图片原始的真实大小
out = out.astype(np.uint8)
if gray:
out = out[..., 0] #如果是灰度图,我们将其降维成两维,因为之前将图像的维度上升到了三维
return out
# sobel filter
def sobel_filter(img, K_size=3):
if len(img.shape) == 3:
H, W, C = img.shape
else:
H, W = img.shape
# Zero padding
pad = K_size // 2
out = np.zeros((H + pad * 2, W + pad * 2), dtype=np.float)
out[pad: pad + H, pad: pad + W] = img.copy().astype(np.float)
tmp = out.copy()
out_v = out.copy()
out_h = out.copy()
## Sobel vertical
Kv = [[1., 2., 1.], [0., 0., 0.], [-1., -2., -1.]]
## Sobel horizontal
Kh = [[1., 0., -1.], [2., 0., -2.], [1., 0., -1.]]
# filtering
for y in range(H):
for x in range(W):
out_v[pad + y, pad + x] = np.sum(Kv * (tmp[y: y + K_size, x: x + K_size]))
out_h[pad + y, pad + x] = np.sum(Kh * (tmp[y: y + K_size, x: x + K_size]))
out_v = np.clip(out_v, 0, 255)
out_h = np.clip(out_h, 0, 255)
out_v = out_v[pad: pad + H, pad: pad + W]
out_v = out_v.astype(np.uint8)
out_h = out_h[pad: pad + H, pad: pad + W]
out_h = out_h.astype(np.uint8)
return out_v, out_h #输出水平方向与竖直方向的sober算子的输出矩阵
def get_edge_angle(fx, fy):
# get edge strength
edge = np.sqrt(np.power(fx.astype(np.float32), 2) + np.power(fy.astype(np.float32), 2))
edge = np.clip(edge, 0, 255)
fx = np.maximum(fx, 1e-10) #maximum输出最大值,这里为了防止我们输出一个零的数值,因此这里特地把0设置为0.00000001
# fx[np.abs(fx) <= 1e-5] = 1e-5
# get edge angle
angle = np.arctan(fy / fx) #这里如果x=0,我们可以认为垂直的导数无限大,因此我们这里除以0.0000001没毛病
return edge, angle #返回边界图的矩阵以及相应的角度矩阵
def angle_quantization(angle):
angle = angle / np.pi * 180 #将角度转化弧度制
angle[angle < -22.5] = 180 + angle[angle < -22.5]
_angle = np.zeros_like(angle, dtype=np.uint8)
_angle[np.where(angle <= 22.5)] = 0 #这里我们设置四个阈值,使得角度判断更正为对0,,45,90,135的判断
_angle[np.where((angle > 22.5) & (angle <= 67.5))] = 45
_angle[np.where((angle > 67.5) & (angle <= 112.5))] = 90
_angle[np.where((angle > 112.5) & (angle <= 157.5))] = 135 #以45度作为一次划分
return _angle
def non_maximum_suppression(angle, edge): #非极大值抑制
H, W = angle.shape
_edge = edge.copy()
for y in range(H):
for x in range(W):
if angle[y, x] == 0:
dx1, dy1, dx2, dy2 = -1, 0, 1, 0 #因为无论哪个方向都有两种定义方式,因此设置两种的梯度方向信息
elif angle[y, x] == 45:
dx1, dy1, dx2, dy2 = 1, 1, -1, -1
elif angle[y, x] == 90:
dx1, dy1, dx2, dy2 = 0, -1, 0, 1
elif angle[y, x] == 135:
dx1, dy1, dx2, dy2 = -1, 1, 1, -1
if x == 0:
dx1 = max(dx1, 0)
dx2 = max(dx2, 0) #这些步骤的作用是避免我的图像边界点加上梯度信息后超出了边界范围,让最小值为0
if x == W - 1:
dx1 = min(dx1, 0)
dx2 = min(dx2, 0)
if y == 0:
dy1 = max(dy1, 0)
dy2 = max(dy2, 0)
if y == H - 1:
dy1 = min(dy1, 0)
dy2 = min(dy2, 0)
if max(max(edge[y, x], edge[y + dy1, x + dx1]), edge[y + dy2, x + dx2]) != edge[y, x]: #根据定义的点的信息来判断当前点与周围点是否为极大值的点,如果不是就舍弃掉。
_edge[y, x] = 0
return _edge
def hysterisis(edge, HT=100, LT=30):
H, W = edge.shape
# Histeresis threshold
edge[edge >= HT] = 255
edge[edge <= LT] = 0 #进行我们自己的二值化处理
_edge = np.zeros((H + 2, W + 2), dtype=np.float32) #因为我们自己设置的图片是经过了3*3的padding,因此每个人边界平均长了2
_edge[1: H + 1, 1: W + 1] = edge
## 8 - Nearest neighbor
nn = np.array(((1., 1., 1.), (1., 0., 1.), (1., 1., 1.)), dtype=np.float32)
for y in range(1, H + 2):
for x in range(1, W + 2):
if _edge[y, x] < LT or _edge[y, x] > HT:
continue
if np.max(_edge[y - 1:y + 2, x - 1:x + 2] * nn) >= HT:
_edge[y, x] = 255 #膨胀一下
else:
_edge[y, x] = 0
edge = _edge[1:H + 1, 1:W + 1] #裁掉多余的部分
return edge
# grayscale
gray = BGR2GRAY(img) #灰度化
# gaussian filtering
gaussian = gaussian_filter(gray, K_size=5, sigma=1.4) #进行高斯处理,高斯的模板为5*5,sigma为1.4(平滑)
# sobel filtering
fy, fx = sobel_filter(gaussian, K_size=3) #进行sobel处理,返回水平方向与垂直方向两个矩阵(提取边缘信息)
# get edge strength, angle
edge, angle = get_edge_angle(fx, fy) #根据水平和竖直方向的矩阵返回边界和角度,此时的角度还为弧度制
# angle quantization
angle = angle_quantization(angle) #将角度进行转化,转化为多少多少度的形式
# non maximum suppression
edge = non_maximum_suppression(angle, edge) #进行非极大值抑制,将每个像素点根据自身以及自己的梯度信息将周围邻域进行比较, 如果不是局部最大值,舍弃掉
#相当于保留最强的边界信息
# hysterisis threshold
out = hysterisis(edge, 100, 30) #二值化
return out
# 霍夫变换实现检测图像中的20条直线
def Hough_Line(edge, img):
## Voting
def voting(edge):
H, W = edge.shape
drho = 1
dtheta = 1
# get rho max length
rho_max = np.ceil(np.sqrt(H ** 2 + W ** 2)).astype(np.int) #ceil为向上取整
# hough table
hough = np.zeros((rho_max, 180), dtype=np.int) #建立hough映射表
# get index of edge
# ind[0] 是 符合条件的纵坐标,ind[1]是符合条件的横坐标
ind = np.where(edge == 255) #找到符合条件的、满足非极大值抑制的实际点
## hough transformation
# zip函数返回元组
for y, x in zip(ind[0], ind[1]):
for theta in range(0, 180, dtheta):
# get polar coordinat4s
t = np.pi / 180 * theta #将角度在转化为弧度制,方便下面公式的计算
rho = int(x * np.cos(t) + y * np.sin(t)) #计算hough空间的纵轴的值
# vote
hough[rho, theta] += 1 #对于符合条件的hough空间进行统计,满足条件个数加一
out = hough.astype(np.uint8)
return out
# non maximum suppression
def non_maximum_suppression(hough):
rho_max, _ = hough.shape
## non maximum suppression
for y in range(rho_max):
for x in range(180):
# get 8 nearest neighbor
x1 = max(x - 1, 0)
x2 = min(x + 2, 180)
y1 = max(y - 1, 0)
y2 = min(y + 2, rho_max - 1)
if np.max(hough[y1:y2, x1:x2]) == hough[y, x] and hough[y, x] != 0: #在hough空间中也进行一次非最大值抑制
pass
# hough[y,x] = 255
else:
hough[y, x] = 0 #注意:hough空间中存储的是[y,x]的个数
return hough
def inverse_hough(hough, img):
H, W, _ = img.shape
rho_max, _ = hough.shape
out = img.copy()
# get x, y index of hough table
# np.ravel 将多维数组降为1维
# argsort 将数组元素从小到大排序,返回索引
# [::-1] 反序->从大到小
# [:20] 前20个
ind_x = np.argsort(hough.ravel())[::-1][:20] #取hough空间中最大的20个的索引作为我们的目标
ind_y = ind_x.copy()
thetas = ind_x % 180 #因为原数据为(最大可能边长*180),因此我们对180取余数即可得到我们的角度值的大小
rhos = ind_y // 180 #因为原数据为(最大可能边长*180)。我们除以180,向下取整,即可得到我们的hough的rhos
# each theta and rho
for theta, rho in zip(thetas, rhos): #通过zip得到每个点的角度和长度
# theta[radian] -> angle[degree]
t = np.pi / 180. * theta #转化为弧度制
line1 = 0
line2 = 0
arctan1 = 0
arectan2 = 0
x1=[]
y1=[]
arctan1_y = 0
line1_end = 0
# hough -> (x,y)
for x in range(W):
if np.sin(t) != 0:
y = - (np.cos(t) / np.sin(t)) * x + (rho) / np.sin(t)
y = int(y)
if y >= H or y < 0:
continue
out[y, x] = [0, 255, 255]
line1 +=1
if len(x1) <= 5:
x1.append((x,y))
if (x1[1][0]-x1[0][0])==0 or (x1[1][1]-x1[0][1])==0:
line1_end = line1
else:
arctan1=(x1[1][1]-x1[0][1])/(x1[1][0]-x1[0][0])
arctan1_y = arctan1 * line1
line1_end = math.sqrt(arctan1_y**2+line1**2)
print('定x的线的长度为 %d' % line1_end)
for y in range(H):
if np.cos(t) != 0:
x = - (np.sin(t) / np.cos(t)) * y + (rho) / np.cos(t)
x = int(x)
if x >= W or x < 0:
continue
out[y, x] = [0, 0, 255]
if len(y1) <= 5:
y1.append((x,y))
line2 +=1
if (y1[1][0] - y1[0][0]) ==0 or (y1[1][1] - y1[0][1])==0:
line2_end = line2
else:
arctan2 = (y1[1][1] - y1[0][1]) / (y1[1][0] - y1[0][0])
arectan2_y = arectan2 * line2
line2_end = math.sqrt(arectan2_y**2+line2**2)
print('定y的线的长度为 %d' % line2_end)
out = out.astype(np.uint8)
return out
# voting
hough = voting(edge)
# non maximum suppression
hough = non_maximum_suppression(hough)
# inverse hough
out = inverse_hough(hough, img)
return out
def otsu_binarization(img, th=128):
H, W = img.shape
out = img.copy()
max_sigma = 0
max_t = 0
# determine threshold
for _t in range(1, 255):
v0 = out[np.where(out < _t)]
m0 = np.mean(v0) if len(v0) > 0 else 0.
w0 = len(v0) / (H * W)
v1 = out[np.where(out >= _t)]
m1 = np.mean(v1) if len(v1) > 0 else 0.
w1 = len(v1) / (H * W)
sigma = w0 * w1 * ((m0 - m1) ** 2)
if sigma > max_sigma:
max_sigma = sigma
max_t = _t
# Binarization
print("threshold >>", max_t)
th = max_t
out[out < th] = 0
out[out >= th] = 255
return out
# Morphology Dilate
def Morphology_Dilate(img, Dil_time=1):
H, W = img.shape
out = img.copy()
# kernel
# MF = np.array(((0, 1, 0),
# (1, 0, 1),
# (0, 1, 0)), dtype=np.int)
MF = np.ones((11, 11))
# each erode
for i in range(Dil_time):
tmp = np.pad(out, (5, 5), 'edge')
# erode
for y in range(5, H):
for x in range(5, W):
if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) >255:
out[y, x] = 255
out_gap = out - img
out_gap = out_gap.astype(np.uint8) #膨胀图减去二值图
cv.imshow('pengzhang - erzhitu',out_gap)
cv.imshow('dilate',out)
def Morphology_Dilate1(img, Dil_time=1):
H, W = img.shape
out = img.copy()
# kernel
# MF = np.array(((0, 1, 0),
# (1, 0, 1),
# (0, 1, 0)), dtype=np.int)
MF = np.ones((11, 11))
# each erode
for i in range(Dil_time):
tmp = np.pad(out, (5, 5), 'edge')
# erode
for y in range(5, H):
for x in range(5, W):
if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) > 255:
out[y, x] = 255
return out
# Morphology Erode
def Morphology_Erode(img, Erode_time=1):
H, W = img.shape
out = img.copy()
# kernel
# MF = np.array(((0, 1, 0),
# (1, 0, 1),
# (0, 1, 0)), dtype=np.int)
MF = np.ones((11,11))
# each erode
for i in range(Erode_time):
tmp = np.pad(out, (5, 5), 'edge')
# erode
for y in range(5, H):
for x in range(5, W):
if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) < 255 * 121:
out[y, x] = 0
out_gap = img-out
out_gap = out_gap.astype(np.uint8) #原图减去复试图
cv.imshow('erode',out)
cv.imshow('yuan tu - fushi ',out_gap)
cv.waitKey(0)
def Morphology_Erode1(img, Erode_time=1):
H, W = img.shape
out = img.copy()
# kernel
# MF = np.array(((0, 1, 0),
# (1, 0, 1),
# (0, 1, 0)), dtype=np.int)
MF = np.ones((11, 11))
# each erode
for i in range(Erode_time):
tmp = np.pad(out, (5, 5), 'edge')
# erode
for y in range(5, H):
for x in range(5, W):
if np.sum(MF * tmp[y - 5:y + 6, x - 5:x + 6]) < 255 * 121:
out[y, x] = 0
return out
def hough_fun(img):
edge = Canny(img)
# Hough
out = Hough_Line(edge, img)
out = out.astype(np.uint8)
cv.imshow('hough',out)
cv.waitKey(0)
def open_and_close(img, Erode_time=1):
img = img[:,:,0]*0.11 + img[:,:,1] *0.59+img[:,:,2] *0.3
img= img.astype(np.uint8)
img = otsu_binarization(img)
H, W = img.shape
out = img.copy()
out_open = Morphology_Erode1(out) #开运算先腐蚀后膨胀
out_open= Morphology_Dilate1(out_open)
out_close = Morphology_Dilate1(out)
out_close = Morphology_Erode1(out_close)
#闭运算先膨胀后腐蚀
cv.imshow('open',out_open)
cv.imshow('close',out_close)
cv.waitKey(0)
if __name__ == '__main__':
msg = "请输入您想要完成的任务(建议您第一步先打开图片)"
title = '第四次作业'
choice = ('打开图片', '退出')
a = g.buttonbox(msg=msg, title=title, choices=choice)
if a == '打开图片':
filename = g.fileopenbox(msg="请打开一个jpg文件")
img = cv.imread(filename)
msg1 = "选择您想要实现的功能"
title1 = '第四次作业'
choice1 = ('图像腐蚀', '直线检测', '数学形态应用','显示原图','显示灰度图','重新选择图片','退出')
q = 1
while q:
b = g.buttonbox(msg=msg1, title=title1, choices=choice1)
if b == '图像腐蚀':
img_gray = img[:, :, 0] * 0.11 + img[:, :, 1] * 0.59 + img[:, :, 2] * 0.3
img_gray = img_gray.astype(np.uint8)
ostu=otsu_binarization(img_gray)
cv.imshow('binary image',ostu)
Morphology_Dilate(ostu)
Morphology_Erode(ostu)
elif b == '直线检测':
hough_fun(img)
elif b == '数学形态应用':
open_and_close(img)
elif b == '显示原图':
cv.imshow('原图',img)
elif b == '显示灰度图':
img_gray1 = img[:, :, 0] * 0.11 + img[:, :, 1] * 0.59 + img[:, :, 2] * 0.3
img_gray1 = img_gray1.astype(np.uint8)
cv.imshow('原图',img_gray1)
elif b == '重新选择图片':
filename = g.fileopenbox(msg="请打开一个jpg文件")
img = cv.imread(filename)
else:
q = 0
opencv第四次作业(图像腐蚀,直线检测,数学形态应用)
最新推荐文章于 2024-09-18 20:39:31 发布