CN100430690C - 利用单数码相机自由拍摄进行物体三维测量的方法 - Google Patents

Abstract

一种利用单数码相机自由拍摄进行物体三维测量的方法，属测试技术。该方法包括测量准备、图像摄取、编码点的识别与定位、相机位姿确定、目标曲线提取、同名曲线的自动匹配优化、目标曲线的三维重建七个步骤。其特征在于：三维测量用一个数码相机、一台电脑、一组编码点和一把标尺，测量时先对被测物体上需要测量的目标曲线进行标记，以利于图像识别；在被测物体周围放置标尺和一组编码点；再手持一个数码相机以自由拍摄方式获得被测物体的一组图像；根据这组图像，自动精确计算各次拍摄时的相机位置与姿态；提供使用者方便的交互手段，实现对标识曲线的半自动提取和不同图像中同名曲线的优化匹配，从而自动计算出所标识曲线结构上的三维点列信息。

Description

利用单数码相机自由拍摄进行物体三维测量的方法

一、技术领域

三维物体测量属于测量、测试技术领域。在国际分类表中的相应代码为G01B。

二、背景技术

三维测量技术在逆向工程、工业检测、质量控制等领域日益广泛的应用需求极大地推动了三维测量技术的迅速发展，出现了基于光学、声学、电磁学以及机械接触原理的各种测量方法，如三坐标测量机、激光扫描仪、结构光测量仪等。其中，三坐标测量机采用机械接触式传感，测量精度高，但一般需要专门的测量室和专用测量台，对测量环境要求很高，测量范围有限，测量效率低，不适合软性物体的测量；激光线扫描测量和结构光照射测量是目前三维几何外形测量的主流方法，通过在物体表面进行激光或结构光照射，可以快速获取模型表面密集的点云数据，但都受扫描范围、物体表面高光反射等的限制，且激光扫描仪和结构光测量仪均价格昂贵。尤其是对于机械产品而言，由于通常含有明显的结构特征，常常需要重点获取对重建被测物体的数字化模型起关键作用的角点、棱边以及模型表面的某些控制线等数据，结构光照射测量法、激光线扫描测量法等获得的都是表面的整体点云或网格数据，一方面输出的数据量十分庞大，另一方面却不能直接显式地获取所需的棱边特征及某些关键截面控制线数据。而且，这些方法输出的测量数据通常在模型的光滑平坦区域的效果比较好，而恰恰在关键的角点和棱边处的测量效果较差。

为了以更加简单的硬件条件和更加灵活方便的方式实现三维几何信息的测量，根据一个数码相机拍摄的多幅图像精确重建物体的位置和形状成为近年来的一个研究热点。其中，德国Gom公司的统，已经能够用单数码相机自由拍摄方式进行较高精度的三维坐标定位。该系统通过在场景中放置一组编码点和长度标尺，并在感兴趣的部位粘贴易于识别的标记点，然后使用者手持一个数码相机自由拍摄多幅图像，要求各幅图像相互之间有一定的重叠，所有图像输入相应的软件系统后，系统一次性自动计算各次拍摄时的相机位置、姿态以及所有标记点的空间坐标。该系统已经商品化并在我国进行销售，但是，该系统目前只能进行特定标记点目标(由一个黑色圆环围绕一个白色圆点构成)的空间坐标定位，一般用于配合其它测量方法进行多视角测量数据的拼合，而不能进行曲线目标的三维测量，不能用于具有复杂几何外形产品的三维数字化模型重建。

三、发明内容

本发明旨在用简单的硬件条件，实现一种面向工业产品测量建模的便于实施、精度较高、成本较低的实用测量方法。为此，本发明通过对被测物体的特征线和进行数字化模型重建所需要的物体表面的某些关键截面控制线进行标记，使其在颜色亮度上明显区别于被测物体本身的颜色，以利于图像识别；在被测物体周围放置一把标尺和一组经特殊设计的编码点；然后手持一个数码相机以自由拍摄方式获得被测物体的一组图像；根据这组图像，自动精确计算各次拍摄时的相机位置与姿态，同时，提供使用者方便的交互手段，实现对标识曲线的半自动提取和不同图像中同名曲线的优化匹配，进而自动计算出所标识曲线结构的三维点列信息。

根据上述方案开发出了使用灵活方便、适用于不同大小物体上三维曲线结构测量的实用系统，可以很好地用于基于实物的机械产品三维数字化模型建立。本发明提出的物体三维测量方法，其特征在于测量只用一个数码相机、一台普通个人电脑，辅以一组编码点和一把标尺，无需复杂的测量硬件系统，也不需要对测量系统进行繁琐的标定；直接显式地产生对物体进行数字化模型重建所需的曲线数据，避免数据冗余，且便于进行高效率模型重建；全部测量数据自动位于一个世界坐标系下，避免了其他测量方法中多次测量的数据需要拼合的难点问题，也避免了多个数据集拼合带来的累计误差。本发明包括测量准备、图像摄取、相机位姿确定、目标曲线提取、同名曲线的自动匹配优化、目标曲线的三维重建等主要步骤。

测量准备与图像摄取

