CN104376552A - 一种3d模型与二维图像的虚实配准算法 - Google Patents

Abstract

本发明为了解决二维图像与三维图形配准的虚实融合中计算效率低和需要预先标定相机内参数的问题，提供一种3D模型与二维可见光图像配准算法。它通过有多视几何关系的多幅图像，根据多视几何理论，用估计加优化的方法求解多个摄像机内参数（焦距）和外参数（外参数包括摄像机视点位置及旋转方向）；再通过GPS位置变换三维恢复坐标下的摄像机视点到3D模型坐标的视点，这个3D模型是根据真实场景预先建好的；最后渲染显示三维模型的二维虚图，它是与二维图像是一致的。本方法可以将视点所观察到的三维场景（虚图）与图片配准，它只需预先测量得到摄像机位置，不需要进行摄像机的预标定，也不需要用户手动设置初始的模型对应于各幅二维图像的摄像机视角。相对传统配准方法，有较高的效率。

Description

一种3D模型与二维图像的虚实配准算法

技术领域

本发明涉及计算机应用技术和计算机视觉中的增强现实领域，具体涉及一种3D模型与二维图像的虚实配准算法。 

背景技术

增强现实是指通过将计算机生成的虚拟场景、文字注释等信息实时、精确地叠加到使用者所观察到的真实世界景象中，对人的视觉***进行延伸和扩充的一种技术。如何实时、精确地计算用户视点相对于真实世界的位置、姿态信息，并利用这些信息将虚拟场景正确地叠加到它所应处的位置上，即虚实配准，是增强现实***中的关键问题之一。 

所以，虚实配准在增强现实中有极大的应用需求。同样，在三维空间进行二维图像与三维图形的虚实配准技术，也有着广泛的应用前景，但是它没有通用成熟配准框架。 

在三维空间进行虚实配准，常用的技术框架是通过多幅图片形成的多目视觉***得到场景的稠密点云，通过稠密点云提取场景的三维几何特征，如平面、直线、特征点、曲线等，再与3D模型中相应的特征进行配准。这种方法的缺点一是生成稠密点云需要的图片多，一般要几十张图片，而在机场应用场景中难以布置如此多指向同一场景的摄像机；二是计算速度慢，根据情况不同，至少在10分钟以上。工程应用中，稠密点云常采用激光雷达高速测量的方法，一分钟可测得上百万的点云数据，但是激光雷达成本高，它主要用于场景3D模型建模，难与视频结合的虚实配准结合使用。并且，现有的二三维配准研究中，基本上都要求用户手动设置初始的配准参数，即模型的初始视角位置，然后再优化求解精确配准位置，如果不手动设置视角初始值，配准问题将变得无法求解。 

目前的解决方法，还有的需要预先标定相机。预标定方法可以提前测量摄像机所有的内参数，但在有的应用环境中（如：机场，山区等）较难测量实际标志物的三维世界坐标，即摄像机的外参数，尤其是摄像机旋转参数；针对可变焦距可旋转的PTZ（Pan-Tilt-Zoom）相机就无法应用。 

  

发明内容

本发明要解决的技术问题是：为了解决二维图像与三维图形配准的虚实融合中计算效率低和需要预先标定相机内参数的问题，提供一种3D模型与二维可见光图像配准算法。此算法不需要预标定相机，而且计算时间在两三分钟以内，可以达到实用的要求。 

本发明解决其技术问题所采用的技术方案是：一种3D模型与二维图像配准算法，它通过有多视几何关系的多幅图像，用估计加优化的方法求解多个摄像机内参数（焦距）和外参数，外参数包括摄像机视点位置及旋转方向；再变换三维恢复坐标下的摄像机视点到3D模型坐标的视点，这个3D模型是根据真实场景预先建好的；最后渲染显示三维模型的二维虚图，它是与二维图像是一致的。 

本发明的创新点在于：1）不需要像以往方法一样预先标定相机的内参数；2）不需要其它激光雷达等复杂设备进行场景的三维测量，也不需要限制相机的个数，只需要简单测量摄像机的位置，而摄像机的焦距和旋转方向靠多视几何关系自动计算出来；3）提出了一个自动配准的算法解决3D模型与二维图像的配准。 

算法具体步骤是： 

