图1 特征点检测结果
Fig.1Results of feature points detection
因为铸件在生产过程中自身的工艺特点，使得成形后的铸件表面存在毛刺、飞边等外观缺陷，影响后续的铸件精加工和装配等。目前，铸件毛刺的检测和剔除主要依靠人工完成，工作量大，效率低，很难满足生产需求，因此，研究基于机器视觉的铸件毛刺自动检测与识别具有重要意义。
机器视觉具有广阔的应用市场，如机器人导航、工业自动化、恢复建筑原貌等[1]。按照工业相机的数量，机器视觉获取图像方法可分为单目视觉、双目视觉和多目视觉[2]。双目和多目视觉不仅需要准确标定相机之间的相对位置，还要保证拍摄时的同步性和稳定性，且拍摄空间大、拍摄过程难、成本较高[3-5]。单目视觉仅需通过一个相机围绕单个目标旋转拍摄一组图像序列，就可根据图像处理算法重建出三维物体或场景。单目视觉无需对相机相对于目标的位置进行标定，拍摄空间小，装置简单，成本低，易于实现[6]。
本文提出一种基于单目图像序列的三维重建方法，通过单目相机对铸件进行多个角度旋转拍摄获得一组图像序列，然后通过云计算对图像序列进行特征点检测与匹配，建立匹配点集合。利用从运动中恢复结构(structure from motion，SFM) 方法对匹配点进行稀疏点云重构和稠密点云重建获得三维点阵，再对点阵进行泊松表面重建，即可获得具有良好几何表面特征和细节的铸件三维重建表面。该方法有助于提高工业机器人对铸件进行精确加工的智能程度。
图像特征点检测与匹配是基于图像序列的三维重建的第一步，本文依据LOWE[7]提出的尺度不变特征变换算法( scale invariant feature transform，SIFT) 进行图像特征点的检测与匹配，并在此基础上提出“两步匹配法”对特征点进行优化，然后利用图像间的极线约束，采用随机采样一致算法(random sample consensus，RANSAC)提纯，得到铸件特征点优化匹配结果[8]。
SIFT特征点检测主要步骤如下：①尺度空间极值检测。搜索所有尺度上的图像位置，通过高斯微分函数来识别潜在的对尺度和旋转不变的兴趣点。②特征点定位。在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。③方向确定。基于图像局部的梯度方向，分配给每个特征点位置一个或多个方向。④特征点描述。采用该特征点邻域的梯度幅值大小和梯度方向来进行描述。以特征点为中心，用高斯窗确定特征点邻域的范围，高斯窗可用圆来代替。计算邻域内每个采样点的梯度幅值和梯度方向，然后将若干个采样点的计算值累加到n×n个统计窗口内，每个窗口用r表示梯度方向数，每个梯度方向的计算值为累加之后的梯度幅值大小。一个特征点有n×n×r维的特征向量描述子。根据文献[7]，当n=4，r=8时，特征点匹配率较高，如果n取更大值，则计算量增大，因此本文一个特征点有4×4×8=128维特征向量描述子。SIFT特征点检测方法不但增强了算法的抗噪能力，还为后续特征点匹配提供了容错性。图1所示为相邻两幅单目图像的特征点检测结果。
图1 特征点检测结果
Fig.1Results of feature points detection
一个特征点由一个128维的特征向量来描述，当相邻图像间的特征向量生成后，就可以进行特征点匹配。常见的特征点匹配方法是用特征点特征向量的欧氏距离来判断两幅图像中特征点的相似程度，即用特征向量的最近邻距离与次近邻距离的比值作为匹配约束条件。匹配原理是:Xa是图像1中的某个特征向量，Xb和Xb′是图像2所有特征向量中与Xa距离最近的两个向量，设Xa与Xb的距离小于Xa与Xb′的距离；如果最近邻距离与次近邻距离之比小于某个阈值T, 则Xb为Xa的匹配向量，表示为dist(Xa,Xb)/dist(Xa,Xb′)<T,保留这一对匹配向量所对应的特征点对；如果匹配向量不满足上述条件，就不保留特征点对，此步骤完成特征点的粗匹配。LOWE[7]推荐阈值T=0.8，但对存在尺度、旋转和亮度变化的两幅图片进行匹配时，T取值在0.4～0.6之间最佳，本研究对匹配点数目要求比较多，所以取T=0.6。粗匹配结果如图2a所示。
上述粗匹配存在的问题是：在匹配过程中由于欧氏距离只计算了两向量之间的距离，而忽略了各向量分量之间的相关性，所以会出现大量误匹配点对，数据冗余量大，后期计算量大。
本文提出“两步匹配法”,在将欧氏距离作为粗匹配约束条件的基础上，进一步通过向量空间余弦相似度匹配来增加匹配约束条件[9]，从而减少误匹配点对数。余弦相似度进一步考虑了上述粗匹配中两个特征向量之间的方向相关性，在向量空间中用这两个特征向量的夹角余弦值衡量特征点对差异的大小。余弦值越接近0，则表明两个特征向量夹角越接近于90°，也就说明此特征点对不相关，是误匹配点，予以剔除。Xa(Xa1,Xa2,…,Xa128)与Xb(Xb1,Xb2,…,Xb128)的余弦相似度为
(1)
首先在粗匹配基础上，判断Xa与Xb这两个特征向量的余弦相似度cosθ，如果大于某一经验阈值k，则进一步确定粗匹配点对为匹配点，否则粗匹配就被判断为误匹配点被剔除。利用余弦相似度匹配法可剔除大量误匹配点，如图2b所示。然后通过RANSAC算法优化提纯，得到最优匹配结果，如图2c所示。利用“两步匹配法”提高了特征点匹配的精度，剔除了大量误匹配点，提高了图像处理效率。
(a)特征点粗匹配结果
(b)余弦相似度匹配结果
(c)RANSAC优化匹配结果
图2 特征点匹配结果
Fig.2Results of feature point match
本文采用SFM方法[10]进行基于单目图像序列的铸件三维稀疏点云重建。SFM方法是一个迭代计算的过程。首先选取具有最多匹配点的图像作为初始图像对，计算得到初始的相机运动信息及铸件结构信息，通过两幅图像之间匹配点对的三角测量过程得到铸件的初始空间三维点云，重建的每个点都对应着被重构铸件的某一部分；然后逐渐增加新的图像进行迭代，SFM方法中每加入一幅新的图像就可以获得该图像所对应相机的运动信息以及铸件的结构信息；最后使用集束调整(bundle adjustment ，BA)[11]方法来进行参数优化[12-13]。
集束调整是一种非线性最小二乘法优化技术，通过估计相机姿态和图像序列之间的相对运动，以迭代方式添加和细化三维点。集束调整使用成本函数来最小化由运动结构三角测量预测的三维点与其相应观察位置之间的重投影误差，直至重建完整的图像集合，获得铸件的最优三维点云重建模型。重投影误差计算式为
(2)
其中，n个三维物点在m张图片中，Xij为第j张图片上看到的第i个物点,vij为第i个物点在第j张图片上的映射判断系数,若vij=1,表明有映射,若vij=0,则无映射。每张图片的相机j用向量aj参数化表示，每个三维物点i用向量bi参数化表示。函数Q(aj,bi)表示物点bi在相机aj下的投影坐标，也就是预测值。函数d表示观测到的图像坐标与预测的图像坐标之间的欧氏距离。
使用SFM方法生成的稀疏点云仅具有SIFT识别的特征点的三维坐标，所得到的点云是低密度的，不包含被重建铸件的更多细节，因此，还需要对稀疏点云进行稠密重建来获得稠密点云匹配。
本文采用多视角密集匹配(CMVS/PMVS)实现稠密点云匹配，其中点云多视角立体匹配(CMVS)[14]通过优化SFM输入的稀疏点云数据来对图像序列进行聚簇分类，以消减重建过程的数据量，提高运算速度和重建精度；然后基于面片模型的多视角立体匹配(PMVS)[15]根据聚簇分类后的图像序列和由SFM算法计算出的各自的投影矩阵，利用稀疏点云重构结果将点云向周围进行漫射，得到定向点云，同时利用局部光度一致性和全局可见性约束经过匹配、扩展、过滤等过程完成稠密重建过程，建立稠密点云模型。
PMVS是一种基于区域增长的方法：首先是匹配，在图像中采用Harris和Difference-of-Gaussian算法搜索图像特征f，在其他图像中搜索与f相匹配的图像特征f′，形成一个匹配子集F，利用点对(f，f′)构成三维点云，由点云再生成稀疏面片，然后将其保存在覆盖图像的网格中；其次是扩展，通过迭代算法在已经生成的稀疏面片邻近的空区域中，通过重复利用点云再生成新的面片，增加重构点的密度，获得稠密点云；最后是滤波，对重构点云通过两步滤波算法，增强可视密度并滤除不匹配误差。扩展和滤波步骤是交替迭代进行的，最终创建了一个致密的定向三维贴片云，获得三维稠密点云模型。
稠密点云模型比较形象地还原铸件模型，但仍是大量孤立的三维空间点集合，且表面还存在一些小的空洞。本文采用泊松曲面重构[16]方法，将对点云数据的曲面重构转化为求解泊松方程，通过求解稠密点云的最佳拟合表面，从而获得铸件连续光滑的三维表面。
基于泊松重构算法的曲面重建原理见图 3。泊松重构算法的实现过程主要包含对全局问题离散化、创建向量场、求解泊松问题、提取等值面4个方面。首先，定义一个函数空间，在此空间里对泊松问题进行离散处理，同时引入深度为D的八叉树对整个函数空间进行分割处理，计算每个点的特征向量以找到所有特征点表面法线所构成的有向点集；其次，定义一个指示函数
(内部为1，外部为 0)，用指示函数的梯度代替点云数据样本，采用三线性插值法来尽量提高八叉树叶子节点的精度，估计向量场V的近似方程；然后，以V为前提条件求解泊松方程▽x=▽·V，从而得到指示函数
之后，选取一个等价值，使等价值满足采样数据点分布在等值面的附近，计算指示函数与等价值相对应的等值面；最后，就可以通过移动立方体算法对等值面进行提取，完成点云的泊松表面重建。泊松表面重建的关键是，使指示函数的梯度▽x等于输入的点云法线向量所建立的向量场V。使指示函数的梯度▽x最佳逼近建立的向量场V，即minx‖▽x-V‖，这样指示函数的梯度场问题就可以转化为散度算子。由于向量场是不可积的，所以把散度算子简化为泊松方程。计算指示函数
向量场V的散度就等于指示函数的梯度场，即
▽x≡▽·▽x=▽·V
(3)
(a)有向点集 (b)指示函数梯度
(c)指示函数 (d)模型表面
图3 曲面重建示意图
Fig.3Diagram of curved surface reconstruction
泊松表面重建不但综合了全部的点云信息，而且无需具体细节就可以很好地过滤冗余信息。
图4为基于单目图像序列铸件三维重建方法的流程图。采用智能手机的摄像头围绕一个带有毛刺的铸件旋转拍摄42张图片，每张图片像素大小为4608×3456。旋转拍摄时只要求获得铸件的表面信息即可，对旋转没有严格要求。不要求以铸件为固定中心，不要求相机按照固定位置或间隔固定角度拍摄。
图4 铸件三维重建流程图
Fig.4 The flow chart of the 3D reconstruction of casting
图5为拍摄的部分图像。图6a为SFM方法迭代计算恢复的相机拍摄位置分布图；图6b为铸件稀疏点云重建结果，稀疏点云只是给出了铸件的一部分空间三维点阵，重建了部分轮廓。图7为铸件的稠密点云重建结果，稠密点云基本上重构了铸件的三维表面，一些明显的毛刺能显示出来，但是三维表面粗糙，外部轮廓局部点缺失明显，还存在一些明显的孔洞。图8所示为铸件的泊松表面重建结果，铸件的三维表面与稠密点云图像相比，重建点密集，表面光滑，表面细节基本修复，铸件上的毛刺的细节也能显示出来。实验结果表明，基于单目无标定拍摄图像序列的铸件三维重建方法能够实现基于机器视觉的铸件毛刺检测，为下一步铸件毛刺的自动配准与剔除提供了可行性。
(a)角度为0° (b)角度为8.32°
(c)角度为15.56° (d)角度为23.20°
(e)角度为33.54° (f)角度为40.37°
图5 图像序列中的部分图像
Fig.5 Some image sequences
(a)相机围绕铸件拍照的位置分布
(b)铸件稀疏点云重建效果
图6 相机的位置分布和铸件稀疏点云重建效果
Fig.6 The location distribution of the camera and the effect of casting sparse point cloud reconstruction
图7 铸件的稠密点云重建结果
Fig.7 Dense point cloud reconstruction results of casting
图8 铸件的泊松表面重建结果
Fig.8 Effect of the casting of poisson surface reconstruction
本文提出并实现了基于单目无标定拍摄图像序列的铸件三维重建方法。对图像序列进行特征点匹配中，通过余弦相似度匹配法增加对特征向量的约束条件，剔除大量误匹配点，压缩冗余数据，提高图像处理效率。以SFM迭代计算的方法从运动中恢复结构获得相机的位置信息以及铸件的稀疏点云，并重建了稠密点云。最后通过泊松表面重建,获得了铸件三维表面的细节特征。实验结果表明，本文提出的方法无需对相机的位置进行标定，所需拍摄空间小，对拍摄环境、镜头位置和图像序列数目也没有严格的要求，适用于机械加工环境下零件三维表面检测和毛刺识别。
[1] 方勇纯. 机器人视觉伺服研究综述[J].智能系统学报，2008, 3(2):109-114.
FANG Yongchun. A Review of Robot Vision Servo [J]. Journal of Intelligent Systems, 2008, 3(2):109-114.
[2] 佟帅，徐晓刚，易成涛，等.基于视觉的三维重建技术综述[J].计算机应用研究，2011, 28(7): 2411-2417.
TONG Shuai, XU Xiaogang, YI Chengtao, et al. A Review of 3D Reconstruction Based on Vision[J]. Computer Application Research, 2011, 28(7): 2411-2417.
[3] 林园胜.立体视觉匹配技术在机器人视觉中的应用研究[D]. 芜湖：安徽工程大学，2015.
LIN Yuansheng. Application of Stereo Vision Matching Technology in Robot Vision [D]. Wuhu:Anhui Polytechnic University, 2015.
[4] 杨阳.基于双目视觉的运动恢复结构[D].成都：电子科技大学，2016.
YANG Yang. Motion Recovery Structure Based on Binocular Vision [D]. Chengdu: University of Electronic Science and Technology, 2016.
[5] 张春秀.基于双目视觉的三维重建 [D]. 天津：天津大学，2010.
ZHANG Chunxiu. 3D Reconstruction Based on Binocular Vision [D]. Tianjin: Tianjin University，2010．
[6] 冉冉. 基于单目视觉的移动机器人目标识别及抓取系统研究[D]. 北京：北京交通大学，2010.
RAN Ran. Research on Target Recognition and Grasping System of Mobile Robot Based on Single Vision [D]. Beijing: Beijing Jiaotong University, 2010.
[7] LOWE D. Distinctive Image Features from Scale Invariant Features [J]. International Journal of Computer Vision, 2004, 60(2): 91-110.
[8] FISCHLER M A，BOLLES R C. Random Sample Consensus： a Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography[J]. ACM Graphics and Image Processing (S0001-0782), 1981, 24(6): 381-395.
[9] 张宇,刘雨东,计钊.向量相似度测度方法[J].声学技术，2009,28(4): 532-536.
ZHANG Yu, LIU Yudong, JI Zhao. Vector Similarity Measure Method [J]. Acoustic Technology, 2009, 28(4): 532-536.
[10] ULLMAN S. The Interpretation of Structure from Motion[C]//Proceedings of the Royal Society of London. London,1979: 405-426.
[11] TRIGGS B, MCLAUCHLAN P, HARTLEY R. Bundle Adjustment—a Modern Synthesis Vision Algorithms: Theory and Practice [J]. Computer Science, 2000, 1883: 153-177 .
[12] 沈武魁. 基于SFM的大幅面壁画图像高保真拼接技术[D].杭州：浙江大学，2011.
SHEN Wukui. SFM-based High Fidelity Splicing Technology for Large Format Mural Image[D]. Hangzhou: Zhejiang University, 2011.
[13] 张黎.基于立体视觉的三维建筑物重建技术研究[D].上海：上海师范大学，2012.
ZHANG Li.Research on 3D Building Reconstruction Based on Stereo Vision[D].Shanghai:Shanghai Normal University, 2012.
[14] 何豫航，岳俊. 基于 CMVS/PMVS 多视觉密集匹配方法的研究与实现[J]. 测绘地理信息，2013，38(3)：21-23.
HE Yuhang,YUE Jun. Research and Implementation of Multi-visual Intensive Matching Method Based on CMVS / PMVS[J]. Mapping Geographic Information，2013, 38(3):21-23.
[15] FURUKAWA Y，PONCE J. Accurate，Dense, and Robust Multi-view Stereopsis[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2010, 32(8):1362-1376.
[16] KAZHDAN M M, BOLITHO M，HOPPE H. Screened Poisson Surface Reconstruction[J].Eurographics Symposium on Geometry Processing, 2006, 32(3):61-70.
开放科学(资源服务)标识码(OSID)：