前言:霍夫变换于1962年被Hough提出用于模式识别领域(此处Hough提出的算法是用斜率(slope)和截距(intercept)进行转换,但是明显若直线斜率趋近于90°则是一个无穷大的数,计算无法表示),霍夫变换提取影像直线的技术于1972年被Richard O.Duda 和 Peter E.Hart两位科学家改进(Richard和Peter将其转换为引垂线的r、来表示并设置归入区间即比较好地解决了这一问题),其算法对于计算机视觉和图像处理做出了较大贡献,后续不断有学者在其算法基础上改进。本文从其论文出发,对原始霍夫变换进行实现,希望对网上部分该算法实现讲解进行修正与补充。
一、数学基础知识简概(并非完全与本次实现相关,但具有启发意义)
(一)平面几何
平面方程与直线方程的表示:
(二)极坐标系
定义:极坐标(二维坐标系统),其指在平面内取一个定点O,叫做极点,引一条射线Ox,叫做极轴,再选定一个长度单位和角度的正方向(通常取逆时针方向)。对于平面内任何一点M,用ρ(或r)表示OM线段的长度,表示从Ox到OM的角度,ρ叫做点M的极径,叫做点M的极角,有序数对(ρ,)就叫做点M的极坐标。
极坐标的唯一性:
此处,-r中的负号实际上表示的是对应极径的相反方向,其极角±180°即可换算为r。
直角坐标系与极坐标系的转换:
二、原理概述
(一)原理
首先,ρ平面(霍夫空间)中的细分数量决定这些点共线的精度,但是若数量过于细分有可能导致共线计数误差,再者计算次数与xy平面中的非背景点的数量n呈线性关系。其次,在Cartesian坐标系转霍夫空间(或直线方程引入和r表示)时,我并不赞同用极坐标系来解释这一说法,只能说思想近似。
(二)对比[1]与[2]的霍夫空间选取的方式:
其实,二者无论是0≤≤Π,还是-Π/2≤≤Π/2,其在ρ上的范围均为-R ≤ρ≤R,那为什么角可以这样取值呢?
对于而言,平面上一条直线(无限延长)总会与坐标轴相交,引入的目的是为了解决近乎垂直线的k无穷大问题,那么聚焦于与直线的关系,可以发现并不严格定义是引垂线与x正半轴的逆时针夹角,而是说其对应直线的斜率表示,那么每一条过原点的引垂线的旋转必然可以落入Π宽内(但是r值仍是独立的嗷)。如果还没理解请看下图:
如果我们要求ρ≥0,那么满足 0≤≤2Π,霍夫空间才能唯一确定。
三、代码部分
# 线特征提取
def lineEXTRACT(self):
img = cv2.merge(self.RGBbands)
# 图像灰度转换
# 判断输入影像是否为单波段
if len(self.RGBbands) == 1:
gray = img
else:
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
# 进行边缘检测提取出二值边缘图(sigma0高斯滤波),canny算子的最小判断和最大接受阈值是自己多次调试后的结果,其结果边缘信息比较丰富、准确
blur = cv2.GaussianBlur(gray, (3, 3), 0)
edge = cv2.Canny(blur, 25, 75)
# 构造投票箱与角度正弦、余弦查询表
# 在图像坐标系中应满足 r = sin(sita)* y + cos(sita)* x
# 我设置r的单位长度为2,sita的单位间隔角度为2°
maxDis = math.sqrt(self.ImageHeight * self.ImageHeight + self.ImageWidth * self.ImageWidth)
# DIS范围为±r(向上取整)
bDIS = int(np.ceil(maxDis / 2))
# 在中间补零位
DisNum = 2 * bDIS + 1
# 为了好计算角度搜索表还是定angle范围为[0,pi],num = 180 / 2 + 1
angleNum = 91
# 角度和弧度转换
zh = np.pi / 180
box = np.zeros((DisNum, angleNum), dtype=np.uint8)
SCsearch = np.zeros((2, angleNum), dtype=np.float64)
for j in range(angleNum):
SCsearch[0][j] = np.sin(j * 2 * zh)
for j in range(angleNum):
SCsearch[1][j] = np.cos(j * 2 * zh)
# 开始进行循环sita计算R进行投票
for i in range(self.ImageHeight):
for j in range(self.ImageWidth):
if edge[i][j] == 255:
for k in range(angleNum):
r = SCsearch[0][k] * j + SCsearch[1][k] * i
if r > 0:
box[bDIS + int(np.round(abs(r) / 2))][k] += 1
else:
box[bDIS - int(np.round(abs(r) / 2))][k] += 1
# 列非极大值抑制(c3)
for j in range(angleNum):
for i in range(0, DisNum-2, 3):
maxC = np.max(box[i:i+3, j:j+1])
for k in range(3):
if not box[i+k][j] == maxC:
box[i+k][j] = 0
# 确定投票认定直线阈值
MAX = np.max(box)
# 多次尝试的结果
thd = int(0.7 * MAX)
result = np.where(box > thd, 1, 0)
# 反算取出符合阈值的直线sita角和r
PiPei = []
for i in range(DisNum):
for j in range(angleNum):
if result[i][j] == 1:
sn = SCsearch[0][j]
cs = SCsearch[1][j]
r = 0
if i > bDIS:
r = (i-bDIS)*2
else:
r = (bDIS-i)*2
PiPei.append((sn, cs, r))
# 线条绘制
rr = np.copy(gray)
gg = np.copy(gray)
bb = np.copy(gray)
# 开始在合成影像上进行匹配点绘制
# 设定判定误差为2
error = 1
for t in PiPei:
a, b, r = t
for q in range(self.ImageHeight):
for w in range(self.ImageWidth):
if abs(r - (a*w+b*q)) < error:
rr[q][w] = 255
gg[q][w] = 0
bb[q][w] = 0
m = []
m.append(rr)
m.append(gg)
m.append(bb)
inputI = cv2.merge(m)
self.graphicRightshow(inputI)
四、结果展示
图一 阈值选取为0.6MAX
图二 阈值选取为0.7MAX
五、结果分析
经过不断调试直线判定阈值,其实不难发现我的代码存在角度分区过细后像素偏分的误差,同时,该误差在最后像素遍历又二次传播,再者Python自带的正弦和余弦计算式自然也携带一定计算误差,因此效果提取较为不好,但是由于初学的心态想记录一下自己的思考历程而作此文章。
我的算法在时间计算效率上其实有较大的改进空间,例如识别直线标记部分。还有对于像素偏分误差还可以从非极大值抑制部分代码进行改进,本人后续也会尝试修正,如果各位机智的网友们有好的建议,欢迎不吝赐教。
此次算法学习,我逐渐加深了对计算数学与计算机实现的Gap理解,同时也对自己的数学思维有了新的启发。最后,感谢各位用心传播知识的博主与伟大数学家和计算机领域的巨人们。
六、彩蛋
祝中国地质大学(武汉)校庆快乐!南望不忘,未来已来!
参考文献:
[1] 数字图像处理:第四版/(美)拉斐尔.C.冈萨雷斯(Rafael C.Gonzalez),(美)理查德.E.伍兹(Richard E.Woods)著;阮秋琦等译. —北京:电子工业出版社,2020.5.
[2]Duda,R.O.;Hart,P.E.(January 1972).Use of the Hough Transformation to Detect LInes and Curves in Pictures.Comm.ACM.15:11-15.doi:10.1145/361237.361242.
[3]Hough,P.V.C.Method and means for recognizing complex patterns.U.S.Patent 3069,654,1962.
[4]https://en.wikipedia.org/wiki/Hough_transform(Wiki)
[5]https://en.wikipedia.org/wiki/Polar_coordinate_system(Wiki)
[6]3.hough transform(霍夫变换直线检测)_哔哩哔哩_bilibili (引自博主:会飞的吴克)