测量准备主要包括三个方面的工作：1)根据数字化模型重建的需要对目标曲线(一般是自然曲面片的边界线、关键的截面控制线等)进行标记，使其在颜色亮度上明显区别于被测物体，以利于图像识别；2)在测量区域内布置若干编码点。每个编码点具有唯一身份编码，易于在不同图像中进行快速可靠的识别。标记点和编码点均可简单地在计算机上产生相应的图案，然后打印生成。把编码点贴到硬纸板、木片等的上面，可以重复使用；3)在测量场景内放置一个有两个编码点的标尺，其上的两个编码点中心的间距已知。放置标尺的目的是获得被测物体的实际尺寸，否则只能得到相差一个比例系数的三维结构。上述准备工作完成后，就可手持数码相机对被测物体进行多角度拍摄。要求各幅图像间具有一定的重叠，即一幅图像至少要与另一幅图像之间有共同可见的5个以上编码点和某些目标曲线。

编码点的识别与定位

编码点的图案设计成由中心白色圆点、中间黑色圆环、外圈圆环组成，其中外圈圆环被分成15个等分，每个等分为黑色或白色，黑色表示二进制码“0”，白色表示二进制码“1”，称为“编码带”。一组编码点中每个点的编码均不相同。根据这种编码，可以在不同图像中可靠地识别其身份，自动建立同名编码点在各图像间的对应关系。正是根据这些同名编码点在图像中的位置和在多幅图像间的对应，才能实现自由拍摄时相机的各个位置和姿态的自动计算。

由于编码点的中心圆形经CCD成像后呈椭圆状，因此，为了识别和定位图像中的编码点，本发明首先采用Canny算子进行图像分割，在图像中提取代表不同区域的轮廓信息，然后根据各轮廓的尺寸、形状、椭圆拟合残差、区域灰度的均值、区域灰度的方差等5个约束条件，对候选编码点目标进行逐步过滤，从而实现编码点目标的提取。

锁定一个编码点目标后，要对其进行解码，即确定其具体是哪个编码点。解码的依据是编码带上各扇区的灰度。本发明采用拟合编码带中间椭圆并对其上的每一象素取线性窗口进行中值滤波的方法，综合考虑了编码带内大多数象素的灰度值，能够消除孤立噪声的影响。对大量的现场拍摄图像的处理结果表明，这一方法对提高编码点身份识别的鲁棒性十分有效。

最后，根据识别出的编码点中心圆形区域内各象素点的灰度值，确定具有亚象素定位精度的编码点中心坐标。

相机位置和姿态确定

根据至少5个同名编码点中心在两幅图像中的象素坐标，首先计算两幅图像之间的基本矩阵。由相机内参数和基本矩阵，可进一步恢复两幅图像对应的相机姿态及这两幅图中共同可见的编码点中心的三维坐标。然后，根据所获得的三维空间点与第三幅图像上编码点之间的同名对应关系，求解第三幅图像对应的相机姿态，进而获得更多的编码点中心三维坐标，再求解下一幅图像对应的相机姿态，如此递增，直到获得所有相机姿态和编码点中心三维坐标。最后采用光束平差法同时对所有相机参数及编码点中心三维坐标进行整体优化，以进一步提高精度。这种递增式方法和全局优化方法相结合的策略，使算法既有较高的效率，又能使相机定位达到较高的精度。

一次读入拍摄到的一组图像，测量系统即自动计算并记录各次拍摄时的相机位置与姿态。确定了各次拍摄时的相机位置和姿态，意味着确定了各幅图像在一个统一的世界坐标系下的位置和姿态，这使得后续算法根据不同图像对重建出的各目标曲线直接位于同一个坐标系中，无需数据拼合。

目标曲线的半自动提取

对于当前选中显示在两个图像窗口中的一个图像对，分别在两幅图像上用鼠标在共同可见的同一条标记曲线附近取点，使这些点的连线大致反映相应的图像曲线轮廓，然后测量系统根据能量优化的原理自动将这样的大致图像曲线轮廓最佳地贴合到图像曲线上。这一半自动提取过程中，使用者只需要在图像曲线附近依次取点即可，简便易行，且由于有较好的初始搜索位置而大大增加了目标曲线提取及后续自动匹配优化的稳定性，同时最佳贴合方法保证目标曲线的提取精度。

同名曲线的自动匹配优化

在一个图像对中提取出一对同名目标曲线后，一个关键的问题就是建立两幅图像间同名目标曲线上各象素点的对应关系。基于图像点的对应，就可以通过立体三角法重建出这些点的空间坐标。

根据立体视觉的基本理论，相机在不同位置和角度获得的一个图像对间的同名点应满足极线约束。对于一个图像对中两个候选同名匹配点v₁和v₂，本发明用在第二幅图像上v₂到v₁的极线的距离，以及在第一幅图像上v₁到v₂的极线的距离之和来度量v₁和v₂的匹配程度。同时，两条同名目标曲线上的点列匹配应满足空间相关性约束，即一条曲线上按顺序排列的点列，一定对应另一条曲线上一个顺序排列的点列。基于这样的分析，本发明首先采用动态规划来获得同名曲线上离散象素点的初始匹配。得到同名图像曲线上点对的初始匹配后，再进一步进行曲线的匹配优化。设两条同名图像曲线分别用参数方程c₁(l)、c₂(l)表示，本发明优化下面的目标函数来达到c₁、c₂曲线上点的精确匹配

