CN110264508B - 一种基于凸四边形原则的消影点估计方法 - Google Patents

Abstract

本发明涉及一种基于凸四边形原则的消影点估计方法，属于视觉传感器技术领域；该方法通过实时对输入图像W进行图像处理获得图像内四方形物体的轮廓图像Q，然后对Q基于消影点定义提取图像中的线段，根据正交的平行直线对在图像平面内成像形成凸四边形与否来计算直线交点，并根据阈值过滤器选择聚类簇数，采用自适应复合聚类的方法来对直线交点进行聚类估计，得到消影点的集合。对比现有技术，本发明减少了估计运算量，保持了估计精度，同时节约了估计时间，不仅适合于任何拍摄、摄像设备的消隐点估计，而且适合应用于需要动态实时估计摄像设备姿态角度的场景。

Description

一种基于凸四边形原则的消影点估计方法

技术领域

本发明涉及一种基于凸四边形原则的消影点估计方法，属于视觉传感器技术领域。

背景技术

欧氏空间中两条平行直线相交于无穷远点,理想的小孔模型下,这两条平行直线经过摄像机得到的投影也为直线,且一般情况下其投影线相交,交点称为消隐点(vanishpoint),它是无穷远点在像平面上的投影。正交的两组平行线可以形成的一组与摄像头关心成直角的消影点，消隐包含了所在直线的方向信息，对其进行分析，可以得到三维空间的结构和方向信息，在三维重建、视觉陀螺仪、机器人视觉以及摄像机标定等领域被广泛运用。

在视觉陀螺仪、机器人视觉领域的实际应用中，面对相机连续姿态变化、以及视角变化、和距离变化导致相机姿态角度估计困难的问题，由于消影点具有稳定、连续、不随特征距离变化的特点，常常采用消影点来估计姿态角，消影点的位置随着姿态的变化而改变。

动态消影点估计技术目前主要依靠的是提取直线特征后采用K均值聚类或者大量直线相交求统计平均方法来获取消影点的位置。这样的统计方法是不稳健的，而且经典Hough变换计算量大，一些噪声比如椒盐噪声没有得到很好的处理。通过简单的k均值聚类并没有解决由于直线交点所引入的假消影点所带来的影响。这些方法也没有考虑到消隐点本身的误差和空间分布对估计精度影响。

针对如何区别由空间结构线相交得到的真消隐点和不是主要结构线与之相交得到的假消隐点的问题，本发明提出一种凸四边形的消影点的估计方法，考虑消影点形成原理和空间分布上进行位置的计算，并通过复合聚类的方法进行均值，减少了噪声或者成像畸变噪声的误差，减少了由于非正态分布的离群点所带来的误差，算法鲁棒性较好，精确度高。

针对基于霍夫空间的消影点检测方法，计算量大而且用于交点个数滤波筛选的过滤器设计困难、噪声处理不彻底的问题，本发明利用建筑物、道路、桥梁、屋顶、街道、墙等人造物体的强边缘信息，提出动态估计消影点位置的一种新方法，采用随机概率采样的方法进行特征提取，大大降低了算法运算量和复杂度，采用针对高斯、椒盐、泊松等噪声的滤波算法，减少了噪声误差。

发明内容

本发明的目的是为解决现有消影点动态估计方法中存在计算量大、过滤器设计困难、噪声处理不彻底、非正态离群点问题，旨在伴随着视角、距离的变化的相机连续姿态变化过程中，动态消影点位置的确定，提出一种基于凸四边形原则的消影点估计方法。

本发明的原理是对具有正交平行直线对的街景、房屋等场景下所拍摄的图像进行消影点提取，从图像边缘信息中，使用局部概率抽样的方法进行霍夫投影，设计了新的过滤器函数来提取图像中的线段，根据正交的平行直线对在图像平面内成像，形成凸四边形与否来计算直线交点，并根据阈值过滤器选择聚类簇数，采用自适应复合聚类的方法来对直线交点进行聚类估计，得到消影点的集合，该方法减少了估计运算量，保持了估计精度，同时节约了时间。

下文中遍历指的是对集合中的每个元素依次进行操作；

下文中出现的图像理解为灰度图像中所有像素点组成的图像矩阵；

图像矩阵中坐标轴原点在图像矩阵的左上角，坐标值都为正值；

超出图像矩阵边界的像素点的灰度值取零；

假设：

W为输入灰度图像中像素灰度值组成的矩阵，图像大小为m×n；

D为差分模板图像中像素值组成的矩阵，大小为m×n；

f为差分模板除燥后的图像中像素值组成的矩阵，大小为m×n；

K为行扫描自适应中值滤波后的图像中像素值组成的矩阵，大小为m×n；

K‘为梯度赋值极大值抑制后的图像中像素值组成的矩阵，大小为m×n；

Q为轮廓图像中像素值组成的矩阵，大小为m×n；

一种基于凸四边形原则的消影点估计方法，包括如下步骤：

步骤1对输入图像W进行图像处理得到只包含轮廓特征的图像Q；

作为优选，所述图像处理包括以下步骤：

步骤1.1对图像W采用多种差分尺度值提取其中的高斯差分图并融合形成差分模板图像D，使用D对W去噪得到去除高斯、泊松、秉性噪声的图像f；

作为优选，本步骤通过以下过程实现：

步骤1.1.1初始化(x,y)＝(1,1),m＝m₀,n＝n₀，

其中，(x,y)为像素坐标m₀,n₀为图像长和宽的初始值，

为差分尺度值；

步骤1.1.2依次从

中选取差分尺度值(σ₁,σ₂)，通过公式(1)计算高斯差分图像上点(x,y)的灰度值；

其中，(σ₁,σ₂)为差分尺度值,D_i(x,y,σ₁,σ₂)为第i对差分尺度值的差分图像上点(x,y)的灰度值，

为i取1,2,3时对应的差分尺度值，满足

和

W(x,y)为输入灰度图像W上点(x,y)的灰度值；

步骤1.1.3通过公式(2)计算D(x,y)；

D(x,y)＝max(|D₁(x,y,σ₁,σ₂)|，|D₂(x,y,σ₁,σ₂)|，|D₃(x,y,σ₁,σ₂)|) (2)

其中D(x,y)为差分模板D上点(x,y)处的灰度值，max为最大值函数；

步骤1.1.4通过公式(3)，计算f(x,y)：

其中f(x,y)是差分模板除燥后的图像f中点(x,y)处的灰度值；

步骤1.1.5如果x≤n且y≤m，则

跳转到步骤1.1.2否则输出差分模板除噪后的图像f和差分模板图像D；

其中，≤表示小于等于关系，&表示与关系，≠表示不等关系，＝＝表示相等关系。

步骤1.2对f采用行扫描的自适应中值滤波滤除椒盐噪声得到图像K；

作为优选，本步骤通过以下过程实现：

步骤1.2.1初始化(x,y)＝(1,1)，e＝e₀，k＝k₀其中(x,y)为像素坐标，e,k分别为动态窗口的大小和差值比例因子，e₀,k₀分别为其初值；

步骤1.2.2如果f(x,y)等于255或者0，则进入步骤1.2.3；否则进入步骤1.2.7；

步骤1.2.3以(x,y)位置的像素点为中心，创建滚动窗口，窗口的大小为e×e的矩形；

步骤1.2.4分别遍历窗口中每一行像素点的灰度值，得到极值与中位值序列，记录到集合Z_lmax＝{z_lmax}，Z_lmin＝{z_lmin}，Z_lmid＝{z_lmid}中；

其中，{}表示集合；

通过公式(4)计算动态窗口全局极值与中位值z_max，z_lmin，z_mid。

其中，max、min为最大值最小值函数；

