CN104504691A - 基于低秩纹理的摄像机位置和姿态测量方法 - Google Patents

Abstract

基于低秩纹理的摄像机位置和姿态测量方法，属于计算机视觉测量技术领域。其特征是：将室内场景中的一些低秩纹理如天花板上的方格纹理，用于摄像机位置和姿态的测量。通过摄像机拍摄低秩纹理的图像，用欧拉角表示摄像机拍摄时的投影矩阵，并通过求解关于低秩纹理成像模型的优化问题求得表示摄像机姿态的欧拉角；再进一步结合场景中纹理的几何关系和摄像机成像原理解出摄像机在场景中的位置。本发明的效果和益处是，能够利用室内场景中的低秩纹理特征实现摄像机姿态和空间位置6个自由度的测量，且在测量过程中，摄像机可以转动，从而使视觉测量的应用范围更广。

Description

基于低秩纹理的摄像机位置和姿态测量方法

技术领域

本发明属于计算机视觉测量技术领域，涉及一种基于低秩纹理的摄像机空间姿态和位置测量方法。

背景技术

室内测量技术是一项应用广泛的技术。随着近年来CCD技术和图像处理技术的不断提高，视觉测量技术以其适用范围广、精度高的特点，受到了研究者越来越多的关注。视觉定位技术采用一个或多个摄像机从场景中获取含有特征的图像，通过不同方法对图像特征进行表达和提取，最终利用这些特征提供的信息求解摄像机在空间中的位置和姿态。图像的特征分为局部特征和全局特征。局部特征主要包含角点、直线等。目前大多数视觉定位方法基于图像局部特征，从图像中提取点或直线等以实现摄像机定位。这类方法的测量精度依赖于局部特征提取的精度，极易受到噪声干扰，且在一些场合中精度往往达不到实际应用的要求。全局特征主要是图像的颜色、纹理等。纹理在室内的天花板、地板或者墙壁等平面上比较容易找到，常具有重复性或对称性，可以利用纹理所具有的这些特征进行摄像机相关测量。但目前使用纹理特征进行摄像机位置和姿态测量的方法较少。

另外，目前各方法多只将摄像机放置在某一个平面内，固定安装摄像机的拍摄角度，摄像机拍摄时不能转动。除此之外，大多数方法仅求解了摄像机在其所放置平面内的二维位置坐标，未求解摄像机在空间中的三维坐标。

发明内容

本发明的目的是提供一种基于低秩纹理的摄像机位置和姿态测量方法。该方法采用安装角度不固定的摄像机拍摄室内场景中一些具有低秩特征的纹理，如天花板上的方格纹理等，通过处理所拍摄图像进行摄像机空间位置和姿态6个自由度的测量。

本发明的技术方案是：

采用摄像机拍摄室内场景中的一些纹理图像，如天花板上的方格纹理等。用欧拉角表示低秩纹理成像模型中的投影矩阵，通过优化方法求解欧拉角表示的低秩纹理成像模型，解出表示摄像机姿态的欧拉角。再进一步利用纹理的几何特点和摄像机的成像关系，求解摄像机在空间中的位置。下面将对技术方案进行详细的介绍，首先介绍低秩纹理及摄像机对低秩纹理的拍摄，然后介绍如何利用纹理的低秩特性求解摄像机姿态，最后介绍求解摄像机空间位置的步骤。

步骤一：获取包含低秩纹理的图像

实际场景中的天花板、地面、墙壁上常可以抽象出如图1所示的纹理。这些纹理具有对称性和重复性。由于这些纹理的图像矩阵的秩较低，我们将这些纹理称为低秩纹理。这里给出低秩纹理的定义：通常二维(2D)的纹理被认为是在平面R²上定义的函数I⁰(x,y)；如果函数族I⁰(x,·)横跨有限低维子空间，则将纹理I⁰定义为低秩纹理。如式(1)所示，对于某较小的整数k：

$r=\mathrm{rank}\left(I^0\right)\stackrel{.}{=}\dim \left(\mathrm{span}\right\{x\left|I^0(x,\·)\right\})\≤k---\left(1\right)$

此时，I⁰被认为是秩-r纹理。通常，当k小于正在被评估的窗口中像素数目的一半，则被认为是低秩的。

本方法首先需要利用标定好的摄像机拍摄实际场景中的低秩纹理。标定好的摄像机是指通过一定的标定方法获得了焦距和主点值的摄像机。拍摄时，摄像机的安装角度可以不固定。摄像机可以一边调整自己的姿态拍摄，一边进行位置和姿态的测量。

步骤二：利用纹理的低秩特征求解摄像机姿态

低秩纹理在摄像机拍摄时存在投影形变和噪声，其所成的图像不再具备低秩特性。通过以下步骤可恢复出该低秩纹理并求得摄像机的拍摄姿态。

1)建立所需的坐标系