1）通过有多视几何关系的多幅图像求解摄像机内参数和外参数。

1-1)  提取SIFT特征描述（它具有尺度和旋转不变性）。 

1-2）再用KD-树进行特征点匹配。 

1-3）采用RANSAC框架通过8点法估计基础矩阵F，同时去掉了噪声匹配点对。 

1-4) 选择共有特征点最多的两幅图像，再通过上步得到的基本矩阵F，按估计的摄像机参数算出本质矩阵E，计算出第一对摄像机之间的相对位置，得到摄像机位置旋转和平移参数。 

1-5) 初始化摄像机内参数，将上述外参数以及内参数、所有特征点的投影坐标和世界坐标作为初始值，用稀疏Levenberg-Marquardt方法进行集束优化。 

1-6) 再选择另一幅图片，根据已算出的空间特征点及特征点跟踪关系用RANSAC方法求解此摄像机投影矩阵，得到摄像机内外参数初始值，将其加入到优化框架下，转步骤1-5）继续优化。 

1-7) 反复执行1-6）直到每一个摄像机（所有的图片）处理完毕，得到所有摄像机的外参数（多幅图片视点的旋转平移的相对位置） 

2）变换三维重建坐标下的摄像机视点到3D模型坐标的视点。

根据预先好的测量摄像机的经纬度坐标找到多个视点的相对位置到3D模型的相对变换关系（旋转、平移、缩放），可以将视点位置变换到3D模型中。我们采用的方法是经纬度转换法，它分为2个步骤： 

2-1) 将测量的摄像机经纬度坐标转换成3D模型坐标系中的XYZ坐标。

2-2) 设摄像机位置的归一化的相对三维坐标为b _i ，经纬度转换后的欧氏坐标为B _i ，则矩阵H是它们的转换的关系。对公式B _i Hb _i=0用最小二乖法求转换矩阵H的初始值，再用Levenberg-Marquardt非线性求得更优的转换矩阵H，H需保持转换后坐标轴的正交性和比例的一致性。若给定多对匹配点可以解出线性的最小二乘解。进一步，需加入旋转、平移、等比例缩放约束再用Levenberg-Marquardt算法进行非线性优化，即得到最终H。 

3）在3D模型中设置变换后的视点，得到相应视点下三维图形渲染显示后的二维视图（如图2~7的b）所示），即对应二维图像的已经配准的虚像。 

本发明的有益效果是，本发明的一种3D模型与二维图像配准算法，它利用多视几何理论恢复3D模型在相对坐标中的视点，只需预先测量得到摄像机位置，将二维真实可见光图像与3D模型投影虚拟图进行自动配准。它不需要进行摄像机的预标定，也不需要用户手动设置初始的模型对应于各幅二维图像的摄像机视角，而是通过有多视几何关系的多幅图像，用估计加优化的方法求解摄像机内参数和外参数，再变换三维恢复坐标下的摄像机视点到3D模型坐标的视点位置及视角。配准后可以进行虚实图像融合，得到更佳的显示效果。一旦焦距或是摄像机内参、旋转方向发生变化，重新将本算法执行一次即可。   

附图说明

图1是一个示例中摄像机位置及世界坐标系下的特征点位置。图中：重建的摄像机位置(*)及世界坐标系下的特征点(+)位置。 

图2是一个示例中第1幅虚实图像对的视点配准结果。 

图3是一个示例中第2幅虚实图像对的视点配准结果。 

图4是一个示例中第3幅虚实图像对的视点配准结果。 

图5是一个示例中第4幅虚实图像对的视点配准结果。 

图6是一个示例中第5幅虚实图像对的视点配准结果。 

图7是一个示例中第6幅虚实图像对的视点配准结果。 

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。 

本发明的算法的使用方式是：在用多台PTZ摄像机进行场景监视的环境中，预先生成监视环境的3D模型，再测量出各个摄像机的经纬度坐标。通过本发明的算法，找到摄像机在3D模型中对应的视点及相关摄像机参数，然后按上述参数进行3D模型渲染，得到配准视图。图片2至7给出了一个监视场景的配准的示例。 

本发明的算法首先要假定采用哪种摄像机模型，本发明中认为摄像机符合小孔成像的标准摄像机模型(Pinhole model)，采用以下透视投影公式描述： 