步骤1.2.5如果z_max≥z_mid≥z_min，则窗口范围固定，进行窗口内像素点的滤波，进入步骤1.2.6；否则e＝e+2，并返回步骤1.2.3；

步骤1.2.6通过公式(5)计算K(x,y)，并更新像素矩阵K中点(x,y)处的灰度值；

其中，f_replace(x,y)＝z_mid，ν＝k(z_max-z_min)，Δ＝|f(x,y)-z_mid|，k为差值比例因子，K(x,y)图像K中点(x,y)的灰度值，≠描述不等关系，＝＝描述相等关系；

步骤1.2.7如果x≤n且y≤m，则

跳转到步骤1.2.2；否则输出的行扫描自适应中值滤波后的图像K；

其中，≤表示小于等于关系，&表示与关系，≠表示不等关系，＝＝表示相等关系。

步骤1.3对K计算梯度幅值矩阵H和梯度方向矩阵

并进行弧度极大值抑制后得到图像K'；

作为优选，本步骤通过以下过程实现：

步骤1.3.1初始化(x,y)＝(1,1)，其中(x,y)为像素坐标；

步骤1.3.2通过公式(6)，计算(x,y)处像素梯度H(x,y),以及估计梯度方向

其中K(x,y)为图像K中点(x,y)的灰度值，atan为反正切函数；

步骤1.3.3如果x≤n且y≤m，则

并跳转到步骤1.3.2；否则，返回梯度幅值矩阵H和梯度方向矩阵

并进入步骤1.3.4；

其中，≤表示小于等于关系，&表示与关系，≠表示不等关系，＝＝表示相等关系；

步骤1.3.4初始化(x,y)＝(1,1)，其中(x,y)为像素坐标；

步骤1.3.5根据

的值，通过公式(7)计算像素点(x,y)处极大值抑制阈值H'(x,y)；

通过公式(8)，使用H'(x,y)对K(x,y)进行极大值抑制，计算K'(x,y)；

其中，K'(x,y)为对图像K进行极大值抑制后的图像中点(x,y)处的灰度值；

步骤1.3.6如果x≤n&y≤m，则

跳转到步骤1.3.5；否则输出的梯度极值筛选后的图像K'；

其中，≤表示小于等于关系，&表示与关系，≠表示不等关系，＝＝表示相等关系。

步骤1.4对K'进行高阈值滤波得到图像K^*，再结合差分模板图像D进行互补连接得到轮廓图像Q。

作为优选，本步骤通过以下过程实现：

步骤1.4.1初始化(x,y)＝(1,1)，TH＝TH₀；其中(x,y)为像素坐标,TH是高阈值滤波的阈值，TH₀是其初始化的值；

步骤1.4.2通过公式(9)，计算K^*(x,y)；

其中，K^*(x,y)为图像K^*中(x,y)点处的灰度值；

步骤1.4.3如果x≤n&y≤m，则

并跳转到步骤1.4.2,；否则，输出图像K^*并进入步骤1.4..4；

其中，≤表示小于等于关系，&表示与关系，≠表示不等关系，＝＝表示相等关系；

步骤1.4.4初始化(x,y)＝(1,1)，i*＝i₀，初始化图像Q中所有像素点的值为0；

其中，(x,y)为像素坐标，i*为正方形搜索框的边长，i₀为其初值；

步骤1.4.5如果K^*(x,y)不为0，则Q(x,y)＝K^*(x,y)并进入步骤1.4.6，否则，跳入步骤1.4.9；

步骤1.4.6初始化(i,j)＝(0,0),其中i,j矩阵搜索框的偏移量；

步骤1.4.7如果K^*(x+i,y+j)不为0，则Q(x+i,y+j)＝K^*(x+i,y+j)，并进入步骤1.4.9；否则，进入步骤1.4.8；

步骤1.4.8如果D(x+i,y+j)不为0，则Q(x+i,y+j)＝D(x+i,y+j)，并进入步骤1.4.9；

步骤1.4.9如果i≤i*&y≤i*则

并返回步骤1.4.7，否则进入步骤1.4.10；

步骤1.4.10如果x≤n&y≤m，则

并跳转到步骤1.4.5；否则，结束循环，输出轮廓图像Q。

步骤2输入边缘轮廓图像Q，通过随机局部霍夫空间中的阈值滤波方法，提取线段，并进行优化更新线段所在直线的参数，得到线段集合Ln，步骤完毕后进入步骤3；

作为优选，本步骤通过以下过程实现：

步骤2.1初始化q＝q₀,g＝g₀,R＝R₀,w*＝w₀,m*＝m₀,δ＝δ₀，Le＝Le₀并在图像Q中得到采样点i(a,b)，在每间隔q个像素点的区间中选择像素点灰度值最大的点作为采样点，并添加到采样点集合I中；

其中，q为采样间隔，g为细分数，g₀为其初值，R为兴趣半径，R₀为其初值，w*,m*为合并精度，w₀,m₀为其初值，δ为搜索精度参数，δ₀为其初值，Le为线段长度阈值，Le₀为其初值，i(a,b)∈I为图像Q中坐标为(a,b)的灰度值，I为采样集合；

步骤2.2随机在采样集合I抽样选择一点i(x,y)，在图像Q以(x,y)位置为中心，以区域半径R＝R₀的圆形为兴趣区域进行遍历，将该区域中的灰度值不等于零的像素点，灰度值置为最大，并添加到兴趣集合的点集Υ中；

步骤2.3通过公式(10)，将兴趣集合Υ中点映射为正弦曲线，计算正弦曲线在霍夫空间中的交点坐标，并添加到集合N^#中；

其中，(a,b)为集合Υ中点的坐标，(β,η)为霍夫空间中的正弦曲线的坐标值；

步骤2.4遍历集合N^#中的元素(ρ₀,θ₀)，初始化(ρ,θ)＝(ρ₀,θ₀)和k＝1，接着全部遍历集合N^#-{(ρ₀,θ₀)}中的元素，通过公式(11)进行依次迭代，计算(ρ,θ,k)，全部迭代完毕后将(ρ,θ,k)添加到重构集合N*中,遍历完成后删除重构集合N*中重复的元素；

其中，(ρ,θ,k)为合并后元素的值，(ρ,θ)为均值坐标，k为迭代过程中均值坐标的个数，(ρ₀,θ₀)为集合N^#中的元素，(ρ₁,θ₁)为集合N^#中除(ρ₀,θ₀)之外的元素，≤表示小于等于关系，&表示与关系，≠表示不等关系，＝＝表示相等关系，≥表示大于等于关系，|表示或的关系，{}表示集合；

步骤2.5在极坐标空间中对横纵坐标的区间进行细分，得到g组细分区间为

[0,Δθ],[Δθ,2Δθ],[2Δθ,3Δθ]…[gΔθ,π],[0,Δρ],[Δρ,2Δρ],[2Δρ,3Δρ]…[gΔρ,ρ_max]；

其中g为细分数，ρ_max为集合N‘参数中的ρ的最大值，Δρ，Δθ分别为细分区间的长度；

步骤2.6初始化p₁＝0，p₂＝0,T＝T₀，其中，T为投票阈值，T₀为其初始值,p₁,p₂为细分区间的序号；

步骤2.7遍历N*中的元素,初始化n＝0，通过式子(12)，计算E(ρ',θ')并添加到集合E中,同时通过式子(13),据细分区间划分N*中的元素，计算V(ρ',θ')并添加到集合V中；

