在日常生活、工业生产中,圆柱类的物体很常见,例如机械行业的心轴、法兰和直钢管等。通过视觉图像中的几何特征实现对圆柱类物体的空间位姿估计和尺寸测量是机器人视觉中的一项重要研究内容^[1-3]。现有的基于圆柱投影的视觉定位研究,多是利用圆柱底面投影椭圆进行位姿估计^[4-5],也有少数是利用圆柱侧面的外轮廓投影直线对进行估计^[1]。其中,由直线轮廓和椭圆轮廓组成的圆柱透视投影轮廓识别精度的高低,将直接影响到圆柱空间位姿估计的精度。

一般情况下,圆柱透视投影轮廓由其上下底面分别投影形成的两个椭圆和侧面投影形成的一对直线段组成。有关独立的直线或椭圆的图像识别和拟合都已有大量的研究和成熟的算法,常见的是采用Hough变换方法进行检测识别。基于Hough变换的直线识别方法,只需识别两个参数^[6]。基于Hough变换的椭圆识别方法,由于需要识别5个参数,采用标准Hough变换法所需的存储空间较大,运算复杂度高,因此现有研究的着眼点主要是利用椭圆的几何特性降低复杂度、提高运算效率^[7-10],或通过概率方法提高运算效率。另一类常用的基于随机抽样一致(random sample consensus,RANSAC)的椭圆识别和拟合方法则对存储空间没有较高要求,并且对含噪声数据集表现出较好的鲁棒性。SONG等^[11]提出了一种包含3点椭圆拟合法的伪RANSAC方法,提高了椭圆识别的实时性和鲁棒性。XIE等^[12]采用了一种改进的RANSAC方法得到更精确的模型估计。 Martelli等^[13]设计了一种基于现场可编程门阵列(field programmable gate array,FPGA)的RANSAC方法以提高椭圆识别的实时性。近年来还出现了一些用于识别椭圆的新方法,如Guerrero等^[14]采用基于扩展Kalman滤波器的方法识别血管轮廓和圆参数,Thévenaz等^[15]通过迭代修正活动轮廓使区域间对比度最大化来得到椭圆轮廓。

然而,在存在遮挡(椭圆有部分不可见)和图像背景含有大量噪声的情况下,分别进行直线识别和椭圆识别的鲁棒性变差,并且缺乏去除噪点即筛选内点的有效方法。如何更有效利用图像数据,是提高投影图像轮廓参数识别精度的关键。Doignon等^[1]在单视图中利用退化的圆锥方程,实现了无需预先将圆柱投影轮廓中的两条直线的数据点进行分组,直接拟合两条直线。Taubin^[16]则提出了估计隐式方程定义的复杂平面曲线、曲面以及空间曲线的方法。在以上研究的基础上,本文提出将圆柱透视投影轮廓用一个高次隐式方程进行描述的方法。这种方法的优势在于可以方便引入投影轮廓中固有的约束,即直线对和椭圆之间的对称和相切约束关系,更有效地利用图像数据,改善曲线拟合精度。

1 隐式方程描述

常见的圆柱类物体,特别是一些机械零件,常常是以圆柱组合体的形式存在。由于组合体中相邻圆柱的直径差异,待识别的圆柱底面存在可见或不可见的可能。一般情况下,存在两种可能：即上下底面均可见,或只有一个底面可见。对于这两种情况,所需建立的高次隐式方程分别为6次和4次。本文先讨论4次隐式方程的情况。

为简化推导过程,先忽略曲线过原点的情况,记直线的一般方程为

$Ax+By+1=0.$

(1)

椭圆的一般方程记为：

$\begin{align} & A{{x}^{2}}+Bxy+C{{y}^{2}}+Dx+Ey+1=0, \\ & {{B}^{2}}-AC < 0. \\ \end{align}$

(2)

对于圆柱只有一个底面可见的情况,可用一个4次隐式方程表示由侧面轮廓投影的直线组和底面轮廓投影的椭圆组成的圆柱投影轮廓曲线,

$\begin{align} & \left( {{q}_{1}}x+{{q}_{2}}y+1 \right)\left( {{q}_{3}}x+{{q}_{4}}y+1 \right)\bullet \\ & \left( {{q}_{5}}{{x}^{2}}+{{q}_{6}}xy+{{q}_{7}}{{y}^{2}}+{{q}_{8}}x+{{q}_{9}}y+1 \right)=0. \\ \end{align}$