首先建立图像坐标系、摄像机坐标系、世界坐标系。如图2所示，图像坐标系以拍摄图像的中心为原点O，横轴OU和纵轴OV分别平行于图像边沿。摄像机坐标系以用于测量的摄像机光心O_c为原点，横轴O_cX_c和纵轴O_cY_c分别平行于图像坐标系OU轴和OV轴，竖轴O_cZ_c垂直于X_cO_cY_c平面。世界坐标系O_wX_wY_wZ_w根据不同的需要选取，一般将世界坐标系的X_wO_wY_w平面建立在低秩纹理特征所在的平面内，且世界坐标系各轴的方向选取使得表达低秩纹理的矩阵的秩最小。

2)建立摄像机成像模型，用欧拉角表示摄像机投影矩阵

选取低秩纹理上的一点P，设其在世界坐标系中坐标为(x_w,y_w,z_w,1)^T，在摄像机坐标系中的坐标表示为(x_c,y_c,z_c)^T，在图像平面上的像p的坐标为(u,v,1)^T。则在针孔成像模型下，P与p之间的关系可以描述为式(2)：

$s\·\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=N\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]=N\left[\begin{array}{cc}{R}_w^c& T_w^c\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]$ 其中 $N=\left[\begin{array}{ccc}f/\mathrm{dx}& 0& u_0\\ 0& f/\mathrm{dy}& v_0\\ 0& 0& 1\end{array}\right],R_w^c=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right],T_w^c=\left[\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right]---\left(2\right)$

式(2)中s为一个非零的比例因子，N为摄像机的内参数矩阵，f为摄像机的焦距；dx,dy分别表示U、V方向的像元尺寸，(u₀,v₀)为摄像机光心在图像坐标系中的坐标。分别表示摄像机相对于世界坐标系的旋转矩阵和平移矩阵。式(2)便是摄像机成像模型。

根据Euler定理，将式(2)中的旋转矩阵用三个欧拉角来表示。选择摄像机各轴的旋转顺序为：先绕Z轴旋转，再绕Y轴旋转，最后绕X轴旋转，顺时针旋转为正，对应的三个欧拉角分别为θ_x,θ_y,θ_z，则投影矩阵另写为：

$R_w^c=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]\left[\begin{array}{ccc}\cos {\mathrm{\θ}}_y\cos {\mathrm{\θ}}_z& -\cos {\mathrm{\θ}}_y\sin {\mathrm{\θ}}_z& \sin {\mathrm{\θ}}_y\\ \cos {\mathrm{\θ}}_x\sin {\mathrm{\θ}}_z+\sin {\mathrm{\θ}}_x\sin {\mathrm{\θ}}_y\cos {\mathrm{\θ}}_z& \cos {\mathrm{\θ}}_x\cos {\mathrm{\θ}}_z-\sin {\mathrm{\θ}}_x\sin {\mathrm{\θ}}_y\sin {\mathrm{\θ}}_z& -\sin {\mathrm{\θ}}_x\cos {\mathrm{\θ}}_y\\ \sin {\mathrm{\θ}}_x\sin {\mathrm{\θ}}_z-\cos {\mathrm{\θ}}_x\sin {\mathrm{\θ}}_y\cos {\mathrm{\θ}}_z& {\sin \mathrm{\θ}}_x\cos {\mathrm{\θ}}_z+\cos {\mathrm{\θ}}_x\sin {\mathrm{\θ}}_y\sin {\mathrm{\θ}}_z& \cos {\mathrm{\θ}}_x\cos {\mathrm{\θ}}_y\end{array}\right]---\left(3\right)$

式(3)便是欧拉角表示的投影矩阵。建立好各坐标系和摄像机成像模型后，摄像机的姿态和位置测量即是要求取摄像机坐标系各轴在世界坐标系中的欧拉角θ_x,θ_y,θ_z，以及摄像机光心O_c在世界坐标系中的坐标

3)建立基于欧拉角的低秩纹理成像模型

用I₀表示原始的低秩纹理图像，I表示实际拍摄的低秩纹理图像。如图3，将世界坐标系的X_wO_wY_w平面建立在低秩纹理特征所在的平面内，且世界坐标系各轴的方向选取使得表达低秩纹理的矩阵的秩最小。将I₀看做时摄像机在沿世界坐标系Z_w轴方向只有平移而没有旋转的位置进行拍摄得到的图像,简称为正对拍摄图像；将I看做摄像机在相对于世界坐标系存在某种旋转和平移的位置处进行拍摄得到的图像，简称为倾斜拍摄图像。根据摄像机模型和成像原理，当对低秩纹理上的一点P(x_w,y_w,z_w)进行正对拍摄时，其图像坐标(u₁,v₁)和世界坐标(x_w,y_w,z_w)之间的关系表示为式(4)；进行倾斜拍摄时，图像坐标(u₂,v₂)和世界坐标(x_w,y_w,z_w)之间的关系表示为式(5):