其中ρ'＝[p₁Δρ,(p₁+1)Δρ],θ'＝[p₂Δθ,(p₂+1)Δθ]，为细分区间，p₁,p₂为细分区间的序号，n为迭代过程中均值坐标在细分区间[p₁Δρ,(p₁+1)Δρ],[p₂Δθ,(p₂+1)Δθ]内的个数，E(ρ',θ')为集合N*中均值坐标在细分区间[p₁Δρ,(p₁+1)Δρ],[p₂Δθ,(p₂+1)Δθ]内的个数，V(ρ',θ')为集合N*中均值坐标在细分区间[p₁Δρ,(p₁+1)Δρ],[p₂Δθ,(p₂+1)Δθ]内的均值坐标集合，∪为集合间的并集运算，∈表示属于符号，

为不属于符号，{}表示集合；

步骤2.8如果E(ρ',θ')＞T，则返回(ρ',θ')并进入步骤2.10,否则进入步骤2.9，其中，T为投票阈值；

步骤2.9如果p₁＜g且p₂＜g则

并跳转到步骤2.7，否则跳入步骤2.14；

步骤2.10将(ρ',θ')代入V中得到点集V(ρ',θ')，通过式子(14)计算搜索方向的直线参数ρ*,θ*；

其中(ρ_m,θ_m)是V(ρ',θ')中的j个点中的第m个点的坐标值，(ρ*,θ*)为极坐标描述直线的解析参数，即所描述方向的法线段的长度和其所在直线倾角；

步骤2.11通过式子(15)计算Ln(ρ*,θ*)，在图像Q中沿着直线，搜索像素灰度值不为0的点的坐标，并添加到集合Ln(ρ*,θ*)中；

其中，(a,b)是搜索直线方向上点的坐标，Ln(ρ*,θ*)为图像Q中在解析参数为ρ*,θ*的直线上的像素点位置的集合，acosθ*+bsinθ*＝ρ*为搜索直线的方程；

步骤2.12通过公式(16)计算更新参数ρ^,θ^；

其中(ρ^{^},θ^{^})为拟合后直线的法线段的长度以及其所在直线倾角，为极坐标描述直线的解析参数；其中，(a_i,b_i)是Ln(ρ*,θ*)中n个元素中的第i个元素，a₀,b₀为笛卡尔坐标系下的直线的斜率和截距，atan为反正切函数；

步骤2.13通过公式(17)更新参数ρ*,θ*，并通过式子(18)将Ln(ρ*,θ*)添加到到线段集合Ln中，并返回步骤2.9；

其中，(a_min,b_min),(a_max,b_max)是Ln(ρ*,θ*)中的横坐标最小和最大对应的两个的坐标值，Le为线段长度阈值；

步骤2.14判断集合I是否为

若为空则结束算法，返回线段集合Ln，若不为空则跳入步骤2.2，继续循环。

步骤3输入线段集合Ln，根据倾角阈值函数，划分直线簇，并基于三角形内一点的空间位置的判断方法，得到组成凸四边形的所有直线簇集合，计算每组直线簇中对边的交点，得到消影点位置的集合M，步骤完成后进入步骤4；

作为优选，本步骤通过以下过程实现：

步骤3.1初始化E₁＝E₁₀，E₂＝E₂₀，T₁＝T₁₀,T₂＝T₂₀，遍历线段集合Ln，通过公式(19)，找到所有满足条件的线段组，并统计所有线段组所在直线的参数，建立直线簇集合L；

L₁,L₂,L₃,L₄满足

L₁,L₂,L₃,L₄∈Ln

其中，E₁，E₂为直线倾角差值的阈值，E₁₀E₂₀为初值，T₁,T₂为聚类族数精度调整参数，T₁₀,T₂₀为其初值，L₁，L₂，L₃，L₄是线段集合Ln中满足要求的一组线段,

为L₁,L₂,L₃,L₄在Ln中对应的索引值，即该点集所在直线的解析参数，

为点集所在直线倾角的差值，∪表示取集合间的并运算，L为所有满足阈值条件的直线簇的解析参数的集合，{}表示集合，&表示与关系；

步骤3.2初始化r＝1，并计算Len＝length(L)；

其中，r为直线簇集合的索引，Len为集合L中直线簇的个数，length为集合长度计算函数；

步骤3.3将公式(20)代入公式(21)中，取

计算(x_A,y_A),(x_B,y_B),(x_C,y_C),(x_D,y_D)；

其中，mL₁,mL₂,mL₃,mL₄为直线簇L(r)中对应直线的标号，X(L*,C*),Y(L*,C*)为横纵坐标的计算函数，L*,C*为其输入参数，ρ_L*,θ_C*为直线参数，A,B,C,D为直线簇按照序号顺序依次相交的交点，(x_A,y_A),(x_B,y_B),(x_C,y_C),(x_D,y_D)为其坐标，L(r)为L中第r组直线簇参数，

为第r组直线簇中对应标号为mL₁,mL₂,mL₃,mL₄的直线的参数；

步骤3.4建立方程组(22)；

其中，u和v是方程组的解，

是二维平面中S₁,S₂.S₃.S₄四个点的坐标；

根据式子(23)依次取值并代入方程组(22)中，

依次求解方程组；若方程组的解都满足u＜0|v＜0|v+u＞1，则表示A,B,C,D四个点组成了凸四边形，跳入步骤3.6，否则进入步骤3.5；

步骤3.5如果r＜Len，则r＝r+1进入步骤3.3，否则进入步骤3.7；

步骤3.6通过公式(24)，将组成凸四边形的直线簇中各个直线的解析参数添加到凸四边形直线簇集合Q_LG中；

其中，{}表示集合，∪表示集合间的并运算；

步骤3.7计算Lo＝length(Q_LG)并初始化N＝2×Lo，N₁＝0；

其中N为聚类族群数，Lo为凸四边形直线簇集合Q_LG的长度，length为集合长度计算函数，N₁为聚类族群数的迭代参数；

步骤3.8遍历集合Q_LG中的元素得

接着全部遍历集合

中的元素，并通过公式(25)迭代更新N₁，所有遍历结束后计算

并进一步判断更新聚类族群数为

其中，

为集合Q_LG中的元素，{}表示集合；

为集合Q_LG中除了元素

之外的元素；

步骤3.9遍历Q_LG中的元素得

通过公式(26)粗略计算消影点的坐标(u₁,v₁),(u₂,v₂)，并添加到消影点集合M中；

其中，(u₁,v₁),(u₂,v₂)是凸四边形直线簇中对边直线得到的交点坐标，对边指的是

所描述的直线对和

所描述的直线对；

步骤4输入消影点位置集合M和复合聚类族群数N，通过自适应复合聚类，得到消影点集合M的N个聚类中心的位置坐标

作为优选，本步骤通过以下过程实现：

步骤4.1随机从消影点集M中的选择N个点的坐标作为聚类族群中心点的坐标初值

开始聚类，并将循环次数初始化为k＝1，初始化γ＝γ₀，α＝α₀,β＝β₀,TG＝TG₀，γ为自适应率，γ₀为其初值,α,β为互补因子，α₀,β₀为其初值，TG为精度阈值，TG₀为其初值；

步骤4.2初始化G＝1；

步骤4.3通过公式(27)，计算M(G)到各个群类中心点的复合距离

其中，M(G)为M中第G个消影点；(X,Y)为M(G)的坐标，

为第k次迭代后第i类族群的中心点的坐标，i取值为1,2,3,...N，α+β＝1为互补因子，

为第k次迭代后点(X,Y)到第i类中心点的欧式距离，

为第k次迭代后(X,Y)到第i类中心点的余弦距离，

为第k次迭代后点(X,Y)到第i类中心点的复合距离；

计算过程中，通过公式(28)，调整互补因子α,β且0≤α,β≤1；

其中,γ为自适应率；

为边缘阈值，d＝min(m,n)，图像大小为m×n。

步骤4.4将距离

按大小排序，将M(G)归类到最近的族群类中；

