第3章
CHAPTER 3


基于平面标志的空间注册








在增强现实中，标志指现实场景中存在的易于识别和定位的物体，其主要作用是辅助相机或者物体的定位跟踪，以较为简单、稳定的方式实现空间注册。在现有增强现实系统中，平面标志是应用最为广泛的一种标志技术，这主要是因为平面的制作、表示和计算都较为简单。平面标志既可以由人工设计制作，也可以是场景中自然存在的平面物体(如书页、海报等)。本章将以人工平面标志为基础，实现一个最为基本的增强现实系统。如无特殊说明，本章所述平面标志皆为人工平面标志，对自然标志物的处理方法将在第5章进行介绍。

3．1系统概述

如图3．1所示为本章将要介绍的简易增强现实系统。其中图3．1(a)是从视频捕获的输入帧，其中包含一个平面标志； 图3．1(b)是在图3．1(a)的基础上绘制了虚拟物体之后的增强现实效果，注意虚拟物体看起来就像被放在了平面标志上一样； 图3．1(c)和图3．1(d)是视频中另外两个时刻的效果图，可以看到，当摄像机运动或者平面标志移动的时候，虚拟物体都会随着平面标志一起运动，看起来像是真的固定在平面标志上的物体一样。这就是本书第1、2章讲到的虚实空间注册所要追求的效果。在这里，即是将虚拟物体所在的局部空间与根据平面标志定义的现实空间进行了注册。



视频演示




图3．1基于平面标志的增强现实系统


为了实现上述虚实融合效果，就需要知道绘制虚拟物体所需的参数。根据本书第1、2章的内容和计算机图形学的相关知识可知，所需的参数主要包含两类： 一是虚拟物体相对于相机的位置和姿态参数，可以表示为虚拟物体从其局部坐标系变换到相机坐标系的旋转矩阵R和平移向量t，根据这个参数可以将虚拟物体嵌入相机坐标系中的对应位置； 二是相机的内参矩阵K，用于确定相机坐标系中一个3D空间点在图像上的对应像素点位置。在OpenGL绘制程序中，R和t将用于决定模型视图变换，而K将决定投影变换。关于参数R、t、K与OpenGL中各变换矩阵的转换关系，将在3．4节进行介绍。

对于定焦相机而言，相机的内参矩阵K是固定不变的，并通过2．4．1节介绍的相机标定方法获得。因此，关键是计算外部参数R和t。由于R和t跟摄像机位置和虚拟物体在场景中的摆放位姿相关，因此需要根据输入视频流进行实时计算。这也是所有增强现实系统需要解决的首要问题。在一般情况下，这需要持续稳定地跟踪摄像机的3D位姿参数，同时还要知道虚拟物体摆放处的部分现实场景的3D几何，二者皆是非常具有挑战性的问题。

在本章中，将主要介绍基于平面标志的R和t的简化计算方法。为此，首先假设算法的目标是把虚拟物体放置到平面标志上，如图3．1所示的效果。这样，虚拟物体在相机坐标系中的位姿，实际上也就是平面标志在相机坐标系中的位姿； 或者更一般地，二者之间只差一个不随时间变化的恒定变换，可以在系统初始化的时候设定，相当于指定虚拟物体和平面标志之间的相对位姿关系。只要知道了平面标志相对于摄像机的位姿参数，就可以容易地得到虚拟物体的位姿参数R和t。

3．2平面位姿估计

计算平面标志在相机坐标系中的位姿的基本方法是，首先获得一组平面标志上的3D空间点与图像上相应像素点的对应关系，即3D2D点对，再以此为约束求解位姿参数。实际上，这也是求解一般3D物体位姿的基本方法。对于一般3D物体而言，根据一组3D2D点对求解其姿态参数可以采用PnP(PerspectivenPoint)方法，相关内容将在4．2节中进行介绍。对于平面物体的位姿而言，虽然也可以采用PnP方法进行计算，但是由于平面的3D位姿存在特殊表示，即单应性矩阵，因此在本节中，先介绍一种专门针对平面的位姿计算方法。这里假设已知一组3D2D点对，如何获取这些点对的方法将在3．3节中介绍。