(3)

根据式(2),式(3)中对应于椭圆的参数应满足约束

$q_{6}^{2}-{{q}_{5}}{{q}_{7}}<0.$

(4)

直线组与椭圆满足相切条件,每条直线与椭圆有且仅有一个交点,即直线方程式(1)和椭圆方程式(2)联立消元后的一元二次方程判别式等于0,于是可得到两个等式约束：

$\beta _{^{i}}^{2}-4{{\alpha }_{i}}{{y}_{i}}=0,i=1,2.$

(5)

其中：

$\begin{align} & {{\alpha }_{i}}={{q}_{5}}-\left( {{q}_{2i-1}}{{q}_{6}} \right)/{{q}_{2i}}+\left( {{q}_{2i-1}}^{2}{{q}_{7}} \right)/q_{2i}^{2}, \\ & {{\beta }_{i}}={{q}_{8}}-\left( {{q}_{2i-1}}{{q}_{9}}+{{q}_{6}} \right)/{{q}_{2i}}+\left( 2{{q}_{2i-1}}{{q}_{7}} \right)/q_{2i}^{2}, \\ & {{\gamma }_{i}}={{q}_{7}}/q_{2i}^{2}-{{q}_{9}}/{{q}_{2i}}+1. \\ \end{align}$

同时,直线组还满足轴对称关系,其对称轴与椭圆的短轴共线,根据这个关系也可得到一个等式约束,

$\frac{{{c}_{1}}+{{c}_{2}}}{1-{{c}_{1}}{{c}_{2}}\pm \sqrt{\left( c_{1}^{2}+1 \right)\left( c_{2}^{2}+1 \right)}}-c=0.$

(6)

其中：c为椭圆短轴倾角的余切,c₁1和c₂是直线组中两条直线各自倾角的余切,即：

$\begin{align} & c=cot\left( \frac{1}{2}ar\tan \frac{{{q}_{6}}}{{{q}_{5}}+{{q}_{7}}}+\frac{\pi }{2} \right), \\ & {{c}_{1}}=-\frac{{{q}_{2}}}{{{q}_{1}}}, \\ & {{c}_{2}}=-\frac{{{q}_{4}}}{{{q}_{3}}.} \\ \end{align}$

因此,得到了1个非线性不等式约束式(4)和3个非线性等式约束式(5)、(6)。

下面对上述推导未涉及的情况进行补充说明：

1) 对于6次隐式方程的情况。隐式方程的左侧仅需在式(3)的基础上再增加一个形为式(2)左侧的因子,相应地约束增加1个非线性不等式约束和3个非线性等式约束即可。

2) 对于直线过原点的情况。只需将图像整体平移后再进行识别即可回避该问题。

2 拟合问题描述和算法

对于任意2维点x=(x₁1,x₂)和参数向量q=(q₁1,…,q₉),定义映射$f:{{}^{2+9}}\to \gg $,

$\begin{align} & f\left( q,x \right)=\left( {{q}_{1}}{{x}_{1}}+{{q}_{2}}{{x}_{2}}+1 \right)\left( {{q}_{3}}{{x}_{1}}+{{q}_{4}}{{x}_{2}}+1 \right)\bullet \\ & \left( {{q}_{5}}{{x}_{1}}^{2}+{{q}_{6}}{{x}_{1}}{{x}_{2}}+{{q}_{7}}{{x}_{2}}^{2}+{{q}_{8}}{{x}_{1}}+{{q}_{9}}{{x}_{2}}+1 \right). \\ \end{align}$

当给定一组参数q时,存在映射φ_q:$${\varphi _q}:{ \gg ^2} \to \gg $ $,

${{\varphi }_{q}}\left( x \right)=f\left( q,x \right).$

再定义集合Z(φ_q)={x|φ_q(x)=0}为φ_q的化零集合,该集合的点组成了以隐式方程φ_q(x)=0定义的曲线。

任意点x到Z(φ_q)的近似最短距离的平方^[16]为

$\delta {{\left( q,x \right)}^{2}}=\frac{{{\varphi }_{q}}{{\left( x \right)}^{2}}}{{{\left\| \nabla {{\varphi }_{q}}\left( x \right) \right\|}^{2}}}=\frac{f{{\left( q,x \right)}^{2}}}{{{\left\| \nabla f\left( q,x \right) \right\|}^{2}}}.$