步骤4.5判断如果G≥length(M)，进入步骤4.6，否则G＝G+1返回步骤4.3；其中，length为集合长度函数计算函数；

步骤4.6通过公式(29)重新计算族群的中心点；

其中X_kj,Y_kj，表示第k次迭代后i类中的第j个数据点，M_ki表示第k次迭代后i类中点的个数，class ki表示第k次迭代后i类的点的集合；

步骤4.7遍历第k次循环后的聚类中心

和第k-1循环后的聚类中心

TG为精度阈值，如果存在

则k＝k+1，并跳转到步骤4.2，否则退出聚类循环，完成运算，进入步骤4.6；其中，<>表示2-范数；

步骤4.8输出当前

的值，作为消影点集合M聚类后N个族群中心的位置。

有益效果

一种基于凸四边形原则的消影点估计方法，相比于现有技术，本方法有如下有益效果：

1.本发明所述方法应用于消影点的估计过程，发现在估计过程中，由于根据解析几何直接计算消影点，并通过聚类来二次估计消影点位置，降低了计算量，而且运算复杂度低，精度得到提高。

2.本发明所述方法能够应用于更多场景比如说街景、房屋等场景下的消影点估计，能够自动识别筛选出图像中的边缘信息，提取出消影点，适用性更强。

3.本发明采用了高斯差分模板来进行对高斯、泊松、秉性等噪声的去噪，相比较LOG算子滤波，在保证去噪效果的前提下，降低了运算量。

4.本发明针对椒盐噪声，设计了行扫描的自适应中值滤波流程，提高了运算速度，较为全面的去除了噪声的影响，去噪效果好。

5.本发明采用了弧度极大值抑制的方法来增强图像的边缘特征，考虑了多方向的因素，增强了抑制结果的鲁棒性和稳定性。

6.本发明采用了基于差分模板图像和高阈值滤波相结合的连接方法，引入了高斯平滑效果，提高了轮廓图像的清晰度。

7.本发明采用局部概率霍夫空间映射来提取线状特征后使用最小二乘优化特征参数，降低了特征提取的计算量并提高了特征筛选的精度。

8.本发明提出了凸四边形方法，根据消影点形成原理和空间分布来设计消影点的筛选的阈值函数，解决了非正态分布的离群点所带来的误差，提高了估计精度。

9.本发明采用自适应复合聚类的方法来对消影点进行聚类估计，可以减少由于噪声或者成像畸变而造成的误差。

附图说明

图1为一种基于凸四边形原则的消影点估计方法的整体流程图。

图2为图1中线状特征提取方法的详细流程图。

图3为图1中凸四边形阈值过滤方法的详细流程图。

图4为图1中自适应复合聚类方法的详细流程图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

在人工环境中如道路、房屋等存在许多空间上的正交平行直线对经过透视投影成像后，空间平行的直线将在成像平面中分别相交形成消影点，而这些消隐点具备物理性质，对于相机或摄像头姿态的估计有着重要的意义。实际应用中，面对相机连续姿态变化、以及视角变化、和距离变化导致相机姿态角度估计困难的问题，由于消影点具有稳定、连续、不随特征距离变化的特点，常常采用消影点来估计姿态角，消影点的位置随着姿态的变化而改变。

本发明提供了一种基于凸四边形原则的消影点估计方法，具体实施包括如下步骤，但本发明的具体实施形式并不局限于此。

数据拍摄于北京理工大学国防科技园XXX实验室，场景包含实验柜、实验台、计算机、桌子、椅子、墙、屋檐等图像边缘较强的人造对象，连续改变相机视角和姿态，拍摄多张图片，图片规格为1200X1500，jpeg格式，将序列图像按照时间顺序，每次输入一张，进行每张图像消影点位置的估计，而且在自动驾驶的行进器等要求消隐点检测算法快速而对精度要求不高、实时性强的条件下，也有很好的应用价值。具体实施步骤如下：

实施例1

本实施例阐述了将本发明“一种基于凸四边形原则的消影点估计方法”应用于从包含实验柜、实验台、计算机、桌子、椅子、墙、屋檐等边缘较强的人造对象的场景所拍摄的图像中降低运算复杂度，高精度估计消影点位置的算法流程，如图1所示。

步骤1输入图像W并选择多种差分尺度值提取原图像的高斯差分图融合形成差分模板用于高斯、泊松、秉性噪声的去噪，输出差分模板除噪后的图像f和差分模板图像D，去除噪声步骤完毕后跳至步骤2。

优选的，本步骤通过以下过程实现：

步骤1.1具体实施的话，首先根据输入图像的分辨率要求，对图像大小进行初始化设m＝1200,n＝1500，W(x,y)为图像W上点(x,y)的灰度值。以及进行图像遍历的初始坐标为(x,y)＝(1,1)，设置差分尺度值

步骤1.2依次选择σ₁，σ₂值，并根据公式(30)<(31),(32)分步卷积依次计算3种差分图像上点(x,y)的灰度值D₁，D₂，D₃；

步骤1.3通过公式(33)计算差分模板图像D上点(x,y)的灰度值D(x,y)；

D(x,y)＝max(|D₁(x,y,0.4,0.3)|，|D₂(x,y,0.6,0.7)|，|D₃(x,y,0.7,0.8)|) (33)

步骤1.4通过公式(3)，计算f(x,y)。其中f(x,y)是f中点(x,y)处的灰度值；

步骤1.5如果x≤1500&y≤1200，则

跳转到步骤1.2否则输出差分模板除燥后的图像f和差分模板图像D；

步骤2输入图像f上，针对椒盐噪声，进行行扫描的自适应中值滤波，输出的行扫描自适应中值滤波后的图像K，步骤完毕后跳转至步骤3。

优选的，本步骤通过以下过程实现：

步骤2.1初始化(x,y)＝(1,1)，e＝5，k＝0.6；

步骤2.2如果f(x,y)等于255或者0，则进入步骤2.3；否则进入步骤2.7；

步骤2.3以(x,y)位置的像素点为中心，创建滚动窗口，窗口的大小为5×5的矩形；

步骤2.4分别遍历窗口中每一行像素点灰度值，得到极值与中位值序列，记录到集合Z_lmaxZ_lmax，Z_lmin，Z_lmid中。通过公式(4)计算动态窗口全局极值与中位值z_max，z_lmin，z_mid；步骤2.5如果z_max≥z_mid≥z_min，则窗口范围固定，进行窗口内像素点的滤波，进入步骤2.6；否则e＝e+2，并返回步骤2.3；

步骤2.6通过公式(34)计算K(x,y)，并更新像素矩阵K中点(x,y)处的灰度值；

其中，f_replace(x,y)＝z_mid，Δ＝|f(x,y)-z_mid|，K(x,y)为图像K中点(x,y)的灰度值。

步骤2.7如果x≤1500&y≤1200，则

并跳转到步骤2.2；否则输出的中值滤波后的图像K；

步骤3输入图像K，计算梯度幅值矩阵H和梯度方向矩阵

并进行弧度极大值抑制，输出K经过极大值抑制后的图像K'步骤完毕后进入步骤4。

优选的，本步骤通过以下过程实现：

步骤3.1初始化(x,y)＝(1,1)；

步骤3.2通过公式(6)，计算(x,y)处像素梯度H(x,y),以及估计梯度方向

步骤3.3如果x≤1500&y≤1200，则

并跳转到步骤3.2；否则，返回梯度幅值矩阵H和梯度方向矩阵

并进入步骤3.4。

步骤3.4初始化(x,y)＝(1,1)，其中(x,y)为像素坐标；

步骤3.5沿着

方向，通过公式(7)计算像素点(x,y)处极大值抑制阈值H'(x,y),通过公式(8)，使用H'(x,y)对K(x,y)进行极大值抑制计算K'(x,y)。其中，K'(x,y)为图像K'中(x,y)点的灰度值；