3．2．1单应性矩阵

回顾2．3节介绍的针孔相机模型，3D空间点X投影到图像上相应像素点x的过程可以表示为
x～~K(Rt)X～(3．1)
式中，变量上方的符号“~”表示相应的齐次坐标，连接两组代数式的符号“~”则表示其中一边乘以一个恰当的数字会使两边相等。对于平面物体而言，可以为其选取一个局部坐标系，使得该平面物体上任意一个3D空间点的Z坐标都等于0，因此有
x～~K(r1r2r3t)X
Y
0
1=K(r1r2t)X
Y
1=HX
Y
1(3．2)
其中
H~K(r1r2t)(3．3)
是一个3×3的可逆矩阵，r1、r2、r3是旋转矩阵R的列向量。

上式中的矩阵H不仅包含物体的旋转和平移信息，同时还与相机内部参数K有关，可以表示3D场景中的平面物体通过透视变换映射到图像平面的过程，称其为3D平面的单应性变换，相应的变换矩阵被称为单应性矩阵。单应性变换其实是一个2D的射影变换，其基本特点是保持直线性，即直线经单应性变换之后仍然为直线。一个单应性变换的典型例子是如图3．2所示的中心投影。π和π′是3D空间的两个平面，O是空间中一点，并且不在这两个平面上。这样，平面π上任意一点x通过与O的连线都与平面π′相交，且交点唯一，即x′。这样就建立了从平面π到π′的一个一一映射，且这个映射显然是可逆并且保持直线性的，因此是一个单应性变换，可以用单应性矩阵表示。



图3．2中心投影


实际上，上述中心投影模型与针孔相机的成像过程也是一致的，中心点O相当于相机的光心，若平面π′与光轴正交，则相当于相机的成像面。因此，平面π到平面π′的映射总是可以用一个单应性矩阵来表示。这也正是式(3．2)所表示的情况，在本章中将使用单应性矩阵来表示和计算平面标志的3D位姿。

根据中心投影模型和单应性变换的性质，可以得出其他一些可以用单应性变换来表示的情况。典型地，如果π是场景中的一个3D平面，则其在不同视点O和O′下所观测到的图像上的像素对应关系也可以表示为一个单应性变换。实际上，设平面π上的点X通过矩阵H1映射到第一张图像上的点x，即x～~H1X～，并有X～~H-11x～； 通过矩阵H2映射到第二张图像上的点x′，即x～′~H2X～，并有X～~H-12x～′，则经过简单的推导，有
x～′~H2H-11x～=Hx～(3．4)
在式(3．4)中，H=H2H-11。显然，矩阵H仍然是3×3的可逆矩阵，因此也是一个单应性矩阵。进一步地，可以得到，任何两个平面都可以通过一个3×3的矩阵建立一一映射。对于不同视频帧中观测到的场景中的同一个平面而言，可以用一个单应性变换来建立两帧之间的像素对应关系，这在图像匹配中有重要作用。

另一个典型的情况是，由旋转相机(相机光心静止，只有绕光心旋转的运动)拍摄到的不同图像之间的像素对应关系，也可以用一个单应性变换来表示。这种情况也可以用如图3．2所示的中心投影来解释，此时平面π和平面π′相当于相机在不同旋转角度下的对应成像面，且由于相机光心和场景3D点都没有发生变化，只有相机成像面在运动，因此也是符合中心投影模型的。注意，对于旋转相机的情况而言，即使所拍摄的3D场景不是平面，也可以用单应性变换来精确表示像素的对应关系。实际上，当视点固定时，基线长度为0，因此有t=0，坐标变换矩阵为
T=R0
0T1=r11r12r130
r21r22r230
r31r32r330
0001(3．5)
因此，透视投影矩阵为
P=K(I0)R0
0T1=K(R0)=(KR0)(3．6)
可见，此时的透视投影矩阵也退化为3×3的矩阵，也就是说，尽管场景中的物体形状各异，远近不同，但是，摄像机在不同姿态下拍摄的图像之间，采用3×3的单应性矩阵就可以表达两幅图像间的映射关系。在使用手机的全景拍照功能时会发现，当相机运动接近一个旋转相机时，全景图的拼接效果最好，其原因就在于，只有旋转相机才能在不恢复场景3D几何结构的情况下，基于图像的单应性变换将不同角度拍摄到的图像进行精确对齐。