$s\·\left[\begin{array}{c}{u}_1\\ v_1\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_u& 0& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cc}{R}_1& T_1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]$ 其中 $R_1=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],T_1=\left[\begin{array}{c}0\\ 0\\ d\end{array}\right]---\left(4\right)$

$s\·\left[\begin{array}{c}{u}_2\\ v_2\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_u& 0& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cc}{R}_2& T_2\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]$ 其中 $R_2=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right],T_2=\left[\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right]---\left(\text{5}\right)$

其中R_i，T_i分别表示正对拍摄和倾斜拍摄时摄像机坐标系相对于世界坐标系之间的旋转矩阵和平移矩阵，i＝1,2。其中d表示正对拍摄时摄像机坐标系与世界坐标系之间的平移深度，r_mn为可用欧拉角表示的旋转矩阵中的对应参数，m,n＝1,2,3；t₁,t₂,t₃分别为平移矩阵中的三个平移分量。

结合(4)和(5)，求得(u₁,v₁)和(u₂,v₂)之间的关系如下：

$\begin{array}{c}{u}_2=\frac{f_d(r_{11}{u}_1+r_{12}{v}_1+\frac{f_d}{d}{t}_1)}{r_{31}{u}_1+r_{32}{v}_1+\frac{f_d}{d}{t}_3}=\frac{f_d(r_{11}{u}_1+r_{12}{v}_1+f_d{t^{\′}}_1)}{r_{31}{u}_1+r_{32}{v}_1+f_d{t^{\′}}_3}\\ v_2=\frac{f_d(r_{21}{u}_1+r_{22}{v}_1+\frac{f_d}{d}{t}_2)}{r_{31}{u}_1+r_{32}{v}_1+\frac{f_d}{d}{t}_3}=\frac{f_d(r_{21}{u}_1+r_{22}{v}_1+f_d{t^{\′}}_2)}{r_{31}{u}_1+r_{32}{v}_1+f_d{t^{\′}}_3}\end{array}---\left(\text{6}\right)$

其中t'＝[t'₁,t'₂,t'₃]＝[t₁/d,t₂/d,t₃/d]为归一化位移矢量。推导(6)时，不失一般性地假设了所拍摄的平面场景为世界坐标系的z＝0平面。

用符号ο表示式(6)中(u₁,v₁)和(u₂,v₂)之间运算，则两像素点之间的关系重写为：

(u₂,v₂)＝το(u₁,v₁)   (7)

其中τ为ο运算依赖的参数r_mn,t₁',t₂',t₃'的集合，m,n＝1,2,3。r_mn为世界坐标系到摄像机之间的旋转矩阵中的参数，已由欧拉角θ_x,θ_y,θ_z表示。由上述推导可知，正对拍摄图像I₀中的任意一点通过ο运算可以转化为倾斜拍摄图像中的一点I。将I₀和I上单个像素点的关系推广到图像所有像素点，并因为随机噪声E，则可推出I₀和I二者的关系如下：

I＝το(I₀+E)   (8)

其中τ为式(7)中ο运算依赖的参数集合。至此，低秩纹理的成像关系模型建立。

4)求解低秩纹理成像模型，解得摄像机的姿态

求解摄像机拍摄姿态即要求解τ。为求解τ，首先将式(8)将上述问题转化为如式(9)所示的优化问题：

其中，||E||₀表示E中非零元的个数。为求解该优化问题，将式(9)近似为：

其中，||I₀||_*表示I₀的核范数，||E||₁表示E的l₁-范数。式(10)中的约束条件是非线性约束，将其根据泰勒公式在点(u,v)处展开以实现约束条件线性化。由于u,v依赖于τ，故展开后有：

I(u(τ)+Δτ,v(τ)+Δτ)＝I(u(τ),v(τ))+▽I_τ(u(τ),v(τ))Δτ   (11)

其中▽I_τ是图像I关于τ中各个参数的雅克比矩阵。经过上述推导，将原优化问题重写为：

将式(12)所表达的优化问题利用迭代方法求解，具体步骤如下：

S1、接受所拍摄的低秩纹理图像作为输入图像，设定τ的初值为τ₀＝(0,0,0,0,0,1)，前三个参数表示欧拉角的初始值，后三个值表示归一化位移的初始值；选定收敛精度ε＞0；选定权值λ＞0；

S2、在输入图像上取含有低秩纹理的矩形窗，所得图像记为I；

S3、对图像I，重复以下迭代步骤，直到目标函数f＝||I⁰||_*+λ||E||₁全局收敛:

将图像I进行归一化，并将归一化的值重新赋给I；

对I求取关于欧拉角和归一化位移的雅克比矩阵；

使用增广拉格朗日乘子法求解式所表述的问题，得到τ的局部最优解τ'：