步骤3.6如果x≤1500&y≤1200，则

跳转到步骤3.5；否则输出的梯度极值筛选后的图像K'；

步骤4对图像K'进行高阈值滤波得到图像K^*，再结合差分模板图像D进行互补连接得到轮廓图像Q，Q(x,y)为轮廓图像中(x,y)点处的灰度值。步骤完成进入步骤5。

优选的，本步骤通过以下过程实现：

步骤4.1初始化(x,y)＝(1,1)，TH＝10。其中(x,y)为像素坐标,TH是高阈值滤波的阈值；

步骤4.2通过公式(9)，计算K^*(x,y),其中，图像K^*中(x,y)点处的灰度值；

步骤4.3如果x≤1500&y≤1200，则

并跳转到步骤4.2,；否则，输出图像K^*并进入步骤4.4；

步骤4.4初始化(x,y)＝(1,1)，i*＝2,图像Q中所有像素点的值为0；步骤4.5如果K^*(x,y)≠0则Q(x,y)＝K^*(x,y),进入步骤4.6,否则跳入步骤4.9。

步骤4.6初始化(i,j)＝(0,0)。

步骤4.7如果K^*(x+i,y+j)≠0，则Q(x+i,y+j)＝K^*(x+i,y+j)，并进入步骤4.9；否则进入步骤4.8，

步骤4.8如果D(x+i,y+j)≠0则Q(x+i,y+j)＝D(x+i,y+j)。并进入步骤4.9

步骤4.9如果i≤2&y≤2则

并返回步骤4.7，else进入步骤4.10

步骤4.10如果x≤1500&y≤1200，则

并跳转到步骤4.5；否则结束循环，输出轮廓图像Q；

本领域技术人员清楚，上述对输入图像进行图像处理以获得图像中物体轮廓的过程不限于所述方法，可以采用现有的任何图像处理方法实现，发明人采用的上述方法是能够精确提取轮廓的优选方法。

步骤5输入边缘轮廓图像Q，通过随机局部霍夫空间中的阈值滤波方法，提取线段，并进行优化更新线段所在直线的参数，得到线段集合Ln，步骤完毕后进入步骤6。

优选的，本步骤通过以下过程实现：

步骤5.1初始化q＝3，g＝100，R＝5，

δ＝0.5，Le＝4，并在图像Q中每间隔3个像素点的区间中选择像素点灰度值最大的点作为采样点，并添加到采样集合I中；

步骤5.2随机在采样集合I抽样选择一点i(x,y)，在图像Q以(x,y)位置为中心，以区域半径R＝3的圆形为兴趣区域进行遍历，将该区域中的灰度值不等于零的像素点，灰度值置为最大，并添加到兴趣集合的点集Υ中；

步骤5.3通过公式(10)，将兴趣集合Υ中点映射为正弦曲线,并计算映射曲线在霍夫空间中的交点坐标，添加到集合N^#中；

步骤5.4遍历集合N^#中的元素(ρ₀,θ₀)，初始化(ρ,θ)＝(ρ₀,θ₀)和k＝1，接着全部遍历集合N^#-{(ρ₀,θ₀)}中的元素，通过公式(11)进行依次迭代，计算(ρ,θ,k)，全部迭代完毕后将(ρ,θ,k)添加到重构集合N*中,遍历完成后删除重构集合N*中重复的元素；具体包括的实例操作子步骤如下：

步骤5.4.A初始化m₁＝1，进入步骤5.4.B；其中，m₁为遍历的索引值；

步骤5.4.B置k＝1，(ρ,θ)＝(ρ₀,θ₀)＝N^#(m₁)，遍历集合N^#-{(ρ₀,θ₀)}并通过式子(35)进行迭代计算(ρ,θ,k)，遍历完毕后停止迭代，并将(ρ,θ,k)添加到重构集合N*，进入步骤5.4.C；其中，N^#(m₁)为集合N^#中第m₁个元素

步骤5.4.C如果m₁＞length(N^#)，则进入步骤5.4.D，否则m₁＝m₁+1，并返回步骤5.4.B；其中，length为集合长度计算函数。

步骤5.4.D删除重构集合N*中重复的元素；

步骤5.5在极坐标空间中对横纵坐标的区间进行细分，得到g组细分区间为

[0,Δθ],[Δθ,2Δθ],[2Δθ,3Δθ]…[gΔθ,π],[0,Δρ],[Δρ,2Δρ],[2Δρ,3Δρ]…[gΔρ,ρ_max]；

步骤5.6初始化p₁＝0，p₂＝0,T＝10；其中，T为投票阈值,p₁,p₂为细分区间的序号。

步骤5.7遍历N*中的元素,初始化n＝0，通过式子(12)，计算E(ρ',θ')并添加到集合E中，同时通过式子(13),根据细分区间划分N*中的元素，计算V(ρ',θ')并添加到集合V中；

步骤5.8如果E(ρ',θ')＞10，则返回(ρ',θ')并进入步骤5.10否则进入步骤5.9，其中，T为投票阈值；

步骤5.9如果p₁＜100且p₂＜100则

并跳转到步骤5.7，否则跳入步骤5.14。

步骤5.10将(ρ',θ')代入V中得到点集V(ρ',θ')，通过式子(14)计算搜索方向的直线参数ρ*,θ*；

步骤5.11通过式子(15)计算Ln(ρ*,θ*)，在图像Q中沿着直线，搜索像素灰度值不为0的点的坐标，并添加到集合Ln(ρ*,θ*)中；

步骤5.12通过公式(16)计算更新参数ρ^{^},θ^{^}；

步骤5.13通过公式(17)更新参数ρ*,θ*，并通过式子(36)将Ln(ρ*,θ*)添加到到线段集合Ln中，并返回步骤5.9；

步骤5.14删除集合I中的点i(x,y)，判断集合I是否为

若为空则结束算法，返回线段集合Ln，若不为空则跳入步骤5.2，继续循环。

步骤6输入线段集合Ln，根据倾角阈值函数，划分直线簇，并基于三角形内一点的空间位置的判断方法，得到组成凸四边形的所有直线簇集合，计算每组直线簇中对边的交点，得到消影点位置的集合M，步骤完成后进入步骤7。

优选的，本步骤通过以下过程实现：

步骤6.1初始化

T₂＝5，E₁，E₂其中直线倾角差值的阈值，T₁,T₂为聚类族数精度调整参数，遍历线段集合Ln，通过公式(37)，找到所有满足条件的线段组，并统计所有线段组所在直线的参数，建立直线簇集合L；

L₁,L₂,L₃,L₄∈Ln

步骤6.2初始化r＝1，并计算Len＝length(L)；

步骤6.3取

将公式(20)的L*,C*依次取值，即将mL₁,mL₂、mL₂,mL₃、mL₃,mL₄、mL₄,mL₁依次代入(20)中的函数X(L*,C*)，Y(L*,C*)中，首先将(X(mL₁,mL₂),Y(mL₁,mL₂)))、(X(mL₂,mL₃),Y(mL₂,mL₃))、(X(mL₁,mL₂),Y(mL₁,mL₂)))、(X(mL₂,mL₃),Y(mL₂,mL₃))中的

通过公式(38)进行值的校正，然后再计算的函数值，得到直线簇交点A,B,C,D的坐标(x_A,y_A),(x_B,y_B),(x_C,y_C),(x_D,y_D)；

步骤6.4通过公式(22)建立方程组,根据式子(23)依次取值，并代入式子(22)中，

依次求解方程组，若方程组的解都满足u＜0|v＜0|v+u＞1，则表示A,B,C,D四个点组成了凸四边形，否则进入步骤6.5；