$\min \underset{0}{\overset{L_1}{\∫}}\frac{\left|c_2{\left(\mathrm{\σ}\left(l\right)\right)}^T{\mathrm{Fc}}_1\left(l\right)\right|}{|\left|e{\mathrm{Fc}}_1\left(l\right)\right||}+\frac{\left|c_2{\left(\mathrm{\σ}\left(l\right)\right)}^T{\mathrm{Fc}}_1\left(l\right)\right|}{\left|\right|c_2{\left(\mathrm{\σ}\left(l\right)\right)}^T{\mathrm{Fe}}^T\left|\right|}\mathrm{dl}---\left(28\right)$

其中， $e=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right];$ F是图像对间的基本矩阵，在各次拍摄时的相机位置和姿态已经确定的情况下时已知的；σ(l)为待求的映射函数，表示c₁曲线上参数为l的点在曲线c₂上的参数值。

目标曲线的三维重建

完成同名曲线上所有象素点的匹配后，由于图像对在各自拍摄时的相机位置和姿态已经自动计算得出，因此可以用双目立体视觉中成熟的三角测量法重建出这些点的空间坐标，从而完成整条曲线的重建。

本发明具有测量硬件简单(一个数码相机、一台普通个人电脑、一把标尺、一组打印生成的编码点)、测量方式十分灵活(自由拍摄)、测量范围不受限制、各角度测量数据自动拼合、测量输出数据无冗余、使用方便、成本低等突出优点，不但能够用于空间点的定位，而且能够根据机械产品测量建模的需要，测量出被测物体上的棱边、特征线、关键截面控制线等三维信息，在逆向工程、产品质量检测等领域有广泛的应用前景。

附图说明

图1：本发明提出的测量方法的基本流程图。

图2：编码点示意图。图2(a)描述编码点结构，即由中心白色圆点、中间黑色圆环、外圈圆环组成，其中外圈圆环被分成15个等分，根据每个扇区分的颜色确定其身份，黑色表示二进制码“0”，白色表示二进制码“1”；图2(b)是三个编码点示例。

图3：实施例测量系统软件图形界面分区示意图。1.菜单区；2.图标工具栏；3.显示图像文件的列表；4.显示当前活动图像对中的一幅图像；5.显示当前活动图像对中的另一幅图像；6.重建出的目标曲线的三维图形显示区。通过点击区域3中的图像文件列表来指定4、5两个窗口内当前显示的图像，用户在4、5两个图像窗口中依照图像上的标识曲线交互勾勒出一对同名曲线的大致形状后，系统即自动计算出该曲线上的三维点列并显示于三维图形区。

具体实施方式

本发明提出的物体三维曲线结构测量方法的实施例说明如下：数码相机采用具有内置闪光灯、分辨率为4256×2848的Nikon手动调焦数码相机，计算机采用主频2.8GHz、内存512MB的Pentium IV微机，测量软件系统在Visual C++6.0平台上实现。

根据图1叙述本发明的具体实施方式及原理：测量之前要进行一定的测量准备，包括对被测物体上需要测量的目标曲线(一般是物体的特征线和进行数字化模型重建所需要的物体表面的某些关键截面控制线)进行标记，使其在颜色亮度上明显区别于被测物体本身的颜色，以利于图像识别；在被测物体周围放置一把标尺和一组经特殊设计的编码点，这组编码点中每个点的编码均不相同，即每个编码点具有身份唯一性。完成上述测量准备工作以后，手持一个数码相机以自由拍摄方式获得被测物体的一组图像，要求各幅图像间具有一定的重叠，即一幅图像至少要与另一幅图像之间有共同可见的5个以上编码点和某些目标曲线。根据这组图像，测量系统自动识别并精确定位各幅图像中的编码点，然后自动精确计算各次拍摄时的相机位置与姿态。使用者经简单的交互完成当前活动图像对中一条标识曲线的半自动提取，测量系统随后自动进行同名图像曲线的优化匹配，进而自动计算这条出所标识曲线结构的三维点列信息。如果尚有未处理的目标曲线，则对下一条目标曲线(可能出现在不同的图像对中)重复上述半自动提取、自动匹配、三维曲线重建的过程，直至完成所有目标曲线的三维重建。以下就图1中的几个主要步骤的具体实施方式做详细说明。

编码点身份识别

本发明的三维测量方法基于对拍摄得到的一组图像进行分析处理，首先就是识别图2中所示的编码点。每个编码点的中心为圆形“目标点”，周围为环状“编码带”，“编码带”按照角度平均分为15份，每24度一份，相当于一个二进制位，取白色为前景色，相应的二进制码为“1”，黑色为背景色，相应的二进制码为“0”。对于每一个编码点，存在15种可能的二进制编码，取这15个二进制数里最小的数对应的十进制数作为编码点的ID。

本发明的编码点自动检测算法主要包括以下三个主要过程：①编码点目标的提取，即在图像中寻找“目标点”；②根据编码点“编码带”上的信息确定编码点的唯一身份，即编码点解码；③编码点中心的亚象素定位。

(1)编码点目标提取

编码点的中心圆形经CCD成像后呈椭圆状。因此，首先采用Canny算子进行图像分割，在图像中提取代表不同区域的轮廓信息，然后采用逐步过滤的方法进行编码点目标提取。首先根据标记点目标的尺寸、形状对可能的目标点进行初步过滤，满足如下条件的封闭轮廓进入进一步的识别过程：