$\begin{array}{ccc}\underset{I_0,E,\mathrm{\Δ\τ}}{\min }{\left|\right|I_0\left|\right|}_{*}+\mathrm{\λ}{\left|\right|E\left|\right|}_1& s.t& I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))+{\▿}_{\mathrm{\τ}}(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\mathrm{\Δ\τ}=I_0+E---\left(13\right)\end{array}$

其中I(u_ij(τ),v_ij(τ))为用像素值(u,v)表示的图像I，u,v均为关于τ的函数；||I(u_ij(τ),v_ij(τ))||_F表示I(u_ij(τ),v_ij(τ))的F-范数；▽I_τ(u_ij(τ),v_ij(τ))为I关于τ的雅克比矩阵；

应用ο运算，将局部最优解τ'对应的投影变换作用在图像I的每一个点上，这样I便转化成了一幅新的图像；将所得的新图像重新赋值给I；

S4、输出全局最优解τ^*，其前三个值即为表示摄像机姿态的欧拉角的值；

通过迭代求解，可求得τ的最优解。从τ中可直接得到表示摄像机姿态的参数θ_x,θ_y,θ_z。需要指出的是，迭代所求的归一化位移值，不能表示摄像机所在位移的真实值，需要通过步骤三才能求解出摄像机的位移。

步骤三：利用成像原理和低秩纹理几何特性求解摄像机在空间中位置

求解摄像机在空间中位置分两步进行。首先，求取世界坐标系原点O_w在摄像机坐标系下的坐标，然后通过坐标变换求得摄像机光心O_c在世界坐标系中的坐标。摄像机光心O_c在世界坐标系中的坐标即是摄像机在空间中位置。具体步骤如下：

将光心O_c在世界坐标系中的坐标用表示，将世界坐标系原点O_w在摄像机系下的坐标用表示。选取低秩纹理所在平面上的两个特征点P1,P2，如图4。两点的选取以方便从图像上提取出这两点的像为原则。对正方形低秩纹理，可选取其顶点。两点之间的长度通过直尺测得，并将其转化为在世界坐标系中的坐标，用表示，其中i＝1,2。设P1,P2在摄像机坐标系中的坐标为从图像上提取出这两个特征点在图像坐标系上的坐标为(u_pi,v_pi)，结合式(2)求得和(u_pi,v_pi)满足如式(14)的关系：

$s\·\left[\begin{array}{c}{u}_{\mathrm{pi}}\\ v_{\mathrm{pi}}\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_d& 0& u_0\\ 0& f_d& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_x^c\\ r_{21}& r_{22}& r_{23}& t_y^c\\ r_{31}& r_{32}& r_{33}& t_z^c\end{array}\right]\left[\begin{array}{c}{x}_i^w\\ y_i^w\\ z_i^w\\ 1\end{array}\right]i=\mathrm{1,2}---\left(14\right)$

由(14)可解得

$t_x^c=(u_{p1}-u_0)(t_z^c+C_1)/f_d-A_1\text{,}$

$t_y^c=(v_{p1}-v_0)(t_z^c+C_1)/f_d-B_1\text{,}$

$t_z^c=\frac{f_d(A_1-A_2)-u_{p1}{C}_1+u_{p2}{C}_2+u_0(C_1-C_2)}{u_{p1}-u_{p2}}---\left(\text{15}\right)$

其中 $A_i=r_{11}{x}_i^w+r_{12}{y}_i^w+r_{13}{z}_i^w,B_i=r_{21}{x}_i^w+r_{22}{y}_i^w+r_{23}{z}_i^w,C_i=r_{31}{x}_i^w+r_{32}{y}_i^w+r_{33}{z}_i^w$ i＝1,2由于摄像机坐标系的坐标(x_c,y_c,z_c)和世界坐标系的坐标(x_w,y_w,z_w)可以通过旋转和平移相互转换，即存在如下关系：

$\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]=R_w^c\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\end{array}\right]+T_w^c=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\end{array}\right]+\left[\begin{array}{c}{t}_x^c\\ t_y^c\\ t_z^c\end{array}\right]---\left(16\right)$

将摄像机光心在摄像机坐标系中的坐标(0,0,0)代入式(16)，求得摄像机光心O_c在世界坐标系下的坐标

$\left[\begin{array}{c}{p}_x^w\\ p_y^w\\ p_z^w\end{array}\right]=-R_w^{c^{-1}}{T}_w^c=-\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right]\left[\begin{array}{c}{t}_x^c\\ t_y^c\\ t_z^c\end{array}\right]---\left(17\right)$

其中R_ij为中参数，i,j＝1,2,3。求得了摄像机光心O_c在世界坐标系下的坐标，即求得摄像机在空间中的位置。

本发明的效果和益处是，利用了室内场景中存在的低秩纹理特征，结合这些纹理特征的几何关系，求取摄像机在空间中的位置和姿态。该方法能够达到较高的精度，且在测量过程中，本方法可以转动摄像机，并实现了摄像机姿态和空间位置的六自由度测量，应用范围更加宽广、灵活。