步骤6.5如果r＜Len，则r＝r+1进入步骤6.3，否则进入步骤6.6；

步骤6.6通过公式(24)，将组成凸四边形的直线簇中各个直线的解析参数添如到凸四边形直线簇集合Q_LG中；

步骤6.7计算Lo＝length(Q_LG)并初始化N＝2×Lo，N₁＝0；

步骤6.8遍历集合Q_LG中的元素得

接着全部遍历集合

中的元素，并全部通过公式(25)迭代更新N₁，遍历完成后更新聚类族群数为

具体包括的实例操作子步骤如下：

步骤6.8.A初始化m₂＝1，进入步骤6.8.B；其中，m₂为遍历的索引值；

步骤6.8.B置

遍历集合

并通过式子(25)迭代更新N₁，遍历完毕后进入步骤6.8.C；

步骤6.8.C如果m₂＞length(Q_LG)，则进入步骤6.8.D，否则m₁＝m₁+1，并返回步骤6.8.B；其中，length为集合长度计算函数；

步骤6.8.D更新聚类族群数为

步骤6.9遍历Q_LG中的元素得

通过公式(26)粗略计算消影点的坐标(u₁,v₁),(u₂,v₂)，并将添加到消影点集合M中；

步骤7输入消影点位置集合M和复合聚类族群数N，通过自适应复合聚类，得到消影点集合M的N个聚类中心的位置坐标

步骤完成后，至此所有步骤结束。

优选的，本步骤通过以下过程实现：

步骤7.1随机从消影点集M中的选择N个点的坐标作为聚类初值，开始聚类，循环次数初始化为k＝1，初始化γ＝0.01，α＝1,β＝0,TG＝0.5；其中，γ为自适应率，α,β为互补因子，TG为精度阈值。

步骤7.2初始化G＝1；

步骤7.3通过公式(27)，计算M(G)到各个群类中心点的复合距离

计算过程中，通过公式(39)，调整互补因子α,β其中0≤α,β≤1；

步骤7.4按照距离

的大小排序，将M(G)归类到最近的族群类中；

步骤7.5判断如果G≥length(M)，进入步骤7.6，否则G＝G+1返回步骤7.3；

步骤7.6通过公式(29)重新计算族群的中心点；

步骤7.7遍历第k次循环后的聚类中心

和第k-1循环后的聚类中心

如果存在

则k＝k+1，并跳转到步骤7.2，否则退出聚类循环，完成运算，进入步骤7.8；

步骤7.8输出当前

的值，作为消影点集合M聚类后N个族群中心的位置。

至此，从步骤1到步骤7完成了本实施例一种基于凸四边形原则的消影点估计方法。

通过考虑如何区别由空间结构线相交得到的真消隐点和不是主要结构线与之相交得到的假消隐点的问题，本发明提出一种凸四边形的消影点的估计方法，从消影点形成原理和空间分布上进行位置的计算，并通过复合聚类的方法进行均值，避免了统计估计方法所带来的巨大的运算量，而且减少了由于非正态分布的离群点所带来的误差，提高了鲁棒性和精确度。

通过考虑常规霍夫空间的消影点检测方法所带来的运算量大的问题，采用随机概率采样并进行局部特征提取的方法进行运算，加快了运算速度，增强了实时性，而且考虑了直线交点个数滤波筛选的过滤器设计问题，提出了阈值筛选的策略，精确了搜索范围，而且采用了归类化处理的思想，批量化处理，提高了运行效率。同时，考虑到了处理过程中椒盐噪声存在所带来的影响，提出了相应的滤波算法，减少了图像的噪声误差，提高了估计精度。

为了说明本发明的内容及实施方法，本说明书给出了上述具体实施例。但是，本领域技术人员应理解，本发明不局限于上述最佳实施方式，任何人在本发明的启示下都可得出其他各种形式的产品，但不论在其形状或结构上作任何变化，凡是具有与本申请相同或相近似的技术方案，均落在本发明的保护范围之内。

Claims (9)

1.一种基于凸四边形原则的消影点估计方法，其特征在于：包括以下步骤：

(1)对输入图像W进行图像处理得到只包含轮廓特征的图像Q；

(2)对Q通过随机局部霍夫空间中的阈值滤波方法提取线段得到线段集合Ln；

(3)对Ln根据倾角阈值函数划分直线簇，并基于三角形内一点的空间位置的判断方法，得到组成凸四边形的所有直线簇集合，计算每组直线簇中对边的交点，得到消影点位置的集合M；

(4)对M和预置的复合聚类族群数N，通过自适应复合聚类方法得到消影点集合M的N个聚类中心的位置坐标

2.根据权利要求1所述的方法，其特征在于：所述图像处理包括去噪、弧度极大值抑制和高阈值滤波，具体如下：

1.1对图像W采用多种差分尺度值提取其中的高斯差分图并融合形成差分模板图像D，使用D对W去噪得到去除高斯、泊松、秉性噪声的图像f；

1.2对f采用行扫描的自适应中值滤波滤除椒盐噪声得到图像K；

1.3对K计算梯度幅值矩阵H和梯度方向矩阵

并进行弧度极大值抑制后得到图像K′；

1.4对K'进行高阈值滤波得到图像K^*，再结合差分模板图像D进行互补连接得到轮廓图像Q。

3.根据权利要求2所述的方法，其特征在于：所述1.1通过以下过程实现：

步骤1.1.1初始化(x,y)＝(1,1),m＝m₀,n＝n₀，

其中，(x,y)为像素坐标，m₀,n₀为图像长和宽的初始值，

为差分尺度值；

步骤1.1.2依次从

中选取差分尺度值(σ₁,σ₂)，通过公式(1)计算高斯差分图像上点(x,y)的灰度值；

其中，(σ₁,σ₂)为差分尺度值,D_i(x,y,σ₁,σ₂)为第i对差分尺度值的差分图像上点(x,y)的灰度值，

为i取1,2,3时对应的差分尺度值，满足

和

W(x,y)为输入灰度图像W上点(x,y)的灰度值；

步骤1.1.3通过公式(2)计算D(x,y)；

D(x,y)＝max(|D₁(x,y,σ₁,σ₂)|，|D₂(x,y,σ₁,σ₂)|，|D₃(x,y,σ₁,σ₂)|) (2)

其中D(x,y)为差分模板D上点(x,y)处的灰度值，max为最大值函数；

步骤1.1.4通过公式(3)，计算f(x,y)：

其中f(x,y)是差分模板除燥后的图像f中点(x,y)处的灰度值；

步骤1.1.5如果x≤n且y≤m，则

跳转到步骤1.1.2否则输出差分模板除噪后的图像f和差分模板图像D；

其中，≤表示小于等于关系，≠表示不等关系，＝＝表示相等关系。

4.根据权利要求2所述的方法，其特征在于：所述1.2通过以下过程实现：

步骤1.2.1初始化(x,y)＝(1,1)，e＝e₀，k＝k₀其中(x,y)为像素坐标，e,k分别为动态窗口的大小和差值比例因子，e₀,k₀分别为其初值；

步骤1.2.2如果f(x,y)等于255或者0，则进入步骤1.2.3；否则进入步骤1.2.7；其中，f(x,y)表示图像f在像素坐标(x,y)处的灰度值；

步骤1.2.3以(x,y)位置的像素点为中心，创建滚动窗口，窗口的大小为e×e的矩形；

步骤1.2.4分别遍历窗口中每一行像素点的灰度值，得到极值与中位值序列，记录到集合Z_lmax＝{z_lmax}，Z_lmin＝{z_lmin}，Z_lmid＝{z_lmid}中；

其中，{}表示集合；

通过公式(4)计算动态窗口全局极值与中位值z_max，z_min，z_mid：

其中，max、min为最大值最小值函数；