(8)

对于曲线拟合问题,式(8)中的q为待拟合的参数,是变量,而x为用于拟合的数据点,是已知的常数。当待拟合的数据点集为D={p₁1,…,p_s},所有点到曲线的距离的平方和可表示为

$\Delta q=\sum\limits_{i=1}^{s}{\delta {{\left( q,{{p}_{i}} \right)}^{2}}.}$

(9)

因此,这个曲线拟合问题可以描述为含约束的非线性优化问题：

$\begin{align} & \Delta q=\sum\limits_{i=1}^{s}{\delta {{\left( q,{{p}_{i}} \right)}^{2}}.} \\ & \min \Delta \left( q \right) \\ & s.t.q_{6}^{2}-{{q}_{5}}{{q}_{7}}<0; \\ & \beta _{i}^{2}-4{{\alpha }_{i}}{{\gamma }_{i}}=0,i=1,2; \\ & \frac{{{c}_{1}}+{{c}_{2}}}{1-{{c}_{1}}{{c}_{2}}\pm \sqrt{\left( c_{1}^{2}+1 \right)\left( c_{2}^{2}+1 \right)}}-c=0. \\ \end{align}$

该优化问题可利用罚函数方法,构造辅助函数将问题转化为无约束问题来求解^[17]。

记:

$\begin{align} & g\left( x \right)=q_{6}^{2}-{{q}_{5}}{{q}_{7}}, \\ & {{h}_{1}}\left( q \right)=\beta _{1}^{2}-4{{\alpha }_{1}}{{\gamma }_{1}}, \\ & {{h}_{2}}\left( q \right)=\beta _{2}^{2}-4{{\alpha }_{2}}{{\gamma }_{2}}, \\ & {{h}_{3}}\left( q \right)=\frac{{{c}_{1}}+{{c}_{2}}}{1-{{c}_{1}}{{c}_{2}}\pm \sqrt{\left( c_{1}^{2}+1 \right)\left( c_{2}^{2}+1 \right)}}-c=0. \\ \end{align}$

构造辅助函数

$F\left( x,\sigma \right)=\Delta \left( q \right)+\sigma \left( \sum\limits_{j=1}^{3}{h_{j}^{2}\left( x \right)+{{\left[ \max \left\{ 0,g\left( x \right) \right\} \right]}^{2}}} \right).$

(10)

其中,σ是一个很大的正数。则原问题可转化为

$minF\left( x,\sigma \right).$

这是一个无约束的非线性优化问题。

现将该曲线拟合算法描述如下：

1) 初值选取。采用常规方法(如Hough变换方法,甚至用手工勾画方法亦可)快速粗略地识别出圆柱图像轮廓的侧面直线组和底面椭圆,计算其相应的参数作为初值,建立参数的初始估计q_est。

2) 内点筛选。根据参数的初始估计q_est,采用式(8)计算原始数据集中每一个点x到Z(f)的近似最短距离平方δ(q,x)²,筛选出δ(q,x)²小于阈值t的数据点作为内点,而大于t的点作为外点被排除。

3) 生成目标函数。用上一步筛选出的所有内点组成待拟合的数据点集D,并代入式(9)生成以q为变量的目标函数。

4) 解算。求解辅助函数式(10)形式的无约束非线性优化问题,得到待拟合的曲线参数q_fit。

3 实验验证

为了验证本文提出的方法,采用仿真数据和实物图像分别进行了实验。

3.1 仿真数据

仿真数据由一系列在Matlab中编程生成的数据点组成。通过给定一组曲线参数的理论值,分别作出一系列属于曲线上的点和属于Gauss白噪声的点,从而得到含有噪声的待拟合点集。

进行拟合时,为了验证本文方法提高参数识别精度的效果,选取了一组与理论值q_t存在一定偏差的初始估计参数,按照本文方法进行拟合,其结果见图 1和表 1。图 1中将曲线的初始估计(估计曲线)和拟合得到的曲线(拟合曲线)对比,估计曲线与理论曲线点存在一定距离,拟合曲线显著地减小了与理论曲线点的距离,提高了曲线参数识别的精度。表 1中的数据显示了曲线参数的初始估计和拟合结果以及分别相对于理论值的误差。其中,相对误差为绝对值,初始估计的相对误差e_re与拟合结果的相对误差e_rf计算方法如下：