附图说明

图1是实际场景中抽象出的低秩纹理图。其中图1(a)为条状低秩纹理图，图1(b)为棋盘格低秩纹理图，图1(c)为正方形低秩纹理。

图2是求解所需的坐标系示意图。其中1为世界坐标系O_wX_wY_wZ_w，用于描述现实世界中点的位置，根据自己需要建立。2为图像坐标系OUV，以CCD中心为原点O，横轴OU和纵轴OV分别平行于CCD的两个垂直边沿。3为摄像机坐标系，以摄像机光心O_c为原点，横轴O_cX_c和纵轴O_cY_c分别平行于图像坐标系OU轴和OV轴，竖轴O_cZ_c垂直于X_cO_cY_c平面。P为场景中低秩纹理上一点，其在图像平面上的像为p。

图3是对正拍摄与倾斜拍摄示意图。其中4为此时根据需要建立的世界坐标系，5为拍摄时的正方形纹理目标。6为正对拍摄时的摄像机坐标系，7为倾斜拍摄时的摄像机坐标系。

图4是利用场景中几何关系一直的低秩纹理求解摄像机空间三维坐标的示意图。图中1,2,3分别同图1中的1,2,3。P1,P2为低秩纹理的两个特征点，其在图像平面上的像为p1，p2。

图5是实施例中实验场景示意图。图中8为实验时所选世界坐标系，9为实验时所选摄像机坐标系，10为测量的出发点。实验时使用摄像机进行拍摄的时候，从10所指的点出发，沿着箭头所指方向，在*位置处拍摄。

具体实施方式

以下结合技术方案和附图详细叙述本发明的具体实施方式。可将利用低秩纹理图像求解摄像机姿态的解法归纳如下：

步骤一：选定可用于测量的低秩纹理。所选择的纹理应具有对称性或者重复性，一般应具有点或直线等特征，以方便从图像中识别。在室内场景中，天花板方格纹理、墙壁方格纹理、地板方格纹理均可选为用于测量的纹理。

选定好低秩纹理以后，在室内场景中放置摄像机，摄像机的位置和角度可以不固定，但需保证摄像机能拍摄到室内场景中的低秩纹理。调整摄像机焦距，使其聚焦。对调整好焦距的摄像机使用标定方法(如张正友标定法)进行标定，得到摄像机焦距大小和图像主点值。标定好摄像机后，摄像机的焦距不能再进行调节。使用标定好的摄像机拍摄室内场景中具有对称重复特征的低秩纹理图像。

步骤二：对拍摄图像进行处理，求解摄像机的拍摄姿态。在所拍摄的图像中选取包含低秩纹理的一块矩形区域作为输入，记为I。选取参数集τ的初始值τ₀＝(0,0,0,0,0,1)，权值λ＞0。为求取表示摄像机姿态的参数τ，对上述输入应用步骤如下：

1、接受输入：图像I,τ的初始值τ₀，权值λ

2、当目标函数f＝||I⁰||_*+λ||E||₁未全局收敛，重复如下循环：

S1:对图像I进行归一化处理，并将值赋给I，并求取I关于τ的雅克比矩阵▽I_τ(u_ij(τ),v_ij(τ))：

$I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\←\frac{I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))}{{\left|\right|I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\left|\right|}_F}$

${\▿}_{\mathrm{\τ}}I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\←\frac{\∂}{\∂\mathrm{\τ}}\left(\frac{I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))}{{\left|\right|I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\left|\right|}_F}\right)$

其中I(u_ij(τ),v_ij(τ))为用像素表示的图像I，||I(u_ij(τ),v_ij(τ))||_F表示I(u_ij(τ),v_ij(τ))的F-范数。

S2:求解如下问题，得到最优解E^*,Δτ^*：