X ≌ K[R | -Rt _c]b

式中，x，b均为齐次坐标，分别为图像坐标和世界坐标。世界坐标系与摄像机坐标系之间的关系用一个旋转矩阵R和一个平移向量t _c 来描述。a _r 和s分别是摄像机的像素纵横比（Aspect ratio）和倾斜度(Skew）；f为等效焦距（像素为单位）；(u₀,v₀)为图像平面坐标中心点，即光轴与图像成像平面的交点。为了进一步简化模型，假设：

(1)     像素纵横是同一个尺度，a _r =1。

(2)     图像坐标X轴与Y轴相互垂直，s=0。 

(3)     图像没有畸变。 

(4)     当焦距变化时，假设光轴和图像平面坐标中心点不变，即(u ₀, v ₀)保持不变。 

由于现在摄像机的制作工艺已经很成熟，这些假设一般情况下均成立。 

在上述假设条件下，本发明的具体实施步骤如下： 

步骤1）通过有多视几何关系的多幅图像求解摄像机内参数和外参数

步骤1-1）提取SIFT特征点

SIFT特征是图像的局部描述量，对尺度缩放、平移、旋转、噪声和亮度变化等具有良好的不变性，对仿射变换、视觉变化也保持一定的稳定性。特征维数较高（128维）。

SIFT算法提取特征点主要分4个步骤: 

（1）首先利用高斯差分（DoG）以及影像金字塔找出尺度空间中的极值点。

（2）求亚像素精度的极值点。因为DoG算子会产生边缘响应，还要去除低对比度的关键点和不稳定的边缘响应点。 

（3）将剩余的极值点定义为特征点后，开启一窗口，利用窗口内的影像梯度和幅值定义出特征点的方向。 

（4）按方向进行旋转，生成SIFT特征向量。 

这个向量就是所求的SIFT特征描述符。SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，如归一化特征向量，则进一步减弱光照变化的影响。 

步骤1-2）特征匹配

这里采用最近邻与次近邻的距离比率作为判断准则。最近邻特征点指与样本特征点欧氏距离最短的特征点，次近邻特征点是具有只比最近邻距离大的欧氏距离（即距离第二小）的特征点。

设两幅图像中检测出的两个特征点向量分别为R _i =(r _i,1, r _i,2,…, r _i,128)和S _i=(s _i,1, s _i,2,… , s _i,128),是它们的欧氏距离为: 

                                                                                           (1)

设S _i ’是次近邻特征点，给定阈值T，当

                                                                                                              (2)

则认为R _i , S _i 匹配，是同一特征，其中，d(R _i ,S _i )为最近邻的特征距离，d(R _i ,S _i ’)为次近邻的特征距离。可以调整T改变检测出的特征点的个数，一般取T=0.8。本步骤采用的均为高清摄像机，图像分辨率达300万，所以特征点较多，可能排除的部分正确匹配对结果影响不大。

一般考虑匹配采用穷举法进行遍历搜索，但SIFT描述符的维数较大，而且一幅300万像素的图像可能有10000个特征，特征点很多，如果采取一一对应查找最近邻域点的时间复杂度为O(n²)，很费时，效率较低。为了提高匹配速度，本发明采用构造特征点K-D树的方法，可以将查找时间复杂度降低到O(logn)。 

其步骤是：为第一幅图像特征描述符分别建立K-D树。虽然建立K-D树比较费时，但在特征点比较匹配时会大大提高效率。依次查找图像每个特征点R _i ，在待第一幅图像的特征描述的K-D树上的最近邻S _i 和次近邻特征点S _i ’，若满足公式(2)则R _i , S _i 就是一对匹配点对。 

步骤1-3）RANSAC框架下8点法估计基本矩阵F

双目视觉***的外极几何***中，对于第一个摄像机像面上任意一个图像点x _i1 ，在第二个摄像机像面上对应的像点x _i2 一定位于该摄像机上的极线l₂上。这就是所谓的“外极”约束。其中，b_i为一个三维空间中的点，O₁和O₂是两个摄像机的光心。e, e’是两个光心的直线与图像平面的交点，称为极点。

外极约束可以用基本矩阵(Fundamental matrix)F公式表示： 

x _i2 ^T  Fx _i1 =0                                                                                                      (3)

其中，x _i1 和x _i2 是图像点的齐次坐标。F因为包含了一个尺度因子，所以是一个秩为2的3×3矩阵，有8个参数。

步骤1-3-1）8点法估计基本矩阵F

根据前面的特征点匹配，得到不同图像帧之间的匹配点对后就可以计算两帧之间的F。令F=(f_ij)，设匹配点对坐标是x _i1 =(u_i1,v_i1,1)^T , x _i2  =(u_i2,v_i2,1)^T，根据公式(3)，可以写成如下形式：

           (4)

只要给定大于8对的点对x _i1 , x _i2 ，基本矩阵F即可由方程(4)解出。通常，可以将f ₁₁~f ₃₃做为参数，而前面的***组成矩阵A。

线性方程组(4)，给定8个以上的点对应，A一般会是一个秩9的矩阵。这样的方程组只有零解。这时需要求解约束条件‖f‖=1下的最小二乘解，即求解下列问题： 

                                                                                                             (5)

为求该最小二乘问题的解A，将A奇异值分解为A=UDV ^T ，取V的最后一个列向量为F的值，即f=V₉，重新排列为3×3矩阵即为F。

我们需要找到一个矩阵去逼近F作为基本矩阵的估计以保证F秩为2，即求解如下最优化问题： 

                                                                                                               (6)

由上述方程得到的解是基本矩阵F的未用归一化坐标的最终估计。为了达到更高的精度，本步骤需对输入匹配点的图像坐标进行了归一化，得到后又解除了归一化。

本步骤用RANSAC估计F时，每一轮迭代中N=8。 

步骤1-3-2）用RANSAC估计F

RANSAC方法可以从有错误的匹配点对集合中估计基本矩阵F的鲁棒结果。其算法的主要步骤如下：

(1)     用1-2）SIFT特征匹配的方法得到初始的匹配点对（其中一般有一些误匹配结果）。

(2)     从第(1)步的匹配点对中，随机抽取8个点对，用1-3-1）基本矩阵F的8点算法计算得到基本矩阵F。 

(3)     用第(2)步结果F计算内点集合S(F)。 

(4)     如果内点集合S(F)中集合残差值err小于原先的最优内点集合S(F)_best，则采用当前内点集合S(F)及对应的F，即将S(F)_best设为S(F)，对应的最优F_best设为F。 

(5)     当迭代次数小于2048次，回到第(2)步迭代。 

(6)     当迭代次数达到2048次，终止迭代估计。F_best是RANSAC迭代得到的最优结果。 

一对匹配点(x _i1 , x _i2 )是否为内点的判断标准是 x _i2 到由x _i1 和F计算出的对应极线l ₂ ^’的距离,记为d(x _i2 , l ₂ ^’)： 

                                                                       (7)

同理，记x _i1 到由x _i2 和F计算出的对应极线l ₁ ^’的距离是d(x _i1 , l ₁ ^’)。

如果d(x _i2 , l ₂ ^’)+ d(x _i1 , l ₁ ^’) <4，则认为(x _i1 , x _i2 )是匹配的内点。 

内点集合S(F)中集合残差值的定义如下： 

                                                                         (8)

其中，x∈S(F)即x是所有的内点匹配对，err是计算所有内点匹配到计算出的极线的距离之和。

步骤1-3-3）将所有内点作为输入再求F更优解

得到上步RANSAC框架下最优的基本矩阵后，为了得到更优的结果，再采用1-3-1）（基本矩阵F的8点算法），将所有匹配的内点（内点个数≥8）做为输入，得到所有内点匹配意义下外极几何约束下的F优化解。

步骤1-3-4）用LM算法优化得到F最终结果

得到上述内点匹配意义下外极几何约束下的F优化解后，再以公式(8)为目标函数，用Levenberg-Marquardt进行非线性优化，得到更优的基本矩阵F。这就是RANSAC框架下8点法估计基本矩阵F的最终结果。

步骤1-4）估计本质矩阵E

如果已知内参数矩阵置，像平面使用归一化坐标x’ _i1 , x’ _i2 ，则称归一化坐标下的基本矩阵为本质矩阵。记为本质矩阵E，设双目视觉对应的摄像机内分别是K ₁，K ₂， 

x _i2 ^T K ₂ ^-T EK ₁ ^-1 x _i1 =0                                                                                   (9)

F=K ₂ ^-T E K ₁ ^-1  =K ₂ ^-T [t]_× R K ₁ ^-1                                                                    (10)

E=[t]_× R                                                                                                 (11)

基本矩阵F包含了两个摄像机的内参数矩阵K和摄像机本质矩阵E之间的相对关系。本质矩阵描述了两幅规范化图像间的极几何，它与基本矩阵的不同之处是它仅与摄像机的外参数有关，即t, R。

步骤1-4-1）估计本质矩阵E

同理，本质矩阵也是秩2的，但它仅有5个自由度，所以求解此E矩阵有5点法、6点法、7点法至8点法。一般由8个图像点对应，从特征点对应用公式(9)可建立关于本质矩阵的线性方程线性求解。由于公式的齐次性，只能获得相差一个非零尺度因子意义下的E。

对于实际图像，常用的方法是将E的最小奇异值置为零，用另外两个奇异值的均值作为两个相等的奇异值，这样的矩阵作E最小二乘意义下的最佳近似。即分解E=Udiag(σ ₁,σ ₂,σ ₃)V ^T ，E的近似矩阵是： 

E=Udiag((σ ₁+σ ₂)/2, (σ ₁+σ ₂)/2,0)V ^T

步骤1-4-2）由F估计本质矩阵E

由公式(10)可以推出基本矩阵F与本质矩阵的关系：

E=K ₂ ^T FK ₁                                                                                       (12)

假设图片的大小为宽度w，高度h，按照本发明设定的***投影模型，则摄像机的内参数初始值为：

                                                                              (13)

本发明采用这种算法在只有2张图片时估计相机的本质矩阵和内参数，虽然这种模型过于理想，与实际情况有所差距，但是经过基于稀疏Levenberg-Marquardt的集束优化，优化结果可以弥补理想模型的误差。

步骤1-4-3）由本质矩阵E估计出R和t

如果第一个摄像机矩阵规定为M ₁ =(I,0)，令E=Udiag(σ,σ,0 ₃)V ^T 是它的一个奇异值分解，则第二个摄像机投影矩阵M ₂ =（R，t）的解可能为四组：

其中， , 

四个解的几何解释中只有一个是合理的，即这个重构点同时位于两个摄像机前方的一对摄像机是合理的。

步骤1-4-4）用三角法解出特征点的世界坐标

由上述方法求解出每一个摄像机的本质矩阵E_i后，结合摄像机内参矩阵K_i，得到摄像机投影矩阵M_i，再解出每一个特征点的世界坐标。

步骤1-5）基于稀疏Levenberg-Marquardt的集束优化摄像机外参数以及内参数、所有特征点的世界坐标

设x _ij 是第i个三维点投影到第j个图像上（由第j个视点看到）的坐标，a _j 是摄像机的第j个内外参数，b _i 是第i个世界坐标系下三维点的坐标。集束优化的目标函数是：

                                                   (14)

Q(a _j , b _i )是第i个点投影到第j幅图像上的投影函数，可以由这个投影函数估计x _ij ；d(Q(a _j , b _i ),x _ij )是由投影函数估计的二维图像中坐标和真实的二维图像坐标x _ij 的欧氏距离；v _ij 表示第i个坐标点是否可以由第j个摄像机看到，即是否在第j幅图像中可见。

如果摄像机内外参数有k个，而三维坐标点需3个坐标，则此集束优化的需求解的参数为mk+3n个。本发明中摄像机的内外参数为k=9个，分别是：摄像机在世界坐标系中的位置3个，表示摄像机平移；旋转分量3个，表示摄像机的朝向；缩放尺度1个；径向畸变参数2个。 

设求解参数的矢量为P，P=(a ₁ ^T , a ₂ ^T ,…a _m ^T , b ₁ ^T , b ₂ ^T ,…b _m ^T ,)^T, 矢量X是图像上特征点的坐标，记为P=(x ₁₁ ^T , x ₁₂ ^T ,…x _1m ^T , x ₂₁ ^T , x ₂₂ ^T ,…x _2m ^T ,…x _n1 ^T , x _n2 ^T ,…x _nm ^T  ) ^T ，对图像上特征点坐标量测值X，有相应的通过投影函数计算得到的估计值， ，

并且相应的误差为：

那么最优化问题变成了求解的最小值, , P ^* 为所求的参数结果。

设雅可比矩阵，将J分块为[A | B]，。由此，LM算法的更新向量δ _p 可以写成(δ _a ^T ,δ _b ^T ) ^T ，则标准的迭代求解步骤J ^T Jδ _p =J ^T ε变成： 

                                                               (15)

实际的特征点量测位置可以组织成一维向量X，

X=(x ₁₁ ^T , x ₁₂ ^T ,…x _1m ^T , x ₂₁ ^T , x ₂₂ ^T ,…x _2m ^T ,…x _n1 ^T , x _n2 ^T ,…x _nm ^T  ) ^T

则参数向量P为P=(a ₁ ^T , a ₂ ^T ,…a _m ^T , b ₁ ^T , b ₂ ^T ,…b _m ^T ,)^T

更新向量δ _p 为δ _p =(δ _a1 ^T ,δ _a2 ^T ,…δ _am ^T , δ _b1 ^T ,δ _b2 ^T ,…δ _bm ^T ,)^T。

定义，则雅可比矩阵J可以写成： 

近似计算Hessian矩阵J ^T J分块为：

         (16)

其中,  ，U和V是对角块矩阵，W是稀疏矩阵。

因此公式(15) (16)引入λ后变成： 

                                                                (17)

实际计算时先U, V, W，再用1+λ乘U和V的对角元素。

计算逆矩阵V ^-1并定义Y=W V ^-1。求解(U-WY ^T)δ _a =ε _a -Yε _b ，求出δ _a ，因为a中的参数个数远少于b中的参数个数。求解δ _b =V ^-1(ε _b -W ^T δ _a )，求出δ _b 。通过增加更新矢量(δ _a ^T ,δ _b ^T ) ^T 更新参数矢量P，并计算新的误差矢量ε。如果新的误差矢量ε小于老的误差，则接受这个新的参数值，把λ减少至1/10，并且返回U, V, W的计算继续，否则终止计算。如果新的矢量ε大于老的误差，把λ增大10倍，并且继续计算实际计算时先U, V, W，再反复进行迭代计算。 

步骤1-6）RANSAC框架下单摄像机投影矩阵估计

上述的步骤已经可以从双目视觉***中求得两个摄像机的相对位置R, t _c,但是对于多目视觉中，如果有更多的摄像机指向了相同的场景，应该充分利用更多的立体视觉信息，得到更精确的视点相对位置。

第一步，通过提取SIFT特征点，得到多视图的特征点对应关系，即 

(x _i1 , x _i2,  x _i3 ,… x _im )。第二步，对于第3个摄像机，再根据(x _i1 , x _i2 )双目视觉得到的三维世界坐标b_i与x _i3 ，利用RANSAC框架求解摄像机投影矩阵M₃，再解得第3个摄像机的求相对位置R ₃, t _c3 。这样就得到了第3个摄像机相对位置的初始值。第三步，将代入基于稀疏Levenberg-Marquardt的集束优化方法即可得到第3个摄像机相对位置的优化结果。第四步，对于加入更多的指向同一场景的摄像机，可以采用继续采用上述几步计算相应摄像机的相对位置。

下面分5个步骤描述本步骤1-6）所应用的方法： 

步骤1-6-1）6点法求得摄像机矩阵M

在经典立体视觉中标定方法是根据一些空间点在世界坐标系下的坐标b与其图像坐标x之间的对应，即x～Mb，建立约束方程（14），再摄像机矩阵M。这里世界坐标的空间点b _i由(x _i1 , x _i2 )已经由上述的双目视觉计算初始值并迭代优化得到。

在基于多幅图像的场景重建问题中，投影矩阵的估计是非常重要的环节。X与b关系表示为如下公式： 

x×Mb=0                                                                                        (18)

不妨设已知n组对应点(b_i, x _i), x _i=(u_i,v_i,1)^T, ，M _i ^T 表示矩阵M的第i行向量。代入公式(18)，消去尺度因子常数后，可得到下述公式：

A _i  p=0                                                                                                           (19)

其中，, 

给定n个点（6个或更多）对应，得到n个形如这样的矩阵A _i，再将这n个矩阵A_i组合起来得到一个3n×12 的矩阵A=(A ₁ ^T , A ₂ ^T , …A _n ^T )，对A作奇异值分解A=UDV ^T ，则V的最后一个列向量p=V ₁₂是方程Ap=0的最小二乘解，再将p写成3×4矩阵的形式就是摄像机矩阵M。

步骤1-6-2）用RANSAC估计M

本发明采用RANSAC估计世界坐标系下的坐标b与其图像坐标x之间的摄像机矩阵M，可以得到更好的结果。

(1)     用1-1）步SIFT特征匹配的方法得到第j幅图片与前j-1幅图片的初始的匹配点对（同样，其中一般有一些误匹配结果）。 

(2)     从第(1)步的匹配点对中，随机抽取6个点对（x_ij,b_i），用求摄像机矩阵M的6点法计算得到第i幅图片对应的摄像机矩阵M_j。 

(3)     用第(2)步结果M_j计算内点集合S(M_j)。 

(4)     如果内点集合S(M_j)中集合中内点的个数大于原先的最优内点集合S(M_j)_best，则采用当前内点集合S(M_j)及对应的M_j，即将S(M_j)_best设为S(M_j)，对应的最优M_jbest设为M_j。 

(5)     当迭代次数小于4096次，回到第(2)步迭代。 

(6)     当迭代次数达到4096次，终止迭代估计。M_jbest是第i个摄像机RANSAC迭代得到的最优结果。 

判断（x _ij, b_i）是否是内点的标准是b_i经过投影到图像上与像素坐标x_ij的欧氏距离，表示为d(x _ij ，M _j b _i )=‖x _ij- M _j b _i ‖, 如果二者之间的距离小于4个像素，则（x _ij, b_i）是内点。 

步骤1-6-3）将所有内点作为输入再求M_j 更优解

得到上步RANSAC框架下最优的摄像机矩阵M_j后，为了得到更优的结果，再采用6点法，将所有匹配的内点（内点个数>6）做为输入，得到所有内点匹配下的M_j最小二乘优化解。

步骤1-6-4）用LM算法优化得到M最终结果

得到上述所有内点匹配下的M_j最小二乘优化解后，再以最大化内点集合中内点个数为目标函数，用Levenberg-Marquardt进行非线性优化，得到更优的摄像机矩阵M_j。这就是RANSAC框架下单摄像机投影矩阵估计M_j的最终结果。

步骤1-6-5）根据投影矩阵M _j 求相对位置

本发明的算法通过投影矩阵M采用RQ分解求解摄像机内参数、外参数。

将M表示成M=[m ₁, m ₂, m ₃, m ₄]，m _i表示M的第i列。令摄像机内参数矩阵为K，外参数矩阵为(R, t _c )，由于此时的摄像机矩阵是欧氏的，所以可写成： 

M=κK[R | -Rt _c] 

首先，按下列公式计算出t_c。

上式中的M ₁₂₃表示矩阵[m ₁, m ₂, m ₃]，m _i表示M的第i列，其它M ₂₃₄等为相同的意义。 

再利用RQ分解，[K R]=RQ((M ₁₂₃))，K是对角元素均为正数的上三角矩阵，R为正交旋转矩阵，M ₁₂₃₌ KR。 

这样，就得到了摄像机外参数(R, t _c )。 

步骤1-7）反复执行直到所有的图片处理完毕

反复执行步骤1-5）~1-6）直到每一个摄像机所对应的图片处理完毕，得到摄像机的内外参数后，优化求出特征点的三维世界坐标后，其一个示例结果如图2~7所示。

步骤2）图像摄像机视点到3D模型的视点配准

根据预先好的测量摄像机的经纬度找到多个视点的相对位置到3D模型的相对变换关系（旋转、平移、缩放），可以将视点位置变换到3D模型中。我们采用的方法是经纬度转换法，它分为2个步骤：

步骤2-1) 将测量的摄像机经纬度坐标转换成模型坐标系中的坐标，采用高斯-克吕格投影（Gauss-Kruger projection）。