3．2．2求解单应性矩阵

为对单应性矩阵进行求解，将单应性矩阵表示为如下形式： 
H=h11h12h13
h21h22h23
h31h32h33(3．7)
因此，求解单应性矩阵即是要求解式(3.7)中的9个矩阵参数，将其表示为一个参数向量h=(h11h12h13h21h22h23h31h32h33)T。注意，由于单应性变换前后都使用齐次坐标，因此，将单应性矩阵乘以一个非零常数，其所表示的变换是不变的。不失一般性，假定h是一个归一化向量，即‖h‖=1，那么，虽然一个单应性矩阵包含9个未知参数，但其自由度实际上为8。由于每个点对可以提供两个约束，所以最少需要4个平面到平面的点对才能唯一确定h。

如果给定一组平面到平面的点对集合，就可以采用直接线性变换(direct linear transform，DLT)来估计满足这组点对约束的最优单应性变换。记(xi,yi)和(x′i,y′i)是一个对应点对，则二者应满足以下约束条件：  
x′i=h11xi+h12yi+h13h31xi+h32yi+h33，
y′i=h21xi+h22yi+h23h31xi+h32yi+h33(3．8)
上述形式对待求解参数h是非线性的。不过，可以在式(3.8)两边同时乘以分母，并将式(3.8)直接变换为如下线性形式： 
h11xi+h12yi+h13-x′ixih31-x′iyih32-x′ih33=0
h21xi+h22yi+h23-y′ixih31-y′iyih32-y′ih33=0(3．9)
将所有特征点对通过上述形式进行排列，可以得到：  
x1y11000-x′1x1-x′1y1-x′1
000x1y11-y′1x1-y′1y1-y′1
x2y21000-x′2x2-x′2y2-x′2
000x2y21-y′2x2-y′2y2-y′2
h11
h12
h13
h21
h22
h23
h31
h32
h33=0
0
0
0
(3．10)
式(3．10)中的系数矩阵记为A，则方程具有Ah=0的形式，是一个齐次线性方程组，一般情况下可基于奇异值分解(singular value decomposition,SVD)进行求解。具体地，如果矩阵A的SVD形式为A=UDVT，且对角阵D的对角线元素按降序排列，则待求解结果是矩阵V的最后一列。由于V是正交矩阵，所以，这样得到的参数向量h也满足‖h‖=1的约束。

DLT方法的一个缺点是，在式(3．8)的两边同时乘以分母，相当于给每一个点对关联的方程乘以一个权重，将导致坐标值较大的点对结果有较大的影响，这显然是不合理的。为了获得更精确的结果，往往以DLT方法所获得的结果为初值，并进一步采用非线性优化的方法来最小化重投影误差。

在非线性估计方法中，最常采用的是非线性最小二乘法。假设{(xi,x′i)，i=1,2,…}是输入点对，f是从xi到x′i的变换，即x′i=f(xi; Θ)，其中Θ是待估计运动模型的参数向量。非线性最小二乘法需要对当前估计参数Θ进行迭代并找到更新的ΔΘ，以最小化目标函数： 
E(ΔΘ)=∑i‖f(xi; Θ+ΔΘ)-x′i‖2
≈∑i‖J(xi; Θ)ΔΘ-ri‖2
=ΔΘT［∑iJTJ］ΔΘ-2ΔΘT［∑iJTri］+∑i‖ri‖2
=ΔΘTAΔΘ-2ΔΘTb+c(3．11)
其中，J为f的雅可比矩阵，且
A=∑iJT(xi)J(xi)b=∑iJT(xi)ri(3．12)
一旦计算出A和b，即可用下式计算ΔΘ： 
AΔΘ=b(3．13)
并相应更新参数向量Θ为Θ+ΔΘ。