P_min≤P≤P_max(1)

1≤P²/4πA≤1.5(2)

式中，P和A分别为封闭轮廓的周长和面积，P_min、P_max分别是轮廓周长的最小和最大阈值。式(1)是对封闭轮廓大小的限定，式(2)则度量了其与圆的接近程度。

对于满足(1)(2)式的封闭轮廓，采用最小二乘方法进行椭圆拟合，剩余误差ε_eli满足给定允差ε_τ方为候选编码点的中心圆，即

ε_eli≤ε_τ    (3)

经过最小二乘模板匹配后，图像中所有椭圆轮廓都已找到。但是，在真实场景中往往存在一些不是标记点目标而具有椭圆形状或者与椭圆形状接近的轮廓也被误认为是标记点目标，但由于本方法采用的标记点目标前景灰度为白色，背景灰度为黑色，两者对比度强烈。这是标记点目标区别于其它非标记点目标的显著特征。因此，要根据标记点的这种灰度特征来进一步排除非编码点目标。由于已经通过式(3)所示的椭圆拟合准则，因此可以确定内部中心椭圆和黑色圆环的区域。记中心白点内部区域的灰度均值为M_I，黑色圆环区域的灰度均值为M_O，则M_I和M_O应满足：

$\left\{\begin{array}{c}{M}_I\≥M_t\\ M_O\≤M_t\\ M_I-M_O\≥\mathrm{\Δ}{M}_t\end{array}\right.---\left(4\right)$

其中，M_t为区分前景灰度与背景灰度的阈值；ΔM_t为前景灰度与背景灰度之差应满足的最小值。

另外，约束中心白点内部区域的灰度方差V_I和黑色圆环区域的灰度方差满足：

$\left\{\begin{array}{c}{V}_I\≤{\mathrm{\δ}}_I\\ V_O\≤{\mathrm{\δ}}_O\end{array}\right.---\left(5\right)$

其中，δ_I，δ_O是允许的最大灰度方差。条件(5)约束编码点的中心必须满足一定的灰度均匀性。满足上述(1)～(5)式，则进入编码点解码过程。

(2)编码点解码

本发明提出的编码点解码算法的具体实现步骤如下：

Step1：拟合编码点中心圆点的外轮廓椭圆(记为椭圆A)、中间黑色圆环的外轮廓椭圆(记为椭圆B)、各白色扇区所在的圆环的外轮廓椭圆(记为椭圆C)。再拟合一个位于椭圆B、C中间的椭圆，其中心及旋转角与椭圆B、C的中心及旋转角相同，长、短轴分别取B、C长、短轴的均值。再采用椭圆绘制的算法，获得椭圆D上各象素点的位置坐标。

Step2：计算A所包围的区域内所有象素灰度的中值作为前景灰度，A和B之间区域内所有象素灰度的中值作为背景灰度。前景灰度与背景灰度的均值作为阈值，用于后续确定编码点各二进制位的码值。

Step3：对于椭圆D上的任一个象素点TD，过椭圆中心做一条射线，记该射线与椭圆B、C的交点为TB和TC。将线段TBTC上所有象素的灰度值进行排序，取最中间的一个象素灰度值作为TD的新灰度值。

Step4：对椭圆D上各点按照公式(6)做逆仿射变换，使椭圆D对应一个单位圆，单位圆上各点的灰度对应于椭圆B的新灰度值。

$X^{\′}=\left[\begin{array}{cc}{a}^{-1}& \\ & b^{-1}\end{array}\right]\left[\begin{array}{cc}\cos \mathrm{\α}& \sin \mathrm{\α}\\ -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right](X-X_o)---\left(6\right)$

式中，X′是与TD对应的单位圆上点的坐标，X是点TD的坐标，X_o是椭圆D的中心O的坐标，a，b分别是椭圆D的长轴和短轴的长度，α是椭圆D的旋转角。

Step5：对单位圆上像素作二值化，取其中的一个边缘点作为起始点。

Step6：从起始点开始，在该单位圆上每隔24°为一个二进制位，计算各位内所有象素点的平均灰度值。若某位的灰度平均值大于阈值，则该位取二进制码为“1”；否则取“0”。从而可以得到编码点的一个二进制编码。找到与该二进制数对应的最小十进制数，该十进制数即为编码点的ID。

上述Step3相当于对点TD在一个线形窗口内进行一次中值滤波，滤波窗口内包含的象素即为线段TBTC上的象素。采用中值滤波后椭圆D上新的灰度值来确定编码点每一位的二进制码，考虑了编码带内所有象素的灰度值，能够消除孤立噪声的影响。对大量的现场拍摄图像的处理结果表明，这一方法对提高编码点身份识别的鲁棒性十分有效。

(3)标记点中心计算

采用式(7)来进行编码点的中心亚象素定位

$\left\{\begin{array}{c}{x}_c=\underset{j}{\mathrm{\Σ}}\underset{i}{\mathrm{\Σ}}i\·I_{i,j}/\underset{j}{\mathrm{\Σ}}\underset{i}{\mathrm{\Σ}}{I}_{i,j}\\ y_c=\underset{j}{\mathrm{\Σ}}\underset{i}{\mathrm{\Σ}}j\·I_{i,j}/\underset{j}{\mathrm{\Σ}}\underset{i}{\mathrm{\Σ}}{I}_{i,j}\end{array}\right.---\left(7\right)$

式中，(x_c，y_c)为编码点中心坐标，I_i，j为中心圆形区域内象素点(i，j)的灰度值。

相机位姿自动确定

在齐次坐标表示下，三维空间点X在摄像机成像平面上的投影x可表示为：

x＝K[R|t]X＝PX    (8)

其中，K是摄像机内参数矩阵；R和t分别为从世界坐标系到摄像机坐标系的旋转变换矩阵和平移变换向量；P为3×4的投影变换矩阵。

设相机在两个不同的位置和朝向对同一场景拍摄了两幅图像，两个相机之间存在旋转矩阵R₁₂和非零平移向量t₁₂，则由极线几何可知两幅图像之间存在以下约束

$x_2^{T}F{x}_1=0---\left(9\right)$

式中x₁和x₂分别是三维空间点X在第一、第二两幅图像上的投影点；F是3×3的基本矩阵，它映射右图像上的一点x₂到左图像相应的对极线Fx₁上。

根据两幅图像之间编码点的对应关系 $x_1^l\↔x_2^l,$ l＝1，L，N，N≥5，我们首先采用MLESAC(Maximum Likelihood Estimation SAmple Consensus)方法计算两幅图像之间的基本矩阵F。

根据相机标示的焦距等内参数，可以构建内参数矩阵K的初值(产品标示参数只作为初值，测量系统在后续过程中会对其做进一步优化)。这样，根据基本矩阵F，可进一步计算两幅图像之间的本质矩阵E

E＝K^TFK    (10)

又根据本质矩阵的定义E＝[t]_xR(式中[g]_x表示向量的反对称矩阵)，利用旋转矩阵的正交性可以很容易推导出

${\hat{E}}^T\hat{E}=\left[\begin{array}{ccc}1-{\hat{t}}_x^2& -{\hat{t}}_x{\hat{t}}_y& {-\hat{t}}_x{\hat{t}}_z\\ -{\hat{t}}_y{\hat{t}}_x& 1-{\hat{t}}_y^2& -{\hat{t}}_y{\hat{t}}_z\\ -{\hat{t}}_z{\hat{t}}_x& -{\hat{t}}_z{\hat{t}}_y& 1-{\hat{t}}_z^2\end{array}\right]---\left(11\right)$

式中 $\hat{E}=E/\sqrt{\mathrm{Tr}\left(E^{T}E\right)/2},$ Tr(g)表示矩阵的迹， ${\hat{t}}_{12}=t_{12}\left|\right|t_{12}\left|\right|$ 是归一化后的平移向量。这样，根据(10)式求出的E矩阵和式(11)可以很容易求得归一化后的平移向量 ${\hat{t}}_{12}={({\hat{t}}_x,{\hat{t}}_y,{\hat{t}}_z)}^T.$ 由于 ${(-\hat{E})}^T(-\hat{E})={\hat{E}}^T\hat{E},$ 因此，归一化得到的矩阵

可能与实际的

相差一个符号。另外，由于矩阵

的每一个元素都是关于向量

的二次项，因此根据式(11)计算出的向量

也具有二义性，即

都满足(11)式。后文将给出解决

和

的二义性的方法。

为了计算第一幅图与第二幅图之间的旋转矩阵R₁₂，定义

$w_i={\hat{E}}_i\×{\hat{t}}_{12},(i=\mathrm{1,2,3})---\left(12\right)$

式中

代表矩阵的各行向量。设r_i是旋转矩阵R₁₂的各行向量，则

r_i＝w_i+w_j×w_k    (13)

式中(i，j，k)是(1，2，3)的循环组合。这样就确定了前两幅视图的相机姿态。

将世界坐标系建立在第一个相机上，根据成像几何关系不难导出空间点X在第一个相机坐标系下的Z方向坐标为

$Z_1=f\frac{{({\mathrm{fr}}_1-x_2{r}_3)}^T{\hat{t}}_{12}}{{({\mathrm{fr}}_1-x_2{r}_3)}^T{x}_1}---\left(14\right)$

进一步可求得另外两个坐标分量

X₁＝x₁Z₁/f，Y₁＝y₁Z₁/f    (15)

X在第二个相机坐标系下的坐标为

X₂＝R₁₂(X₁-t₁₂)    (16)

由于

和

的二义性，可能产生四对不同的

根据拍摄时的实际情况，只有当采用某个对重建的所有点均同时在两个相机前面，即所有点(这里是指两个视图共同可见的各个编码点中心)的Z₁和Z₂都为正时，才表明重建结果正确，相应的一组

即为正确解。

由于上述重建算法中，两个相机之间基线长度未知，因此只能得到归一化的平移向量

从式(14)容易看出，重建出来的场景与实际场景相差一个固定的比例因子。为此，在场景中放置标尺，标尺上两个标记点之间的距离已知，从而可以确定一个比例因子，得到被测物体的实际尺寸。

在两视图相机姿态确定和编码点中心三维坐标重建的基础上，需要进一步依次确定其它各次拍摄时的相机姿态。对第j幅图像进行处理时，要求该图像中存在至少6个已在前面的步骤中重建出了三维坐标的编码点，即已知图像点和空间点的对应 $X_i\↔x_i,$ i＝1，L L，L≥6。将这些约束带入投影方程(8)中可得

x_i＝P_jX_i，i＝1，L L，L≥6    (17)

由于每组 $X_i\↔x_i$ 的对应产生两个线性方程，因此，根据(17)式采用最小二乘法即可求解出第j幅图像的投影矩阵P_j中的11个未知元素。

将3×4的投影矩阵P_j表示成

P_j＝K[R_j|t_j]＝[KR_j|Kt_j]＝[M|p₄]    (18)

其中，M是矩阵P_j前3×3的子矩阵，p₄表示矩阵P_j的第四列。从式(18)中很容易确定出平移向量t_j＝K^-1p₄。

由于内参数矩阵是上三角的且旋转矩阵是正交的，因此，对矩阵进行QR分解就可得旋转矩阵R_j。在估计出拍摄第j幅图像时相机相对于世界坐标系的外部姿态参数R_j和t_j的基础上，我们进一步通过光学三角形法，重建出第j幅图像中新出现的且可在前j-1幅图像中找到同名匹配的编码点的三维空间点坐标。至此，就完成了对应第j幅图像的相机姿态确定和编码点三维坐标计算。然后继续下一幅图像的处理，直至处理完所有图像。

(3)相机姿态优化

由于图像噪声等因素的存在，三维空间点X_i经投影矩阵P_j变换后的像点，与实际识别出的X_i在第j幅图像中的像点坐标x_ij并不重合。为了进一步提高系统精度，本发明基于光束平差算法，以再投影误差最小建立目标函数

$\underset{\mathrm{ij}}{\mathrm{\Σ}}d{(P_j{X}_i,x_{\mathrm{ij}})}^2\→\min ---\left(19\right)$

对前面求出的相机参数和三维空间点坐标进行全局优化。具体求解采用LM(Levenberg-Marquardt)算法。由于是以前面求出的已经比较接近真实值的X_i和P_j作为初值，因此，全局优化可以较快收敛。

目标曲线的半自动提取

本发明开发的物体曲线结构三维测量软件系统图形界面示意图如附图3所示，该图的左侧显示所有拍摄得到的图像文件的列表，该图的右侧下部为当前活动图像对的两个图像显示窗口，通过点击左侧的图像文件列表来指定两个窗口内当前显示的图像，该图的右侧上部为重建出的目标曲线的三维图形显示区。

本发明采用能量优化的基本思想将在两个活动图像上交互勾勒出的目标曲线的大致轮廓通过后台算法自动“贴合”到相应的目标曲线上。

具体实现时，首先连接用户输入的目标曲线上的点构成一条折线(封闭情况下为多边形)，然后本发明采用计算机图形学中直线段光栅扫描转换的DDA算法快速得到折线经过的所有象素点，从这些象素点中按固定间隔(实施例中每隔两个象素取一个)进行

采样，记为v_i，i＝(0，1，L，n)。这里需要进一步利用Canny边缘检测算子自动检测出的图像边缘信息，因此记检测出的边缘点集为P。另记v_ij，j＝(1，L，8)为点v_i的八邻域象素点，同时为了叙述方便，记v_i0＝v_i。

本发明建立每个点v_i及其八邻域点的如下能量函数：

E(v_ij)＝[αE_tension(v_ij)+βE_bend(v_ij)+γE_img(v_ij)+δE_attr(v_ij)]    (20)

其中，E_tension(v_ij)，E_bend(v_ij)，E_img(v_ij)，E_attr(v_ij)分别为v_ij点处的拉伸能量、弯曲能量、图像能量和边缘点引力产生的能量，α，β，γ，δ分别为各能量项的权值，用以调节各能量项的比重。为了平衡各项的影响，各能量项全都归一化到区间[0，1]：

$E_{\mathrm{tension}}\left(v_{\mathrm{ij}}\right)=\frac{|\stackrel{\‾}{d}-|v_{\mathrm{ij}}-v_{i-1}\left|\right|}{\underset{0\≤j\≤8}{\max }\left\{\right|\stackrel{\‾}{d}-|v_{\mathrm{ij}}-v_{i-1}|\left|\right\}},$ $\stackrel{\‾}{d}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}|v_i-v_{i-1}|$

$E_{\mathrm{bend}}\left(v_{\mathrm{ij}}\right)=\frac{{|v_{i-1}-2v_{\mathrm{ij}}+v_{i+1}|}^2}{\underset{0\≤j\≤8}{\max }\left\{{|v_{i-1}-{2v}_{\mathrm{ij}}+v_{i+1}|}^2\right\}}$

$E_{\mathrm{img}}\left(v_{\mathrm{ij}}\right)=\frac{\underset{0\≤j\≤8}{\min }\left(E_{\mathrm{im}g}\left(v_{\mathrm{ij}}\right)\right)-E_{\mathrm{img}}\left(v_{\mathrm{ij}}\right)}{\underset{0\≤j\≤8}{\max }\left(E_{\mathrm{img}}\left(v_{\mathrm{ij}}\right)\right)-\underset{0\≤j\≤8}{\min }\left(E_{\mathrm{img}}\left(v_{\mathrm{ij}}\right)\right)}$

$E_{\mathrm{attr}}\left(v_{\mathrm{ij}}\right)=\frac{|v_{\mathrm{ij}}-p_{\mathrm{ij}}|}{\underset{0\≤j\≤8}{\max }\left\{\right|v_{\mathrm{ij}}-p_{\mathrm{ij}}\left|\right\}}$

其中， $E_{\mathrm{img}}=-{|\▿I(x,y)|}^2,$ E_attr(v_ij)中的p_ij∈P是与v_i距离最近的边缘点。本发明在拉伸能、弯曲、图像能之外增加边缘点引力能E_attr(v_i)的目的是进一步促使点列向目标曲线收敛。p_ij的搜索限定在以v_i为中心的一个窗口内进行。如果该窗口内没有任何边缘点，E_attr(v_ij)＝0，j＝0，L，8。

通过迭代的方法使点列v_i，i＝(0，1，L，n)向使式(20)能量极小的位置移动，最终锁定在图像特征附近，形成光滑的目标曲线点列。再用均匀B样条曲线来拟合这些点，以备后续曲线匹配使用。对于开曲线，约束其两端点始终保持位置不变，从而防止开曲线的点列收缩为一点产生退化。

同名曲线的自动匹配优化

(1)匹配程度度量

相机在不同位置和角度获得的一个图像对间同名点应满足的一个基本约束是(9)式描述的极线约束。在相机内参数和相对位置及姿态已知的情况，图像对间的基本矩阵F是已知的。对于一个图像对中两个候选匹配点v₁和v₂，可以用在第二幅图像上v₂到v₁的极线的距离

$D_2(v_1,v_2)=\frac{\left|v_2^T{\mathrm{Fv}}_1\right|}{\left|\right|eF{v}_1\left|\right|}---\left(21\right)$

以及在第一幅图像上v₁到v₂的极线的距离

$D_1(v_1,v_2)=\frac{\left|v_2^T{\mathrm{Fv}}_1\right|}{\left|\right|v_2^T{\mathrm{Fe}}^T\left|\right|}---\left(22\right)$

度量v₁和v₂的匹配程度，其中， $e=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right].$ 本发明正是根据式(21)和式(22)的极线约束条件，建立同名曲线间点对匹配的优化目标函数。

(2)图像曲线重采样

为了得到拟合出来的均匀B样条曲线p(u)，u∈[0，1]上的各象素点，取离散参数间距Δu＝1/L，L为各型值点的累加弦长，得到均匀B样条曲线上的离散象素点v_i，i＝0，1，...，N，构造分段线性插值曲线满足：

$\left(23\right),\left\{\begin{array}{cc}c\left(l_i\right)=v_i,& i=0,...,N\\ c\left(l\right)=\frac{l_{i+1}-1}{l_{i+1}-l_i}{v}_i+\frac{l-l_i}{l_{i+1}-l_i}{v}_{i+1},& l_i\&le;l<l_{i+1}\end{array}\right.$

其中l₀＝0.0， $l_i=\underset{j=1}{\overset{i}{\mathrm{\&Sigma;}}}|v_j-v_{j-1}|.$ 记 $L=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}|v_j-v_{j-1}|.$

不失一般性，将两条同名图像曲线中象素点较多的一条记为c₁(l)，较少的那条记为c₂(l)。c₁(l)上的离散象素点记为v_k ⁽¹⁾，0≤k≤N₁；c₂(l)上的离散象素点记为v_j ⁽²⁾，0≤j≤N₂。寻求c₁(l)，l∈[0，L₁]上的各象素点在c₂(l),l∈[0，L₂]上的匹配点，可以实现其在c₂(l)上的亚象素匹配，从而实现高精度三维曲线重建。

(3)匹配优化

首先采用动态规划来初始匹配对应曲线上的离散象素点。同名曲线一个对应点对的累积代价函数定义为

$C(v_k^{\left(1\right)},v_j^{\left(2\right)})=D(v_k^{\left(1\right)},v_j^{\left(2\right)})+\underset{m\&Element;G_{\mathrm{kj}}}{\min }C(v_{k-1}^{\left(1\right)},v_m^{\left(2\right)})---\left(24\right)$

其中 $D=(v_k^{\left(1\right)},v_j^{\left(2\right)})=D_1(v_k^{\left(1\right)},v_j^{\left(2\right)})+D_2(v_k^{\left(1\right)},v_j^{\left(2\right)}),$ G_kj表示在v_k ⁽¹⁾与v_j ⁽²⁾已经匹配的情况下m的所有可能取值。由于两条同名图像曲线上的象素点个数不等，动态规划匹配的结果会出现多对一的情况，即较长的c₁曲线上的多个象素点对应较短的c₂曲线上的同一个象素点。

得到图像曲线上点对的初始匹配后，再进行曲线的匹配优化。根据式(21)和式(22)，本发明优化下面的目标函数来达到c₁、c₂曲线上点的精确匹配

$\min \underset{0}{\overset{L_1}{\&Integral;}}\frac{\left|c_2{\left(\mathrm{\&sigma;}\left(l\right)\right)}^T{\mathrm{Fc}}_1\left(l\right)\right|}{\left|\right|{\mathrm{eFc}}_1\left(l\right)\left|\right|}+\frac{\left|c_2{\left(\mathrm{\&sigma;}\left(l\right)\right)}^T{\mathrm{Fc}}_1\left(l\right)\right|}{\left|\right|c_2{\left(\mathrm{\&sigma;}\left(l\right)\right)}^T{\mathrm{Fe}}^T\left|\right|}\mathrm{dl}---\left(25\right)$

其中σ(l)为待求的映射函数，表示c₁曲线上参数为l的点在曲线c₂上的参数值。将式(25)的积分形式改写为求和的形式

(26)式水，l₀，l₁，L，分别是v_k ⁽¹⁾，0≤k≤N₁对应c₁(l)上的参数值。以上述动态规划法产生的粗匹配为迭代初始条件，采用共轭梯度法即可求解式(26)的最小化问题。

曲线三维重建

完成同名曲线上所有象素点的匹配后，由于相机内参数已知，且图像对在各自拍摄时的相机位置和姿态已经自动计算得出，即相对于世界坐标系的旋转矩阵和平移向量也已知，因此可以用双目立体视觉中成熟的三角测量法重建出这些点的空间坐标，具体计算式可表达为

$\left.\begin{array}{c}{x}^{\left(j\right)}=f\frac{R_{11}^{\left(j\right)}X+R_{12}^{\left(j\right)}Y+R_{13}^{\left(j\right)}Z+T_x^{\left(j\right)}}{R_{31}^{\left(j\right)}X+R_{32}^{\left(j\right)}Y+R_{33}^{\left(j\right)}Z+T_z^{\left(j\right)}}\\ y^{\left(j\right)}=f\frac{R_{21}^{\left(j\right)}X+R_{22}^{\left(j\right)}Y+R_{23}^{\left(j\right)}Z+T_y^{\left(j\right)}}{R_{31}^{\left(j\right)}X+R_{32}^{\left(j\right)}Y+R_{33}^{\left(j\right)}Z+T_z^{\left(j\right)}}\end{array}\right\},j=\mathrm{1,2}---\left(27\right)$