$\begin{array}{ccc}(I_0^{*},E^{*},{\mathrm{\Δ\τ}}^{*})\←\underset{I_0,E,\mathrm{\Δ\τ}}{\min }{\left|\right|I_0\left|\right|}_{*}+\mathrm{\λ}{\left|\right|E\left|\right|}_1& s.t& I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))+{\▿}_{\mathrm{\τ}}(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\mathrm{\Δ\τ}=I_0+E\end{array}$

S3:更新参数τ:τ←τ+Δτ^*

3、输出:优化问题最终的最优解τ^*，并取τ^*的前3个值作为欧拉角的解

上述步骤中S2所需求解的问题需要用增广拉格朗日乘子法求解，其具体算法步骤归纳如下：

1、接受输入:I经过ο运算后的图像Iοτ，I关于τ的雅克比矩阵▽I，权值λ＞0，这里的τ为每一轮迭代的局部最优解。

2、代入初始值：Y₀＝0,E₀＝0,Δτ₀＝0,μ₀＞0,ρ＞1,k＝0.其中Y₀是拉格朗日乘子，E₀是E的初始值，▽τ₀是τ的初始更新步长，μ₀是增广拉格朗日函数中增广项系数μ的初始值，ρ为μ的更新因子，k为迭代次数；

3、依次按照下列步骤进行迭代计算直至||Iοτ-I₀-E+▽IΔτ||_F收敛：

$I_{k+1}^0\←U_k{S}_{{\mathrm{\μ}}_k^{-1}}\left[{\mathrm{\Σ}}_k\right]V_k^T$

μ_k+1＝ρμ_k

上述各式中，SVD(·)表示对括号内表达式进行奇异值分解；Σ_k为第k次迭代的缩放因子；Y_k、Y_k+1分别是拉格朗日乘子的第k次和第k+1次迭代值；E_k、E_k+1是随机噪声的第k次和第k+1次迭代值；▽τ_k、▽τ_k+1是τ的第k次和第k+1次迭代步长；μ_k、μ_k+1是增广拉格朗日函数中增广项系数的第k次和第k+1次迭代值；

4、输出：低秩图像的局部最优解I₀,随机噪声的局部最优解E,，步长的局部最优解Δτ。

补充说明，如果拍摄低秩纹理时，摄像机相对于世界坐标系的偏转角度不大，可以假设图像中的低秩纹理在投影前后面积不变，也即可在步骤二中的优化问题 $\begin{array}{ccc}\underset{I_0,E,\mathrm{\Δ\τ}}{\min }{\left|\right|I_0\left|\right|}_{*}+\mathrm{\λ}{\left|\right|E\left|\right|}_1& s.t& I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))+{\▿}_{\mathrm{\τ}}(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\mathrm{\Δ\τ}=I_0+E\end{array}$ 中添加条件▽S(τ)＝0，以使得求解结果更加准确。经过上述步骤，表示摄像机姿态的欧拉角参数得到求解。

步骤三：对拍摄图像进行特征提取，将纹理在场景中的几何特征和其在图像坐标系下的几何特征通过成像关系对应起来，求取摄像机光心在世界坐标系下的坐标，也即摄像机在空间中的位置。计算摄像机光心坐标的具体步骤如下：

1、在真实的低秩纹理上选取两个用于测量的特征点。选取特征点时以方便从图像上提取其像为原则，正方形纹理一般选取正方形纹理的顶点。用直尺测量出纹理上特征点之间的真实长度，并根据真实长度得出特征点在世界坐标系中的坐标。这里，特征点至少选取一对，也可以多选取几对。选取多对时，分别利用每一对特征点求取摄像机的位置，再将整合各对特征点的求取结果后的结果作为摄像机的真实位置。

2、利用图像处理方法，从图像上提取出所选纹理特征点在图像上对应的像素位置，得到特征点的像素坐标。所用的图像处理方法可根据纹理的不同等因素选择，一般用Harris角点提取算法提取角点，Hough变换提取直线，用Canny等算子提取边缘。对选择顶点为特征点的正方形纹理，可用Hough变换提取直线，求取直线的交点，便得到了特征点的像素坐标。

3、将步骤二中所求取的τ中的欧拉角参数带入式(3),计算出投影矩阵将所选特征点在世界坐标系中的坐标以及其在像素坐标系下的坐标带入式(15)和(17),计算摄像机在空间中的位置

实施例：

以下给出利用天花板上的正方形纹理进行实际测量的实施例。实验采用MCV1000SAM-HD相机，光学尺寸为1/2英寸，像元大小为5.2×5.2μm²,图像分辨率为1280*1024，接口为千兆网口，所使用镜头为4～10mm的变焦镜头。实验中天花板纹理为边长600mm的正方形，天花板距离地面3330mm，摄像机不固定角度地安装距离地面425mm的三脚架上。实验前进行摄像机内参数标定，得到摄像机焦距值和光心在图像平面上的坐标。实验时，实验场景示意图如图5，摄像机依照图5中标号1到24所示位置移动,对天花板进行拍摄。选取天花板方格的顶点为特征点，依次通过边缘提取、图像形态学开运算、形态学闭运算、Harris角点提取算法对拍摄的图像中的天花板纹理特征点提取，得到特征点的像素坐标，并结合其在世界坐标系中的坐标计算每幅图片对应的摄像机光心位置。测量所得摄像机在各拍摄点的欧拉角、各拍摄点在世界坐标系中的坐标测量值(x_i,y_i,z_i)与各拍摄点的标准参考值如表1所示。表中数据为多次测量的平均值。i表示测量点的序号，d_i,i-1表示序号i点和序号i-1点之间的测量距离，计算式如式(18),d⁰ _i,i-1表示序号i点和序号i-1点之间的真实距离,由直尺测量得到。绝对值误差为d_i,i-1与d⁰ _i,i-1之差的绝对值。

$d_{i,i-1}=\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2+{(z_i-z_{i-1})}^2}---\left(18\right)$