$\begin{align} & {{e}_{re,i}}=\left| \frac{{{q}_{est,i}}-{{q}_{t,i}}}{{{q}_{t,i}}} \right|\times 100%, \\ & {{e}_{rf,i}}=\left| \frac{{{q}_{fit,i}}-{{q}_{t,i}}}{{{q}_{t,i}}} \right|\times 100%, \\ \end{align}$

图 2对比了参数的初始估计的相对误差e_re与拟合结果的相对误差e_rf。可以看到,各参数的e_rf都小于e_re,表明拟合结果的参数精度得到了提高。

3.2 实物图像

本文提出的方法在真实油气井口管件的特征提取中得到了应用。现场拍摄的图像为图 3a,其中感兴趣的是位于最下方的法兰和连接法兰的管道,即图 3a中矩形框标记的部分。其特点是背景中的大量碎石形成噪点,并且圆柱投影轮廓中的椭圆有部分不可见。为了便于进行定量分析,采用了由一大一小两个同轴圆柱组合而成的实物模型,采集其图像进行实验。为了引入噪声,背景采用表面粗糙的平板。拍摄的图像如图 3b所示。

对图 3b中的两个圆柱体图像,先用Canny算子得到含有噪点的轮廓点图像,再分别用直线和椭圆识别算法获得轮廓曲线的初始估计。其中：采用概率Hough直线方法^[6]获取轮廓中的直线部分,采用文^[18]的方法快速获取轮廓中的椭圆部分。

通过本文方法进行轮廓曲线拟合时,由于图 3b中的圆柱均为上下底面的椭圆轮廓都可见,故采取6次隐式方程及相应的约束描述轮廓曲线并拟合,拟合结果如图 4a—4h所示。其中：图 4a和4e分别是上下两个圆柱的含噪轮廓点图像; 图 4b和4f是拟合得到的曲线,其中的矩形框是示意局部放大的区域; 图 4c和图 4c是初始估计曲线的局部放大图; 图 4d和4h是拟合得到曲线的局部放大图。对比图 4c和4d可以看到,图 4c中的椭圆初始估计与图像中的轮廓存在偏差,经过本文方法拟合得到的图 4d中该偏差显著减小了。对比图 4g和4h,肉眼可见图 4c中的椭圆初始估计与直线的初始估计相交,椭圆偏大,经本文方法拟合得到的图 4h中轮廓曲线的椭圆部分与直线部分相切,并且与图像中的轮廓重合度提高。

本文方法旨在通过曲线拟合得到精度较高的轮廓,从而为圆柱类物体的3维重建提供更为有效的输入,使空间位置和尺寸重建精度得到提高。由于图像中物体的轮廓真实值较难获取,本文采用从拟合得到的轮廓曲线重建的物体尺寸与真实的物体尺寸进行对比的方式,对曲线拟合的精度进行定量评价。

图 3b采用的实物模型中：从上至下第1个圆柱的半径为39.9 mm,高度为70.0 mm; 第2个圆柱的半径为54.6 mm,高度为56.3 mm。已知空间圆的半径,可以利用单目图像中的椭圆重建计算出空间圆心的3维坐标^[19]。计算圆柱上下两个顶面圆心的距离,可得到圆柱的高度。表 2用圆柱高度的理论值作为参考,对比了曲线拟合前后两个圆柱高度分别重建的结果。其中：第1个圆柱的高度重建计算结果,经曲线拟合后绝对误差约减小为拟合前的20.30%,第2个圆柱的高度重建计算结果绝对误差则约减小为拟合前的41.67%。

4 结 论

针对遮挡和背景含噪点较多条件下对圆柱透视投影轮廓的提取存在鲁棒性差、精度低的情况,本文提出一种含约束的隐式方程曲线对其进行描述及拟合的方法。由于合理利用了组成投影轮廓的直线和椭圆之间的固有约束,该方法能在用常规方法分别提取的直线和椭圆轮廓存在一定误差的情况下改善结果的精度。仿真实验和实物图像验证反映出了该方法具有较强的鲁棒性并且能有效改善轮廓曲线的精度。

本文提出的方法适用于对圆柱透视投影轮廓提取的鲁棒性和精度要求较高的场合,例如采用2维图像特征进行视觉定位或测量等。该方法在进行内点筛选时阈值的取值会对运算的收敛性产生较大影响,如何更有效地进行内点筛选,值得进一步研究。