其中，f是相机的焦距，x^(j)和y^(j)分别是同名匹配点在第j幅图像上的象素坐标的两个分量，R.^(j)是第j次拍摄时相机相对于世界坐标系的旋转矩阵的各个分量，T.^(j)是第j次拍摄时相机相对于世界坐标系的各个平移分量，根据(27)式中的4个方程用最小二乘法即可解出该点的空间坐标的3个未知分量(X，Y，Z)。对同名曲线上每一个匹配点对执行上述求解过程，就完成了整条目标曲线点列的三维测量。完成所有目标曲线的三维测量后，可以进一步根据这些点列构造出曲线的参数方程及模型的曲面方程。

Claims (3)

1.一种利用单数码相机自由拍摄进行物体三维测量的方法，其特征在于，包括测量准备、图像摄取、编码点的识别与定位、相机位姿确定、目标曲线提取、同名目标曲线的自动匹配优化、目标曲线的三维重建七个步骤，具体步骤在于：首先以被测物体的特征线和进行数字化模型重建所需要的物体表面的关键截面线为目标曲线，对这些目标曲线进行标记，使其在颜色亮度上明显区别于被测物体本身的颜色，以利于图像识别；在被测物体周围放置一把标尺和一组编码点；然后手持一个数码相机以自由拍摄方式获得被测物体的一组图像，该组图像中的每一幅图像至少与另一幅图像之间有共同可见的5个以上编码点；根据这组图像，自动精确计算各次拍摄时的相机位置与姿态；然后使用者只需要分别在两幅图像上用鼠标在共同可见的同一条目标曲线附近取点，使这些点的连线大致反映相应的图像曲线轮廓，然后测量软件系统迭代地使图像曲线轮廓的能量达到最小，从而自动将交互勾勒的大致图像曲线轮廓最佳地贴合到图像曲线上；采用动态规划法来获得同名目标曲线上象素点的初始匹配，再通过非线性优化方法使同名目标曲线上所有匹配点对到各自极线的距离之和达到最小，实现图像对上同名目标曲线点列的最优匹配；进而根据立体三角原理计算出点列的三维坐标。

2.根据权利要求1所述的利用单数码相机自由拍摄进行物体三维测量的方法，其特征在于，所述编码点的识别，是设计的一组编码点中每个点的编码均不相同，采用尺寸、形状、椭圆拟合残差、区域灰度的均值、区域灰度的方差5个约束条件，对候选编码点目标进行逐步过滤；编码点解码过程中综合考虑编码带内大多数象素的灰度值，来达到消除噪声影响、提高编码点身份识别的鲁棒性的效果；采用目标区域内灰度加权的方法进行编码点中心亚象素定位。

3.根据权利要求1所述的利用单数码相机自由拍摄进行物体三维测量的方法，其特征在于，所述相机位姿确定，是根据编码点的身份唯一性，自动建立同名编码点在各图像间的对应关系；根据至少5个同名编码点中心在两幅图像中的象素坐标，恢复两幅图像对应的相机姿态及这两幅图中共同可见的编码点中心的三维坐标，进而采用递增式方法和全局优化方法相结合的策略，实现所有图像对应的相机位置和姿态的自动求解。

Images