表1 利用天花板纹理测量摄像机位置和姿态的实验结果

Claims (2)

1.基于低秩纹理的摄像机位置和姿态测量方法，其特征在于，步骤如下：

步骤一：获取包含低秩纹理的图像

用标定好焦距和主点位置的摄像机对场景中的低秩纹理进行拍摄，获得包含低秩纹理的图像；

步骤二：利用拍摄的低秩纹理图像求解摄像机拍摄姿态

在CCD中心建立图像坐标系，在低秩纹理平面上建立世界坐标系，在摄像机光心处建立摄像机坐标系，并基于上述坐标系推导出低秩纹理上一点在三个坐标系中的坐标之间的关系，如式(2):

$s\·\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=N\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\end{array}\right]=N\left[\begin{array}{cc}{R}_w^c& T_w^c\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]$ 其中 $N=\left[\begin{array}{ccc}f/\mathrm{dx}& 0& u_0\\ 0& f/\mathrm{dy}& v_0\\ 0& 0& 1\end{array}\right],R_w^c=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right],T_w^c=\left[\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right]---\left(2\right)$

式(2)中(u,v)为低秩纹理上一点的像素坐标，(x_c,y_c,z_c)为该点在摄像机坐标系中的坐标，(x_w,y_w,z_w)为该点在世界坐标系中的坐标；s为一个非零的比例因子，N为摄像机的内参数矩阵，f为摄像机的焦距；dx,dy分别表示U、V方向的像元尺寸，(u₀,v₀)为摄像机光心在图像坐标系中的坐标,也即主点坐标；分别表示摄像机相对于世界坐标系的旋转矩阵和平移矩阵；

用欧拉角表示摄像机拍摄低秩纹理时的投影矩阵；并基于式(2)将低秩纹理成像模型修改为基于欧拉角的低秩纹理成像模型，将摄像机姿态的求解问题转化为一个优化问题:其中I₀表示原始的低秩纹理图像，I表示低秩纹理的实际拍摄图像，τ表示投影变换所依赖的参数集，包含6个参数,分别为摄像机坐标系各轴绕世界坐标系各轴旋转的三个欧拉角θ_x、θ_y、θ_z和表征摄像机光心相对于世界坐标系原点的归一化位移量t′₁,t′₂,t′₃；符号о表示将投影变换作用在I上的一种运算，E表示随机噪声，||E||₀表示随机噪声E中非零元的个数，λ为E的权值；

利用增广拉格朗日乘子法求解该优化问题，解出表示摄像机姿态的欧拉角参数；求解欧拉角的具体步骤如下：

S1、接受所拍摄的低秩纹理图像作为输入图像，设定τ的初值为τ₀＝(0,0,0,0,0,1)，其前三个参数表示欧拉角的初始值，后三个值表示归一化位移的初始值；选定收敛精度ε＞0；选定权值λ＞0；

S2、在输入图像上取含有低秩纹理的矩形窗，所得图像记为I；

S3、对图像I，重复以下迭代步骤，直到目标函数f＝||I⁰||_*+λ||E||₁全局收敛:

将图像I进行归一化，并将归一化的值重新赋给I；

对I求取关于欧拉角和归一化位移的雅克比矩阵；

用增广拉格朗日乘子法求解式(12)所表述的问题，得到τ的局部最优解τ':

$\begin{array}{ccc}\underset{I_0,E,\mathrm{\Δ\τ}}{\min }{\left|\right|I_0\left|\right|}_{*}+\mathrm{\λ}{\left|\right|E\left|\right|}_1& s.t& I(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))+\▿I_{\mathrm{\τ}}(u_{\mathrm{ij}}\left(\mathrm{\τ}\right),v_{\mathrm{ij}}\left(\mathrm{\τ}\right))\mathrm{\Δ\τ}=I_0+E\end{array}---\left(13\right)$

其中||I₀||_*表示I₀的核范数，||E||₁表示随机噪声E中非零元的个数的l1-范数，I(u_ij(τ),v_ij(τ))为用像素值(u,v)表示的图像I，且u,v均是关于τ的函数，为I关于τ的雅克比矩阵，Δτ为τ的增量，i,j分别表示I的行序号和列序号；

应用о运算，将局部最优解τ'对应的投影变换作用在图像I的每一个点上，这样I便转化成了一幅新的图像；将所得的新图像重新赋值给I；

S4、输出全局最优解τ^*，其前三个值即为表示摄像机姿态的欧拉角的值；

步骤三：利用成像原理和低秩纹理几何关系求解摄像机在空间中位置

在低秩纹理上选取两个特征点，测量其在空间中的几何关系，计算其在世界坐标系下的坐标，同时从摄像机拍摄的图像上提取特征点对应的像素点，得到特征点的像素坐标；利用针孔成像原理推导出摄像机位置的求解公式：