可见，采用非线性最小二乘法的关键是计算变换的雅可比矩阵J。对单应性矩阵而言，Θ=h，f为式(3．8)的形式，其雅可比矩阵为
J(xi)=f(xi)h=1Dixiyi1000-x′ixi-x′iyi-x′i
000xiyi1-y′ixi-y′iyi-y′i(3．14)
其中Di=h20xi+h21yi+h22。注意，式（3.14）中Di、x′i、y′i的计算都依赖h，可以基于h的当前值进行计算。

3．2．3从单应性矩阵分解旋转和平移

对于场景中的一个平面物体而言，在获得其相对于图像的单应性矩阵H之后，还可以进一步分解得到该平面物体在相机坐标系中的旋转和平移参数。假设相机的内参矩阵K已知，则根据式(3．3)可得
（r1r2t）~K-1H(3．15)
不妨设r1、r2、t分别为矩阵K-1H的3个列向量。又由于旋转矩阵R是正交矩阵，因此可以得到r3=r1×r2，这样便得到了R、t的初始估计。

不过，由于对应点误差和计算误差的存在，按照上述方法直接计算得到的r1、r2不一定是正交的单位向量，因此不能保证R是一个旋转矩阵。为此，首先对H乘以一个系数λ：
(r1r2t)=λK-1H(3．16)
其中λ的取值应使r1、r2尽量接近单位向量。记λ1=1/‖K-1h1‖，λ2=1/‖K-1h2‖，其中h1、h2分别为H的第一个和第二个列向量，则可以取λ=0．5(λ1+λ2)。注意这里再次利用了单应性变换的齐次性，对H而言，乘以λ并不会改变其所表示的变换。

为了最终能够获得一个正交的旋转矩阵，就必须以上述方法所得结果为初值，求解一个离该结果最近的正交矩阵作为最终的旋转矩阵： 
minR‖R-Q‖2Fs．t．RTR=I(3．17)
其中Q=(r1r2r1×r2)是根据式(3．15)计算得到的结果，‖·‖F是矩阵的Frobenius范数，即矩阵各元素的平方和再开方。上述问题的最优解可以通过对矩阵Q的SVD得到，记Q=USVT，其中S为对角阵，可以证明最优解为R*=UVT。

3．3平面标志的检测

在3．2节中，假定已知一组平面标志到图像上的对应点对，就可以求解出平面标志在相机坐标系中的3D位姿。为了求解单应性矩阵，最少需要4组点对。对于如图3．1所示的正方形或矩形平面标志而言，一般选择正方形或矩形的4个角点建立对应。4个角点在平面标志上的坐标可以根据平面标志的尺寸提前设定，并在系统运行过程中保持不变。因此，主要问题是获取平面标志的4个角点在图像上的对应像素坐标，这也是对平面标志进行检测的主要目的。

根据平面标志设计特点的不同，对其进行检测的方法也会有所不同。对于如图3.1所示的正方形或矩形平面标志而言，其检测方法一般包含两个关键步骤： 第一步是检测出图像中的候选四边形区域，并定位四边形的角点和边； 第二步是根据每个候选四边形区域内部的图像内容，识别该四边形是否为平面标志，以及平面标志的类型(对包含多种标志的情况而言)。下面以平面标志检测库ARUCO中采用的检测算法为例，进一步说明上述两个步骤的具体实现过程。如图3．3所示为ARUCO中的平面标志检测算法工作的主要过程。



图3．3平面标志的检测过程