步骤1.2.5如果z_max≥z_mid≥z_min，则窗口范围固定，进行窗口内像素点的滤波，进入步骤1.2.6；否则e＝e+2，并返回步骤1.2.3；

步骤1.2.6通过公式(5)计算K(x,y)，并更新像素矩阵K中点(x,y)处的灰度值；

其中，f_replace(x,y)＝z_mid，ν＝k(z_max-z_min)，Δ＝|f(x,y)-z_mid|，k为差值比例因子，K(x,y)图像K中点(x,y)的灰度值，≠表示不等关系，＝＝表示相等关系；

步骤1.2.7如果x≤n且y≤m，则

跳转到步骤1.2.2；否则输出的行扫描自适应中值滤波后的图像K；其中，m和n表示图像f的长和宽；

其中，≤表示小于等于关系，&表示与关系，≠表示不等关系，＝＝表示相等关系。

5.根据权利要求2所述的方法，其特征在于：所述1.3通过以下过程实现：

步骤1.3.1初始化(x,y)＝(1,1)，其中(x,y)为像素坐标；

步骤1.3.2通过公式(6)，计算(x,y)处像素梯度H(x,y),以及估计梯度方向

其中K(x,y)为图像K中点(x,y)的灰度值，atan为反正切函数；

步骤1.3.3如果x≤n且y≤m，则

并跳转到步骤1.3.2；否则，进入步骤1.3.4；

其中，≤表示小于等于关系，≠表示不等关系，＝＝表示相等关系，m和n表示图像K的长和宽；

步骤1.3.4初始化(x,y)＝(1,1)，其中(x,y)为像素坐标；

步骤1.3.5根据

的值，通过公式(7)计算像素点(x,y)处极大值抑制阈值H'(x,y)；

通过公式(8)，使用H'(x,y)对K(x,y)进行极大值抑制，计算K'(x,y)；

其中，K'(x,y)为对图像K进行极大值抑制后的图像中点(x,y)处的灰度值；

步骤1.3.6如果x≤n&y≤m，则

跳转到步骤1.3.5；否则输出的梯度极值筛选后的图像K′。

6.根据权利要求2所述的方法，其特征在于：所述1.4通过以下过程实现：

步骤1.4.1初始化(x,y)＝(1,1)，TH＝TH₀；其中(x,y)为像素坐标,TH是高阈值滤波的阈值，TH₀是其初始化的值；

步骤1.4.2通过公式(9)，计算K^*(x,y)；

其中，K^*(x,y)为图像K^*中(x,y)点处的灰度值；K'(x,y)为对图像K进行极大值抑制后的图像中点(x,y)处的灰度值；

步骤1.4.3如果x≤n&y≤m，则

并跳转到步骤1.4.2；否则，进入步骤1.4.4；

其中，≤表示小于等于关系，≠表示不等关系，＝＝表示相等关系，m和n表示图像K^*的长和宽；

步骤1.4.4初始化(x,y)＝(1,1)，i*＝i₀，初始化图像Q中所有像素点的值为0；

其中，(x,y)为像素坐标，i*为正方形搜索框的边长，i₀为其初值；

步骤1.4.5如果K^*(x,y)不为0，则Q(x,y)＝K^*(x,y)并进入步骤1.4.6，否则，跳入步骤1.4.9；

步骤1.4.6初始化(i,j)＝(0,0)，其中i,j矩阵搜索框的偏移量；

步骤1.4.7如果K^*(x+i,y+j)不为0，则Q(x+i,y+j)＝K^*(x+i,y+j)，并进入步骤1.4.9；否则，进入步骤1.4.8；

步骤1.4.8如果D(x+i,y+j)不为0，则Q(x+i,y+j)＝D(x+i,y+j)，并进入步骤1.4.9；其中，D(x+i,y+j)为差分模板D上点D(x+i,y+j)处的灰度值；

步骤1.4.9如果i≤i*&y≤i*则

并返回步骤1.4.7，否则进入步骤1.4.10；

步骤1.4.10如果x≤n&y≤m，则

并跳转到步骤1.4.5；否则，结束循环，输出轮廓图像Q。

7.根据权利要求1所述的方法，其特征在于：所述步骤(2)通过以下过程实现：

步骤2.1初始化q＝q₀,g＝g₀,R＝R₀,w*＝w₀,m*＝m₀,δ＝δ₀，Le＝Le₀并在图像Q中得到采样点i(a,b)，在每间隔q个像素点的区间中选择像素点灰度值最大的点作为采样点，并添加到采样点集合I中；

其中，q为采样间隔，g为细分数，g₀为其初值，R为兴趣半径，R₀为其初值，w*,m*为合并精度，w₀,m₀为其初值，δ为搜索精度参数，δ₀为其初值，Le为线段长度阈值，Le₀为其初值，i(a,b)∈I为图像Q中坐标为(a,b)的灰度值，I为采样集合；

步骤2.2随机在采样集合I抽样选择一点i(x,y)，在图像Q以(x,y)位置为中心，以区域半径R＝R₀的圆形为兴趣区域进行遍历，将该区域中的灰度值不等于零的像素点，灰度值置为最大，并添加到兴趣集合的点集Υ中；

步骤2.3通过公式(10)，将兴趣集合Υ中点映射为正弦曲线，计算正弦曲线在霍夫空间中的交点坐标，并添加到集合N^#中；

其中，(a,b)为集合Υ中点的坐标，(β,η)为霍夫空间中的正弦曲线的坐标值；

步骤2.4遍历集合N^#中的元素(ρ₀,θ₀)，初始化(ρ,θ)＝(ρ₀,θ₀)和k＝1，接着全部遍历集合N^#-{(ρ₀,θ₀)}中的元素，通过公式(11)进行依次迭代，计算(ρ,θ,k)，全部迭代完毕后将(ρ,θ,k)添加到重构集合N*中，遍历完成后删除重构集合N*中重复的元素；

其中，(ρ,θ,k)为合并后元素的值，(ρ,θ)为均值坐标，k为迭代过程中均值坐标的个数，(ρ₀,θ₀)为集合N^#中的元素，(ρ₁,θ₁)为集合N^#中除(ρ₀,θ₀)之外的元素，&表示与关系，|表示或的关系，{}表示集合；

步骤2.5在极坐标空间中对横纵坐标的区间进行细分，得到g组细分区间为[0,Δθ],[Δθ,2Δθ],[2Δθ,3Δθ]…[gΔθ,π],[0,Δρ],[Δρ,2Δρ],[2Δρ,3Δρ]…[gΔρ,ρ_max]；

其中g为细分数，ρ_max为集合N‘参数中的ρ的最大值，Δρ，Δθ分别为细分区间的长度；

步骤2.6初始化p₁＝0，p₂＝0,T＝T₀，其中，T为投票阈值，T₀为其初始值,p₁,p₂为细分区间的序号；

步骤2.7遍历N*中的元素,初始化n＝0，通过式子(12)，计算E(ρ',θ')并添加到集合E中,同时通过式子(13),据细分区间划分N*中的元素，计算V(ρ',θ')并添加到集合V中；其中，集合V为集合N*中均值坐标在所有细分区间内的均值坐标集合；