$\begin{array}{c}{t}_x^c=(u_{p1}-u_0)(t_z^c+C_1)/f_d-A_1,\\ t_y^c=(v_{p1}-v_0)(t_z^c+C_1)/f_d-B_1,\\ t_z^c=\frac{f_d(A_1-A_2)-u_{p1}{C}_1+u_{p2}{C}_2+u_0(C_1-C_2)}{u_{p1}-u_{p2}}\end{array}---\left(15\right)$

$\left[\begin{array}{c}{p}_x^w\\ p_y^w\\ p_z^w\end{array}\right]={-{R_w^c}^{-1}T}_w^c=-\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right]\left[\begin{array}{c}{t}_x^c\\ t_y^c\\ t_z^c\end{array}\right]---\left(17\right)$

其中表示世界坐标系原点O_w在摄像机坐标系下的坐标；表示摄像机光心O_c在世界坐标系中的坐标；(u_pi,v_pi)表示特征点的像素坐标；(u₀,v₀)为摄像机主点，f_d为摄像机像的焦距； $A_i=r_{11}{x}_i^w+r_{12}{y}_i^w+r_{13}{z}_i^w,B_i=r_{21}{x}_i^w+r_{22}{y}_i^w+r_{23}{z}_i^w$ $C_i=r_{31}{x}_i^w+r_{32}{y}_i^w+r_{33}{z}_i^w,$ 且i,j＝1,2；R_mn为投影矩阵的逆阵中参数，参数r_mn为投影矩阵中参数，m,n＝1,2,3,如式(3)；

将欧拉角的值、特征点在世界坐标系中的坐标和特征点的像在图像坐标系中的坐标带入公式(15)和(17)，计算得出摄像机的空间位置。

2.根据权利要求1所述的基于低秩纹理的摄像机位置和姿态测量方法，其特征在于，步骤二中о运算的步骤为将低秩纹理I₀上一点(u₁,v₁)通过式(6)，转化成图像I上的对应点(u₂,v₂)：

$\begin{array}{c}{u}_2=\frac{f_d(r_{11}{u}_1+r_{12}{v}_1+\frac{f_d}{d}{t}_1)}{r_{31}{u}_1+r_{32}{v}_1+\frac{f_d}{d}{t}_3}=\frac{f_d(r_{11}{u}_1+r_{12}{v}_1+f_d{t^{\′}}_1)}{r_{31}{u}_1+r_{32}{v}_1+f_d{t^{\′}}_3}\\ v_2=\frac{f_d(r_{21}{u}_1+r_{22}{v}_1+\frac{f_d}{d}{t}_2)}{r_{31}{u}_1+r_{32}{v}_1+\frac{f_d}{d}{t}_3}=\frac{f_d(r_{21}{u}_1+r_{22}{v}_1+f_d{t^{\′}}_2)}{r_{31}{u}_1+r_{32}{v}_1+f_d{f^{\′}}_3}\end{array}---\left(6\right)$

其中t'＝[t'₁,t'₂,t'₃]＝[t₁/d,t₂/d,t₃/d]为归一化位移矢量，t₁,t₂,t₃为摄像机真实的位移矢量；d为摄像机的平移深度，f_d为摄像机焦距；参数r_mn为投影矩阵中参数且m,n＝1,2,3,如式(3)：

$R_w^c=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_y{\cos \mathrm{\θ}}_z& -{\cos \mathrm{\θ}}_y{\sin \mathrm{\θ}}_z& \sin {\mathrm{\θ}}_y\\ {\cos \mathrm{\θ}}_x{\sin \mathrm{\θ}}_z+\sin {\mathrm{\θ}}_x{\sin \mathrm{\θ}}_y{\cos \mathrm{\θ}}_z& {\cos \mathrm{\θ}}_x{\cos \mathrm{\θ}}_z-{\sin \mathrm{\θ}}_x{\sin \mathrm{\θ}}_y{\sin \mathrm{\θ}}_z& -{\sin \mathrm{\θ}}_x\cos {\mathrm{\θ}}_y\\ {\sin \mathrm{\θ}}_x{\sin \mathrm{\θ}}_z-{\cos \mathrm{\θ}}_x{\sin \mathrm{\θ}}_y{\cos \mathrm{\θ}}_z& {\sin \mathrm{\θ}}_x{\cos \mathrm{\θ}}_z+{\cos \mathrm{\θ}}_x{\sin \mathrm{\θ}}_y{\sin \mathrm{\θ}}_z& {\cos \mathrm{\θ}}_x{\cos \mathrm{\θ}}_y\end{array}\right]---\left(3\right)$

其中θ_x，θ_y，θ_z为摄像机坐标系各轴绕世界坐标系各轴旋的三个欧拉角。