获得视频输入帧以后，首先将其转换成灰度图像，然后进行边缘检测或阈值化处理，获得如图3．3(b)所示的边缘图像。边缘图像是一张二值图像，其中白色像素表示强边缘所在的位置。由于平面标志一般都被设计为高对比度的黑白图案，因此组成平面标志的区域边缘可以被较为稳定地检测到。相比于Canny边缘检测，阈值化处理一般速度更快，且运动模糊会更稳定，缺点是对边缘的定位精度相对较低。图中所示图像为用OpenCV中的局部自适应阈值函数adaptiveThreshold计算获得的结果。

获得边缘图像后，可以进一步提取候选的四边形，结果如图3．3(c)所示。这可通过两个步骤完成。第一步是用轮廓跟踪算法找出边缘图中的闭合轮廓，并将轮廓边缘的像素连接为一个多边形。轮廓跟踪可以采用SuzukiAbe算法，该算法可采用OpenCV中的findContours函数来实现。第二步，为了从闭合轮廓中找出四边形，可以对表示闭合轮廓的多边形进行简化，将近似为直线的轮廓点用一条线段表示。如果简化后的轮廓只包含四条边，则可以认为是一个候选的四边形区域。多边形简化可以采用DouglasPeucker算法，该算法可用OpenCV中的approxPoly函数实现。

最后，基于图像内容，剔除候选四边形中的非平面标志区域，并识别每个平面标志的类别ID。这主要通过图3．3(d)的二值编码来完成。对每个候选四边形来说，都应先将其4个顶点与图3．3(d)所示的平面标志的标准形状建立对应。注意，由于无法仅根据形状确定具体的对应关系，因此需要对所有的4种可能逐一进行尝试。对每一种可能而言，都应根据4组点对应，采用3．2．2节介绍的方法求得从图像到标准形状的单应性变换H，再根据H将候选四边形内部的图像内容校正到标准形状。如果某个候选四边形是一个平面标志，且点对应正确，则经过H变换后的结果将与相应平面标志的标准模板对齐。为了度量对齐程度，应将经H变换后的图像区域转换成图3．3(d)所示的二值编码。在系统初始化时，对系统中所有可能的平面标志，采用同样的方式计算其图像的二值编码并建立索引表。为了确定一个候选四边形区域是否为平面标志，可以将其二值编码在索引表中进行检索，找出与之最相似的一个平面标志。如果相似度足够高，则认为该候选是一个平面标志，并赋予相应平面标志的类别ID。

可见，平面标志的特殊形状和颜色极大地简化了对其进行检测的过程，只需用基本的图像处理技术即可获得较为稳定的检测结果。在视频中，物体的运动一般是连续的，在视频的相邻帧之间不会有太大的差异。利用这个特点，还可以基于上一帧的平面标志位姿参数，跟踪出当前帧的平面标志位姿参数。相比于检测，跟踪只是进行局部搜索，因此一般来说可以更快、更精确。不过，当物体或者相机运动较快的时候，跟踪也容易因视频的连续性假设不满足而失败。因此，实际系统中往往需要检测和跟踪相互配合才能获得较好的效果。

3．4虚拟物体绘制

为了获得最终的绘制输出，首先需要将相机内参矩阵K和物体位姿参数R和t转换成图形，并绘制引擎的相应参数。这里以OpenGL为例，说明相应的转换方法。

在OpenGL中，相机内参矩阵K的作用等价于投影变换矩阵。回顾2．3．1节，如果不考虑成像面的倾斜角，则内参矩阵K可以写为如下形式： 
K=αx0x0
0αyy0
001(3．18)
不过，上述形式的内参矩阵K与OpenGL要求的投影变换矩阵是不一致的，无法直接设置。在OpenGL中，投影变换矩阵P是一个4×4的矩阵，除了表示3D到2D的投影关系外，还包含对深度值的变换。从K到P的变换可以采用如下形式：  
P=2αxw01-2x0w0
0-2αyh1-2y0h0
00-Z1+Z0Z1-Z0-2Z1Z0Z1-Z0
00-10(3．19)
其中，w和h分别为图像的宽和高，Z0和Z1分别为OpenGL中近裁剪面和远裁剪面所对应的深度值(具体可参考gluPerspective函数)。