其中ρ'＝[p₁Δρ,(p₁+1)Δρ],θ'＝[p₂Δθ,(p₂+1)Δθ]，为细分区间，p₁,p₂为细分区间的序号，n为迭代过程中均值坐标在细分区间[p₁Δρ,(p₁+1)Δρ],[p₂Δθ,(p₂+1)Δθ]内的个数，E(ρ',θ')为集合N*中均值坐标在细分区间[p₁Δρ,(p₁+1)Δρ],[p₂Δθ,(p₂+1)Δθ内的个数，V(ρ',θ')为集合N*中均值坐标在细分区间[p₁Δρ,(p₁+1)Δρ],[p₂Δθ,(p₂+1)Δθ]内的均值坐标集合，∪为集合间的并集运算，∈表示属于符号，

为不属于符号，{}表示集合；

步骤2.8如果E(ρ',θ')＞T，则返回(ρ',θ')并进入步骤2.10,否则进入步骤2.9，其中，T为投票阈值；

步骤2.9如果p₁＜g且p₂＜g则

并跳转到步骤2.7，否则跳入步骤2.14；

步骤2.10将(ρ',θ')代入V中得到点集V(ρ',θ')，通过式子(14)计算搜索方向的直线参数ρ*,θ*；

其中(ρ_m,θ_m)是V(ρ',θ')中的j个点中的第m个点的坐标值，(ρ*,θ*)为极坐标描述直线的解析参数，即所描述方向的法线段的长度和其所在直线倾角；

步骤2.11通过式子(15)计算Ln(ρ*,θ*):在图像Q中沿着直线，搜索像素灰度值不为0的点的坐标，并添加到集合Ln(ρ*,θ*)中；

其中，(s,t)是搜索直线方向上点的坐标，Ln(ρ*,θ*)为图像Q中在解析参数为ρ*,θ*的直线上的像素点位置的集合，scosθ^*+tsinθ^*＝ρ^*为搜索直线的方程；

步骤2.12通过公式(16)计算更新参数ρ^,θ^；

其中，(ρ^,θ^)为拟合后直线的法线段的长度以及其所在直线倾角，为极坐标描述直线的解析参数；其中，(a_i,b_i)是Ln(ρ*,θ*)中n个元素中的第i个元素，a₀,b₀为笛卡尔坐标系下的直线的斜率和截距，atan为反正切函数；

步骤2.13通过公式(17)更新参数ρ*,θ*，并通过式子(18)将Ln(ρ*,θ*)添加到线段集合Ln中，并返回步骤2.9；

其中，(a_min,b_min),(a_max,b_max)是Ln(ρ*,θ*)中的横坐标最小和最大对应的两个的坐标值，Le为线段长度阈值；

步骤2.14删除集合I中的点i(x,y)，判断集合I是否为

若为空则结束算法，返回线段集合Ln，若不为空则跳入步骤2.2，继续循环。

8.根据权利要求1所述的方法，其特征在于：所述步骤(3)通过以下过程实现：

步骤3.1初始化E₁＝E₁₀，E₂＝E₂₀，T₁＝T₁₀，T₂＝T₂₀；

遍历线段集合Ln，通过公式(19)，找到所有满足条件的线段组，并统计所有线段组所在直线的参数，建立直线簇集合L；

其中，E₁，E₂为直线倾角差值的阈值，E₁₀ E₂₀为初值，T₁,T₂为聚类族数精度调整参数，T₁₀,T₂₀为其初值，L₁，L₂，L₃，L₄是线段集合Ln中满足要求的一组线段，

为L₁,L₂,L₃,L₄在Ln中对应的索引值，

为点集所在直线倾角的差值，∪表示取集合间的并运算，L为所有满足阈值条件的直线簇的解析参数的集合，{}表示集合，&表示与关系；

步骤3.2初始化r＝1，并计算Len＝length(L)；

其中，r为直线簇集合的索引，Len为集合L中直线簇的个数，length为集合长度计算函数；

步骤3.3将公式(20)代入公式(21)中，取

计算(x_A,y_A),(x_B,y_B),(x_C,y_C),(x_D,y_D)；

其中，mL₁,mL₂,mL₃,mL₄为直线簇L(r)中对应直线的标号，X(L*,C*),Y(L*,C*)为横纵坐标的计算函数，L*,C*为其输入参数，ρ_L*,θ_C*为直线参数，A,B,C,D为直线簇按照序号顺序相交的交点，(x_A,y_A),(x_B,y_B),(x_C,y_C),(x_D,y_D)为其坐标，L(r)为L中第r组直线簇参数，

为第r组直线簇中对应标号为mL₁,mL₂,mL₃,mL₄的直线的参数；

步骤3.4建立方程组(22)；

其中，u和v是方程组的解，

是二维平面中S₁，S₂，S₃，S₄四个点的坐标；

根据式子(23)依次取值并代入方程组(22)中，

依次求解方程组；若方程组的解都满足u＜0|v＜0|v+u＞1，则表示A,B,C,D四个点组成了凸四边形，跳入步骤3.6，否则进入步骤3.5；

步骤3.5如果r＜Len，则r＝r+1进入步骤3.3，否则进入步骤3.7；

步骤3.6通过公式(24)，将组成凸四边形的直线簇中各个直线的解析参数添加到凸四边形直线簇集合Q_LG中；

其中，{}表示集合，∪表示集合间的并运算；

步骤3.7计算Lo＝length(Q_LG)并初始化N＝2×Lo，N₁＝0；

其中N为聚类族群数，Lo为凸四边形直线簇集合Q_LG的长度，length为集合长度计算函数，N₁为聚类族群数的迭代参数；

步骤3.8遍历集合Q_LG中的元素得

接着全部遍历集合

中的元素，并通过公式(25)迭代更新N₁，所有遍历结束后计算

并进一步判断更新聚类族群数为

其中，

为集合Q_LG中的元素，{}表示集合；

为集合Q_LG中除了元素

之外的元素；

步骤3.9遍历Q_LG中的元素得

通过公式(26)粗略计算消影点的坐标(u₁,v₁),(u₂,v₂)，并添加到消影点集合M中；

其中，(u₁,v₁),(u₂,v₂)是凸四边形直线簇中对边直线得到的交点坐标，对边指的是

所描述的直线对和

所描述的直线对。

9.根据权利要求1所述的方法，其特征在于：所述步骤(4)通过以下过程实现：

步骤4.1随机从消影点集M中的选择N个点的坐标作为聚类族群中心点的坐标初值

开始聚类，并将循环次数初始化为k＝1，初始化γ＝γ₀，α＝α₀,β＝β₀,TG＝TG₀，γ为自适应率，γ₀为其初值,α,β为互补因子，α₀,β₀为其初值，TG为精度阈值，TG₀为其初值；

步骤4.2初始化G＝1；

步骤4.3通过公式(27)，计算M(G)到各个群类中心点的复合距离

其中，(X,Y)为M(G)的坐标，

为第k次迭代后第i类族群的中心点的坐标，i取值为1,2,3,...N，α+β＝1为互补因子，

为第k次迭代后点(X,Y)到第i类中心点的欧式距离，

为第k次迭代后(X,Y)到第i类中心点的余弦距离，

为第k次迭代后点(X,Y)到第i类中心点的复合距离，M(G)为M中第G个消影点；

计算过程中，通过公式(28)，调整互补因子α,β且0≤α,β≤1；

其中,γ为自适应率；

为边缘阈值，d＝min(m,n)，图像大小为m×n；

步骤4.4将距离

按大小排序，将M(G)归类到最近的族群类中；

步骤4.5判断如果G≥length(M)，进入步骤4.6，否则G＝G+1返回步骤4.3；

步骤4.6通过公式(29)重新计算族群的中心点；

其中X_kj,Y_kj，表示第k次迭代后i类中的第j个数据点，M_ki表示第k次迭代后i类中点的个数，class ki表示第k次迭代后i类的点的集合；

步骤4.7遍历第k次循环后的聚类中心

和第k-1循环后的聚类中心

TG为精度阈值，如果存在

则k＝k+1，并跳转到步骤4.2，否则退出聚类循环，完成运算，进入步骤4.8；其中，<>表示2-范数；

步骤4.8输出当前

的值，作为消影点集合M聚类后N个族群中心的位置。

Images