步骤2-2) 对公式(20)用最小二乖法求转换矩阵H的初始值，再用Levenberg-Marquardt非线性求得更优的转换矩阵H，H保持转换后坐标轴的正交性和比例的一致性，它仅由平移，旋转，缩放的7个参数确定。设摄像机位置的归一化的相对三维坐标为b _i ，经纬度转换后的欧氏坐标为B _i ，则它们与转换矩阵H的关系是： 

B _i Hb _i=0                                                    (20)

若给定多对匹配点可以解出线性的最小二乘解。进一步，需加入旋转、平移、等比例缩放约束再用Levenberg-Marquardt算法进行非线性优化，即得到最终H。

步骤3）在3D模型中设置变换后的视点并进行渲染

在3D模型中设置变换后的视点，得到相应视点下三维图形渲染显示后的二维视图（如图2~7中的b）图所示），即对应二维图像的已经配准的虚像。这时三维图形渲染显示后的二维视图与图像已经配准，算法执行完毕。

经过本步骤的集束优化后，得到的所有摄像机内外参数和世界坐标系下的特征点位置，就是最后的结果。 

  

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。

  

Claims (5)

1.一种3D模型与二维图像配准算法，其特征是具体步骤如下：

1）通过有多视几何关系的多幅图像求解摄像机内参数和外参数