根据R和t计算OpenGL的模型视图变换较为简单。但需要注意的是，由于在图像坐标系一般假设y轴向下，x轴向右，导致引出的z轴向里(远离视点)； 而在OpenGL中，图像坐标系的y轴向上，z轴向外，因此，在图像坐标系中计算R和t转换到OpenGL坐标系时，需要将y坐标和z坐标反向，相应的模型视图变换为
T=1000
0-100
00-10
0001Rt
01(3．20)

采用上述公式获得投影矩阵P和模型视图变换矩阵T之后，就可以采用glLoadMatrix函数设置变换矩阵到绘制管线，再进一步完成绘制。为了将虚拟物体直接绘制到图像中，可以将背景图像作为纹理贴图先进行绘制。

最后，总结一下本章所述的简易增强现实系统。其基本流程是，对每一帧输入图像，首先用3．3节的方法检测平面标志，定位出平面标志的4个角点； 其次用3．2节的方法，通过点对估计平面标志的3D位姿； 最后用3．4节的方法绘制虚拟物体。

上述过程主要解决的是虚拟物体在现实空间的注册问题，其中使用了特殊设计的平面标志作为辅助，并将虚拟物体的位姿与平面标志的局部坐标系相绑定。平面标志的优点是其特征点非常明显，易于检测。本书接下来的几章将讨论更一般的空间注册方法，其3D特征点坐标的获取各具特色，第4章介绍由各类传感器通过实时测量获取3D坐标的方法，第5章通过对自然图像上特征点进行检测和匹配来获取3D点坐标，第6章则采用视觉方法通过对视频图像序列上特征点的匹配来重构特征点的3D坐标。

扩展阅读

增强现实系统开发涉及图像视频处理、计算机视觉、计算机图形学等方面的专业知识，同时为了满足实时性，对代码的计算效率有较高的要求。实际项目中往往都基于已有的增强现实SDK实现底层核心功能。ARToolKit是一个较为早期的增强现实SDK，最早发布于1999年，具有开源、免费的优势，可以支持Windows、Linux和OS X，但不支持移动平台。近年来，随着手机、平板电脑等移动设备的普及，增强现实的主流应用逐渐转移到了移动计算平台。目前常用的增强现实SDK有高通公司的Vuforia、谷歌公司的ARCore、苹果公司的ARKit，以及视辰信息科技公司的EasyAR等。这些工具包都支持增强现实的基本核心功能，如相机标定、平面标志物的识别跟踪、基于SLAM技术的相机跟踪与场景3D重建，以及图形绘制等。关于这些工具包的更多信息都可以通过搜索引擎找到，感兴趣的读者可以自行查阅。

习题

1． 已知X=(X,Y,Z)T是物体3D模型上的一个点，物体在相机坐标系中的旋转和平移参数分别为R、t，相机内参矩阵为K。求X在图像上的投影点坐标x=(x,y)T。

2． 如果单应性矩阵H表示某3D平面到图像的变换，则H主要包含哪些信息？

3． π是一个3D平面，I1和I2分别是平面π在不同视角下观测到的图像，证明： 平面π在I1和I2中的投影像素坐标可以通过一个单应性矩阵进行精确对应，即存在单应性矩阵H，使得X∈π，有x～1~Hx～2，其中x～1和x～2分别是X在I1和I2中投影像素的齐次坐标。

4． 推导单应性矩阵的雅可比矩阵J，即式(3．14)。

5． 理解从K到OpenGL投影变换矩阵P的变换公式，即式(3．19)。

6． 根据本章所述基本原理，实现一个如图3．1所示的简易增强现实系统。