1-1)提取SIFT特征描述（它具有尺度和旋转不变性）;

1-2)再用KD-树进行特征点匹配;

1-3) 采用RANSAC框架通过8点法估计基础矩阵来去掉噪声匹配点对;

1-4) 选择共有特征点最多的两幅图像，再通过上步得到的基本矩阵F，按估计的摄像机参数算出本质矩阵E，估计一对摄像机之间的相对位置，得到摄像机位置旋转和平移参数;

1-5) 初始化摄像机内参数，将上述外参数以及内参数、所有特征点的投影坐标和世界坐标作为初始值，用稀疏Levenberg-Marquardt方法进行集束优化;

1-6) 再选择另一幅图片，根据已算出的空间特征点及特征点跟踪关系用RANSAC方法求解此摄像机投影矩阵，得到摄像机内外参数初始值，将其加入到优化框架下，转步骤1-5）继续优化；

1-7) 直到所有的图片处理完毕，得到所有摄像机的外参数（多幅图片视点的旋转平移的相对位置）

2）变换三维重建坐标下的摄像机视点到3D模型坐标的视点

根据预先好的测量摄像机的经纬度坐标找到多个视点的相对位置到3D模型的相对变换关系（旋转、平移、缩放），可以将视点位置变换到3D模型中

3）在3D模型中设置变换后的视点，得到相应视点下三维图形渲染显示后的二维视图，即对应二维图像的已经配准的虚像。

2.如权利要求1所述的配准算法，其特征在于：所监视场景的3D模型预先生成，所有摄像机的位置也于运行算法前测量得到。

3.如权利要求1所述的配准算法，其特征在于：在步骤1-6）中，采用通过投影矩阵M采用RQ分解求解摄像机内参数、外参数；即：将M表示成m _i表示，M=[m ₁, m ₂, m ₃, m ₄]的第i列；令摄像机内参数矩阵为K，外参数矩阵为(R, t _c )，所以此时的摄像机矩阵M可写成：M=κK[R | -Rt _c] 按下列公式计算出平移量t_c：

上式中的M ₁₂₃表示矩阵[m ₁, m ₂, m ₃]，m _i表示M的第i列，其它M ₂₃₄等为相同的意义；再利用RQ分解，[K R]=RQ((M ₁₂₃))，K是对角元素均为正数的上三角矩阵，R为正交旋转矩阵，M ₁₂₃₌ KR；这样，就得到了摄像机外参数(R, t _c )。

4.如权利要求1所述的配准算法，其特征在于：将测量的摄像机经纬度坐标转换成模型坐标系中的坐标时，采用高斯-克吕格投影（Gauss-Kruger projection）。

5.如权利要求1所述的配准算法，其特征在于：用最小二乖法求B _i Hb _i=0中转换矩阵H的初始值，其中，b _i 为摄像机位置的归一化的相对三维坐标，B _i 为经纬度转换后的模型坐标系中的坐标，；再用Levenberg-Marquardt非线性求得更优的转换矩阵H，在非线性迭代中H保持转换后坐标轴的正交性和比例的一致性。