CN101882306B - 一种表面凹凸物体照片的高精度拼接方法 - Google Patents

Abstract

本发明公开了一种表面凹凸物体照片的高精度拼接方法，根据需要拼接的第一图像和第二图像的重叠区域，利用三维扫描得到表面凹凸物体上若干采样点在世界坐标系中的坐标，用最小二乘法拟合得到最佳重投影平面，建立最佳重投影坐标系；把世界坐标系中的每个采样点的坐标转化为在最佳重投影坐标系中的坐标，结合相机内外参数计算出采样点的理想图像坐标和实际图像坐标之间的误差；通过插值得出重叠区域中所有像素的理想图像坐标和实际图像坐标间的误差；对第一图像和第二图像进行修正，得到消除误差的第一图像和第二图像并进行拼接。本发明通过相机成像所遵守的针孔模型模拟无穷大焦距镜头的成像，生成的照片就消除了表面凹凸透视投影带来的拼接误差。

Description

一种表面凹凸物体照片的高精度拼接方法

技术领域

本发明涉及计算机数字图像处理领域，尤其涉及一种表面凹凸物体照片的高精度拼接方法，通过三维扫描所获取的深度信息来消除这种普通相机照片由于透视投影所带来的凹凸表面的仿射对应误差。

背景技术

目前的图像配准算法都是基于这样的假设，即设(u，v)和(U，V)分别为同一物点在不同图像中的像点(对应点)，则一般认为满足如下隐函数关系

F(u，v，U，V)＝0                 (1)

注意：上式中函数F是一个向量函数，即可以写成如下的分量形式：

$\left\{\begin{array}{c}f(u,v,U,V)=0\\ g(u,v,U,V)=0\end{array}\right.---\left(2\right)$

如果能从隐函数解出显函数，那么就可以得到如下形式：

$\left\{\begin{array}{c}u=h(U,V)\\ v=i(U,V)\end{array}\right.---\left(3\right)$

或者

$\left\{\begin{array}{c}U=j(u,v)\\ V=k(u,v)\end{array}\right.---\left(4\right)$

至于h和i或者j和k具体是什么，不同的论文有不同的提法，但都逃不出仿射变换(包括平移变换、旋转变换、相似(比例)(缩放)变换等和透视(投影)变换的组合)。于是求解泛函h和i或者j和k就简化为了求解已知函数形式的系数问题。

从图1可以看到，由于相机成像遵守的是透视投影(因为不存在焦距无限大的镜头)，物体表面凹凸程度Δh会带来拼接误差。

首先计算相机光心到物体间的距离H。设期望的获取分辨率为R(PPI，Pixel Per Inch)，则折算成每毫米像素数为：

$\frac{R}{25.4}---\left(5\right)$

另一方面，设f_35mm为等效的全画幅135相机(感光元件大小为36毫米×24毫米)的焦距，若感光元件大小小于全画幅135相机的感光元件大小，设小的倍数为n(对角线距离的比值)，则相机的实际焦距f和等效135焦距存在如下关系，

f_35mm＝nf                    (6)

设水平方向的感光元件的分辨率为R_x，则感光元件每毫米像素数为：

$\frac{R_x}{\frac{36}{n}}=\frac{{\mathrm{nR}}_x}{36}---\left(7\right)$

根据透镜的物象放大比例关系，可得：

$\frac{\frac{{\mathrm{nR}}_x}{36}}{H}=\frac{\frac{R}{25.4}}{v}---\left(8\right)$

由于相机光心到物体间的距离远大于像距v，因此像距约等于焦距f，所以：

$\frac{\frac{{\mathrm{nR}}_x}{36}}{H}=\frac{\frac{R}{25.4}}{f}---\left(9\right)$

整理并将公式(6)代入得：

$H=\frac{\frac{{\mathrm{nR}}_x}{36}f}{\frac{R}{25.4}}=\frac{{127\mathrm{nfR}}_x}{180R}=\frac{{127f}_{35\mathrm{mm}}{R}_x}{180R}---\left(10\right)$

以水平方向的拼接为例，设p为相邻图像的重叠率(如完全重叠则p＝1)。

对于左图像，偏差为：

$E_l=\frac{\mathrm{\Δh}}{H+\mathrm{\Δh}}{R}_l=\frac{H}{H+\mathrm{\Δh}}\×\frac{\mathrm{\Δh}}{\frac{{127f}_{35\mathrm{mm}}{R}_x}{180R}}{R}_l=\frac{H}{H+\mathrm{\Δh}}\×\frac{180\mathrm{R\Δh}}{{127f}_{35\mathrm{mm}}}\×\frac{R_l}{R_x}---\left(11\right)$

类似的，对于右图像，偏差为：

$E_r=\frac{H}{H+\mathrm{\Δh}}\×\frac{180\mathrm{R\Δh}}{127f_{35\mathrm{mm}}}\×\frac{R_r}{R_x}---\left(12\right)$

因此，总偏差为：

$E_b=E_l+E_r=\frac{H}{H+\mathrm{\ΔH}}\×\frac{180\mathrm{R\Δh}}{127f_{35\mathrm{mm}}}\×\frac{R_l+R_r}{R_x}=\frac{H}{H+\mathrm{\Δh}}\frac{180\mathrm{R\Δh}}{127f_{35\mathrm{mm}}}\×(1-p)---\left(13\right)$

考虑到表面凹凸程度Δh远小于拍摄时光心和物体间的距离H，则：

$E_b=\frac{H}{H+\mathrm{\Δh}}\×\frac{180\mathrm{R\Δh}}{127f_{35\mathrm{mm}}}\×(1-p)\≈\frac{180\mathrm{R\Δh}}{127f_{35\mathrm{mm}}}\×(1-p)---\left(14\right)$

对于300PPI(即R＝300)的获取精度，设凹凸的最大值Δh为1毫米，135相机的等效焦距f_35mm为100毫米，图像的重叠率p为0.2，则误差为：

$\frac{180\×300\×1}{127\×100}\×(1-0.2)=3.4---\left(15\right)$

即表面每凹凸1毫米，误差就为3.4个像素，而实际凹凸远不止1毫米，如果按照5毫米计算，就是17个像素，完全不可接受！因此，必须有一种方法能够有效消除这种误差。

从公式(14)可以看出，误差和焦距成反比，和分辨率及表面凹凸程度成正比。因此，据此至少可以得到如下两个推论：

在高精度获取书画等表面起伏文物时此问题会突显。一般对于书画的获取要求至少为300PPI，而书画表面凹凸一般有数毫米，假设获取时的常用135相机的等效焦距f_35mm为100毫米，此时误差为17个像素，必须首先消除再拼接。

由于误差和焦距成反比，因此在拍摄中应该尽量采用长焦镜头。如果镜头焦距达到无穷大(此时实际上已经是平行投影了，这是由于实际拍摄时分辨率R是一定的，根据公式(10)，焦距f_35mm越大，要保持分辨率R不变，拍摄时的距离H就得越远，焦距无穷大时距离H也是无穷大，透视投影即成为平行投影。同时由于拍摄距离H是无穷大，此时有限大的表面凹凸程度Δh可以忽略不计，当作0，即表面没有凹凸的物体)，就不存在这种误差，但实际上此类镜头是不存在的。

发明内容

本发明提供了一种表面凹凸物体照片的高精度拼接技术，根据凹凸程度计算误差，然后在相邻图像的拼接前首先纠正这种误差然后再拼接，从而有效提高拼接精度。可以应用于表面凹凸物体的高精度获取，如书画等文物的高精度数字化等。

本发明所述的表面凹凸物体一般指书画、照片(印出来的非电子版的)等需要数字化存档的，大体上为平面结构，但表面略带凹凸(凹凸起伏的高度一般不超过1厘米)。

一种表面凹凸物体照片的高精度拼接方法，包括如下步骤：

(1)根据需要拼接的两幅图像(这两幅图像分别为第一图像和第二图像)中的重叠区域，利用三维扫描(例如采用结构光三维扫描仪)得到表面凹凸物体上若干采样点上的五元组(x_w，y_w，z_w，u_r，v_r)；

三维扫描时至少要求覆盖表面凹凸物体上与重叠区域相对应的部分，以便于误差的校正。

五元组中，(x_w，y_w，z_w)为采样点在三维空间(世界坐标系)中的坐标，(u_r，v_r)为采样点在第一图像中的坐标；

所述的采样点一般均匀分布在扫描区域中。

(2)根据所有采样点在三维空间中的坐标，用最小二乘法拟合得到最佳重投影平面；

以最佳重投影平面作为一个新的坐标系的xoy平面，建立坐标系，得到最佳重投影坐标系；

(3)计算世界坐标系和最佳重投影坐标系之间的旋转矩阵R′；

(4)把世界坐标系中的每个采样点的坐标左乘旋转矩阵R′，得到每个采样点在最佳重投影坐标系中的坐标(x_w’，y_w’，z_w’)，将步骤(1)中的五元组(x_w，y_w，z_w，u_r，v_r)更新，得到新五元组(x_w’，y_w’，z_w’，u_r，v_r)；

(5)根据新五元组(x_w’，y_w’，z_w’，u_r，v_r)，利用相机标定完成最佳重投影坐标系和图像坐标系之间所有参数的标定，得到拍摄第一图像时相机的内外参数；

(6)根据新五元组(x_w’，y_w’，z_w’，u_r，v_r)和拍摄第一图像时相机的内外参数，计算采样点的理想图像坐标(即假设需要拼接的两幅图像中的物体不存在凹凸部分)和实际图像坐标之间的误差；

(7)对步骤(6)的结果进行插值，得出所述的重叠区域中所有像素的理想图像坐标和实际图像坐标之间的误差，进行插值时可以采用现有技术中的算法。

由于采样点对应的图像坐标点(u_r，v_r)个数低于重叠区域的图像像素数，如果待纠正的图像点(u，v)没有对应的三维扫描数据点从而无法计算误差可用周围一些有三维数据点并且可以计算误差的点所计算出来的误差通过插值算法估算得出。

(8)利用重叠区域中所有像素的理想图像坐标和实际图像坐标之间的误差对第一图像进行修正，得到消除表面凹凸透视投影误差的第一图像，该图像是保证相邻图像间满足仿射变换的，参与后续拼接相对误差更小；

同理，根据步骤(1)～(8)，也可以得到消除表面凹凸透视投影误差的第二图像，区别仅在于(u_r，v_r)为采样点在第二图像中的坐标；而计算采样点的理想图像坐标和实际图像坐标之间的误差时，利用的是拍摄第二图像时相机的内外参数(根据步骤(5)标定)。

实际上，步骤(7)和步骤(8)可以通过纹理映射完成，即根据(x_w’，y_w’，z_w’，u_r，v_r)可以算出误差，进而得到消除误差后的(u_i，v_i)，将(u_r，v_r)和(u_i，v_i)分别作为输入图像和输出图像中的关键点坐标进行纹理映射完成误差纠正。

(9)将消除表面凹凸透视投影误差的第一图像和第二图像进行拼接。

步骤(2)中，基于最小二乘法的最佳重投影平面的拟合。该步骤的目的是在满足纠正表面凹凸透视投影误差的前提下尽可能接近原始图像(纠正量尽可能小，在不同的投影平面上误差大小是不同的)，因为图像的拉伸程度越大，对于细节的损失也就越大。因此需要在达到目的的前提下尽可能减少副作用。而根据公式(14)可以看到，误差(即纠正量)与凹凸程度成正比，因此选择恰当的坐标系，使得凹凸尽量小就可以保证误差所带来的纠正量尽可能小，而所谓凹凸就是点到某个平面的距离，凹凸尽量小可以通过最小二乘法实现，即最小化点到某个平面的距离的平方和。

步骤(3)中，最佳重投影平面与世界坐标系的对应关系的计算。在最小化凹凸时，对应着至少一个最佳重投影平面，即以这个平面衡量的凹凸是最小的，这一步主要是确定这个最佳重投影平面与世界坐标系的对应关系，因为后面计算误差纠正量时需要两者的关系。

本发明进行不平整表面的仿射变换误差的纠正，主要是根据相机成像所遵守的针孔模型，根据表面的凹凸信息，分别计算图像上每个像素按照实际镜头焦距拍摄所在的坐标和按照无穷大焦距拍摄(平行投影)所在的坐标，这两个坐标的差值就是误差，通过减掉误差就可以达到消除由于表面凹凸透视投影所带来的拼接误差问题。当然，实际上每个像素所对应物点的准确凹凸程度的获取很困难同时也没必要。通过获取一些特征点像素或者采样点像素上的准确凹凸程度，可以获得这些点上的拼接误差，而其余点的拼接误差可以通过插值方法得到，而通过纹理映射，就可以达到插值的目的。

上述方法可以推广到需要更多图像拼接的场合，既可以两张两张先校正后拼接，也可以统一先校正成平行投影成像的图像然后一起拼接。

本发明通过相机成像所遵守的针孔模型，模拟无穷大焦距镜头的成像，生成的照片就消除了表面凹凸带来的拼接误差。

附图说明

图1表面凹凸物体透视投影所带来的相邻图像的拼接误差原理图。

具体实施方式

下面结合实例详细描述本发明的一种具体实施方式。

(1)根据需要拼接的第一图像和第二图像中的重叠区域，利用结构光三维扫描仪得到表面凹凸物体上若干采样点上的五元组(x_w，y_w，z_w，u_r，v_r)；

这里的第一图像和第二图像通过拍摄表面凹凸物体得到，且有部分重叠区域，三维扫描时至少要求覆盖表面凹凸物体上与重叠区域相对应的部分，以便于误差的校正。

五元组中，(x_w，y_w，z_w)为采样点在三维空间(世界坐标系)中的坐标，(u_r，v_r)为采样点在第一图像中的坐标；。

(2)基于最小二乘法的最佳重投影平面的拟合

假设最佳重投影平面已经存在，则可用

AX+BY+CZ+D＝0                (16)

表示。此时，任意一点(x_i，y_i，z_i)到该平面的距离为：

$d_i=\left|\frac{{\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+D}{\sqrt{A^2+B^2+C^2}}\right|---\left(17\right)$

利用最小二乘法，可得所有点距离的平方和为：

$F=\underset{i}{\mathrm{\Σ}}{d}_i^2=\underset{i}{\mathrm{\Σ}}\frac{{({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+D)}^2}{A^2+B^2+C^2}---\left(18\right)$

上式中如果A、B、C和D都同时乘以一个非零的数并不会改变结果，因此，可以并且必须保证：

A²+B²+C²＝1                  (19)

且：

C≥0

                             (20)

此时，目标函数F可以简化为：

$F=\underset{i}{\mathrm{\Σ}}{({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+D)}^2---\left(21\right)$

如果令F取最小值，就可以得到最佳重投影平面，可用数学语言形式化地表示如下：

$\min \left(\underset{i}{\mathrm{\Σ}}{({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+D)}^2\right)$

$s.t.\left\{\begin{array}{c}{A}^2+B^2+C^2=1\\ C\≥0\end{array}\right.---\left(22\right)$

令

$\left\{\begin{array}{c}A=\cos \mathrm{\β}\sin \mathrm{\α}\\ B=\sin \mathrm{\β}\sin \mathrm{\α}\\ C=\cos \mathrm{\α}\\ 0\≤\∝\≤\frac{\mathrm{\π}}{2}\\ 0\≤\mathrm{\β}\≤2\mathrm{\π}\end{array}\right.---\left(23\right)$

则限制条件公式(19)和公式(20)等价为公式(23)，代入公式(21)可得：

$F=\underset{i}{\mathrm{\Σ}}{(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α}+D)}^2---\left(24\right)$

将F对α、β和D分别求偏导数得：

$\frac{\∂F}{\∂\mathrm{\α}}=2\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β\α}\sin +y_i\sin \mathrm{\β\α}\sin +z_i\cos \mathrm{\α}+D)(x_i\cos \mathrm{\β}\cos \mathrm{\α}+y_i\sin \mathrm{\β\α}\cos -z_i\sin \mathrm{\α})---\left(25\right)$

$\frac{\∂F}{\∂\mathrm{\β}}=2\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α}+D)(-x_i\sin \mathrm{\β}\sin \mathrm{\α}+y_i\cos \mathrm{\β}\sin \mathrm{\α})---\left(26\right)$

$\frac{\∂F}{\∂D}=2\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β\α}\sin +z_i\cos \mathrm{\α}+D)$

$=2[\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β\α}\sin +z_i\cos \mathrm{\α})+\mathrm{nD}]---\left(27\right)$

其中n表示待拟合的点的个数。取最值时，上述三个偏微分应该全部为零，即：

$\left\{\begin{array}{c}\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α}+D)(x_i\cos \mathrm{\β}\cos \mathrm{\α}+y_i\sin \mathrm{\β}\cos \mathrm{\α}-z_i\sin \mathrm{\α})=0\\ \underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α}+D)(-x_i\sin \mathrm{\β}\sin \mathrm{\α}+y_i\cos \mathrm{\β}\sin \mathrm{\α})=0\\ \underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α})+\mathrm{nD}=0\end{array}\right.---\left(28\right)$

其中由公式(28)中第三式可得：

$D=-\frac{\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α})}{n}---\left(29\right)$

因此，方程组(28)变为：

$\left\{\begin{array}{c}\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α}+D)(x_i\cos \mathrm{\β}\cos \mathrm{\α}+y_i\sin \mathrm{\β}\cos \mathrm{\α}-z_i\sin \mathrm{\α})=0\\ \underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α}+D)(-x_i\sin \mathrm{\β}\sin \mathrm{\α}+y_i\cos \mathrm{\β}\sin \mathrm{\α})=0\\ D=\frac{\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α})}{n}\end{array}\right.---\left(30\right)$

上述方程组的解析解很难求得，不过幸运的是，由于是在计算机上求解，因此可以用数值解代替。由方程组(30)中第二式可得：

$\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α}+D)(-x_i\sin \mathrm{\β}\sin \mathrm{\α}+y_i\cos \mathrm{\β}\sin \mathrm{\α})$

$=\underset{i}{\mathrm{\Σ}}(x_i\sin \mathrm{\α}\cos \mathrm{\β}+y_i\sin \mathrm{\α}\sin \mathrm{\β}+z_i\cos \mathrm{\α}+D)(-x_i\sin \mathrm{\α}\sin \mathrm{\β}+y_i\sin \mathrm{\α}\cos \mathrm{\β})---\left(31\right)$

将方程组(30)中第三式代入得：

$\underset{i}{\mathrm{\Σ}}(x_i\sin \mathrm{\α}\cos \mathrm{\β}+y_i\sin \mathrm{\α}\sin \mathrm{\β}+z_i\cos \mathrm{\α}-\frac{\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α})}{n})(-x_i\sin \mathrm{\α}\sin \mathrm{\β}+y_i\sin \mathrm{\α}\cos \mathrm{\β})$

$=\underset{i}{\mathrm{\Σ}}(x_i\sin \mathrm{\α}\cos \mathrm{\β}+y_i\sin \mathrm{\α}\sin \mathrm{\β}+z_i\cos \mathrm{\α})(-x_i\sin \mathrm{\α}\sin \mathrm{\β}+y_i\sin \mathrm{\α}\cos \mathrm{\β})$

$-\underset{i}{\mathrm{\Σ}}\left(\frac{\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α})}{n}\right)(-x_i\sin \mathrm{\α}\sin \mathrm{\β}+y_i\sin \mathrm{\α}\cos \mathrm{\β})$

$=\underset{i}{\mathrm{\Σ}}(x_i\sin \mathrm{\α}\cos \mathrm{\β}+y_i\sin \mathrm{\α}\sin \mathrm{\β}+z_i\cos \mathrm{\α})(-x_i\sin \mathrm{\α}\sin \mathrm{\β}+y_i\sin \mathrm{\α}\cos \mathrm{\β})$

$-\frac{\underset{i}{\mathrm{\Σ}}(x_i\cos \mathrm{\β}\sin \mathrm{\α}+y_i\sin \mathrm{\β}\sin \mathrm{\α}+z_i\cos \mathrm{\α})}{n}\underset{i}{\mathrm{\Σ}}(-x_i\sin \mathrm{\α}\sin \mathrm{\β}+y_i\sin \mathrm{\α}\cos \mathrm{\β})$

$=\underset{i}{\mathrm{\Σ}}(x_i{y}_i{\sin }^2\mathrm{\α}{\cos }^2\mathrm{\β}+y_i^2{\sin }^2\mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\β}+y_i{z}_i\cos \mathrm{\α}\sin \mathrm{\α}\cos \mathrm{\β}-x_i^2{\sin }^2\mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\β}-x_i{y}_i{\sin }^2\mathrm{\α}{\sin }^2\mathrm{\β}-x_i{z}_i\sin \mathrm{\α}\cos \mathrm{\α}\sin \mathrm{\β})$

$-\frac{(\cos \mathrm{\β}\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}+\sin \mathrm{\β}\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}+\underset{i}{\mathrm{\Σ}}{z}_i\cos \mathrm{\α})(\cos \mathrm{\β}\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}-\sin \mathrm{\β}\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α})}{n}$

$={\cos }^2\mathrm{\β}\underset{i}{\mathrm{\Σ}}{x}_i{y}_i{\sin }^2\mathrm{\α}+\sin \mathrm{\β}\cos \mathrm{\β}\underset{i}{\mathrm{\Σ}}(y_i^2-x_i^2){\sin }^2\mathrm{\α}-{\sin }^2\mathrm{\β}\underset{i}{\mathrm{\Σ}}{x}_i{y}_i{\sin }^2\mathrm{\α}+\cos \mathrm{\β}\underset{i}{\mathrm{\Σ}}{y}_i{z}_i\cos \mathrm{\α}\sin \mathrm{\α}-\sin \mathrm{\β}\underset{i}{\mathrm{\Σ}}{x}_i{z}_i\sin \mathrm{\α}\cos \mathrm{\α}$

$-\frac{{\cos }^2\mathrm{\β}\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}+\sin \mathrm{\β}\cos \mathrm{\β}[{\left(\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\right)}^2-{\left(\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\right)}^2]-{\sin }^2\mathrm{\β}\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}}{n}---\left(32\right)$

$-\frac{\cos \mathrm{\β}\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{z}_i\cos \mathrm{\α}-\sin \mathrm{\β}\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{z}_i\cos \mathrm{\α}}{n}$

$={\cos }^2\mathrm{\β}(\underset{i}{\mathrm{\Σ}}{x}_i{y}_i{\sin }^2\mathrm{\α}-\frac{\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}}{n})+\sin \mathrm{\β}\cos \mathrm{\β}(\underset{i}{\mathrm{\Σ}}(y_i^2-x_i^2){\sin }^2\mathrm{\α}-\frac{{\left(\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\right)}^2-{\left(\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\right)}^2}{n})$

${+\sin }^2\mathrm{\β}(\frac{\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}}{n}-\underset{i}{\mathrm{\Σ}}{x}_i{y}_i{\sin }^2\mathrm{\α})+\cos \mathrm{\β}(\underset{i}{\mathrm{\Σ}}{y}_i{z}_i\cos \mathrm{\α}\sin \mathrm{\α}-\frac{\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{z}_i\cos \mathrm{\α}}{n})$

$+\sin \mathrm{\β}(\frac{\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{z}_i\cos \mathrm{\α}}{n}-\underset{i}{\mathrm{\Σ}}{x}_i{z}_i\sin \mathrm{\α}\cos \mathrm{\α})$

$=P{\cos }^2\mathrm{\β}+Q\sin \mathrm{\β}\cos \mathrm{\β}+{R\sin }^2\mathrm{\β}+S\cos \mathrm{\β}+T\sin \mathrm{\β}$

其中：

$P=\underset{i}{\mathrm{\Σ}}{x}_i{y}_i{\sin }^2\mathrm{\α}-\frac{\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}}{n}---\left(33\right)$

$Q=\underset{i}{\mathrm{\Σ}}(y_i^2-x_i^2){\sin }^2\mathrm{\α}-\frac{{\left(\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\right)}^2-{\left(\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\right)}^2}{n}---\left(34\right)$

$R=\frac{\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}}{n}-\underset{i}{\mathrm{\Σ}}{x}_i{y}_i{\sin }^2\mathrm{\α}---\left(35\right)$

$S=\underset{i}{\mathrm{\Σ}}{y}_i{z}_i\cos \mathrm{\α}\sin \mathrm{\α}-\frac{\underset{i}{\mathrm{\Σ}}{y}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{\mathrm{\αz}}_i\cos \mathrm{\α}}{n}---\left(36\right)$

$T=\frac{\underset{i}{\mathrm{\Σ}}{x}_i\sin \mathrm{\α}\underset{i}{\mathrm{\Σ}}{z}_i\cos \mathrm{\α}}{n}-\underset{i}{\mathrm{\Σ}}{x}_i{z}_i\sin \mathrm{\α}\cos \mathrm{\α}---\left(37\right)$

而公式(32)所组成的方程可简化为：

Pcos²β+Qsinβcosβ+Rsin²β+Scosβ+Tsinβ＝0            (38)

可以按照如下方法求解，由方程(38)可得：

Pcos²β+R(1-cos²β)+Scosβ+(Qcosβ+T)sinβ＝0           (39)

由此可解得：

$\sin \mathrm{\β}=-\frac{{P\cos }^2\mathrm{\β}+R(1-{\cos }^2\mathrm{\β})+S\cos \mathrm{\β}}{Q\cos \mathrm{\β}+T}---\left(40\right)$

所以：

${(-\frac{P{\cos }^2\mathrm{\β}+R(1-{\cos }^2\mathrm{\β})+S\cos \mathrm{\β}}{Q\cos \mathrm{\β}+T})}^2+{\cos }^2\mathrm{\β}=1---\left(41\right)$

即：

[Pcos²β+R(1-cos²β)+Scosβ]²＝(1-cos²β)(Qcosβ+T)²    (42)

这是一个关于cosβ的四次方程(如果四次方系数正好为0，那么也可能不到四次)，而四次及以下方程具有求根公式，利用它可以求解得出cosβ进而进一步求解出sinβ。当然，由于进行过两边平方，有可能产生增根，因此需要把β代入原来的二次三角方程检验，以去除可能的增根。

总之，方程组(30)可以按照如下方法求解：

首先写一个单独的函数，根据α，利用方程组(30)后面两式求解β和D。该函数又可以具体分为如下几步：

1.根据α，利用公式(33)～公式(37)求得常数P、Q、R、S和T；

2.将P、Q、R、S和T代入公式(42)，求解关于cosβ的四次方程；

3.将cosβ和不确定符号的sinβ代入二次三角方程(38)检验，得到β；

4.将α和β代入方程组(30)的第三式得到D。

其次，由于

可以把这个区间分成若干个(比如等分成1000个)不相交的子区间，根据第一步的函数求出每个子区间分界点上的α、β和D，然后代入方程组(30)第一式左边，得到一个值g。如果子区间两端点的g值相反，则利用二分法可得这个子区间上唯一的一组解α、β和D。通过计算所有的可能子区间，最终可以找到所有的解α、β和D。

最后，对于所有的解α、β和D，根据A＝cosβsinα、B＝sinβsinα和C＝cosα分别求得A、B和C，并代入公式(24)可以计算出距离的平方和F。选取使得F最小的那组A、B、C和D就是期望的解，或者说这组解唯一地确定了一个最佳的重投影平面。

(3)最佳重投影坐标系与世界坐标系的对应关系的计算

在上一步中，已经利用世界坐标系中的点坐标数据(x_w，y_w，z_w)拟合出一个最佳重投影平面

AX+BY+CZ+D＝0，            (16)

这些点到这个平面的距离的平方和最小，但是这时，点到平面的距离并不是以zw的形式直接的表示出来，本步骤的目的，就是通过旋转变换，使得点到平面的距离为在新世界坐标系下的z的绝对值。为此，实际上只要通过旋转变换把平面方程(16)中的系数A和B消去即可，即通过如下的旋转变换，最终使得平面方程变成

C′Z+D′＝0，              (43)

的形式，进而得到如下形式：

$Z=-\frac{D^{\′}}{C^{\′}}---\left(44\right)$

也就是建立一个新的坐标系，使得重投影平面成为这个新坐标系的xoy平面。令

$\left[\begin{array}{c}X\\ y\\ Z\end{array}\right]=R\left[\begin{array}{c}{X}^{\′}\\ Y^{\′}\\ Z^{\′}\end{array}\right],---\left(45\right)$

其中，[X Y Z]^T和[X’Y’Z’]^T分别为世界坐标系中和最佳重投影坐标系中的坐标，而

$R=\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}1& 0& 0\\ 0& \frac{C}{\sqrt{B^2+C^2}}& \frac{B}{\sqrt{B^2+C^2}}\\ 0& -\frac{B}{\sqrt{B^2+C^2}}& \frac{C}{\sqrt{B^2+C^2}}\end{array}\right]\left[\begin{array}{ccc}\frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}& 0& \frac{A}{\sqrt{A^2+B^2+C^2}}\\ 0& 1& 0\\ \frac{-A}{\sqrt{A^2+B^2+C^2}}& 0& \frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}\end{array}\right]---\left(46\right)$

$=\left[\begin{array}{ccc}\frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}& 0& \frac{A}{\sqrt{A^2+B^2+C^2}}\\ -\frac{\mathrm{AB}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}& \frac{C}{\sqrt{B^2+C^2}}& \frac{B}{\sqrt{A^2+B^2+C^2}}\\ -\frac{\mathrm{AC}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}& -\frac{B}{\sqrt{B^2+C^2}}& \frac{C}{\sqrt{A^2+B^2+C^2}}\end{array}\right],$

将公式(46)代入公式(45)可得：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\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}^{\′}\\ Y^{\′}\\ Z^{\′}\end{array}\right]$

$=\left[\begin{array}{ccc}\frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}& 0& \frac{A}{\sqrt{A^2+B^2+C^2}}\\ -\frac{\mathrm{AB}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}& \frac{C}{\sqrt{B^2+C^2}}& \frac{B}{\sqrt{A^2+B^2+C^2}}\\ -\frac{\mathrm{AC}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}& -\frac{B}{\sqrt{B^2+C^2}}& \frac{C}{\sqrt{A^2+B^2+C^2}}\end{array}\right]\left[\begin{array}{c}{X}^{\′}\\ Y^{\′}\\ Z^{\′}\end{array}\right]---\left(47\right)$

$=\left[\begin{array}{c}\frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}{X}^{\′}+\frac{A}{\sqrt{A^2+B^2+C^2}}{Z}^{\′}\\ -\frac{\mathrm{AB}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}{X}^{\′}+\frac{C}{\sqrt{B^2+C^2}}{Y}^{\′}+\frac{B}{\sqrt{A^2+B^2+C^2}}{Z}^{\′}\\ -\frac{\mathrm{AC}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}{X}^{\′}-\frac{B}{\sqrt{B^2+C^2}}{Y}^{\′}+\frac{C}{\sqrt{A^2+B^2+C^2}}{Z}^{\′}\end{array}\right],$

代入公式(16)得：

$A(\frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}{X}^{\′}+\frac{A}{\sqrt{A^2+B^2+C^2}}{Z}^{\′})$

$+B(-\frac{\mathrm{AB}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}{X}^{\′}+\frac{C}{\sqrt{B^2+C^2}}{Y}^{\′}\frac{B}{\sqrt{A^2+B^2+C^2}}{Z}^{\′}),---\left(48\right)$

$+C(-\frac{\mathrm{AC}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}{X}^{\′}-\frac{B}{\sqrt{B^2+C^2}}{Y}^{\′}+\frac{C}{\sqrt{A^2+B^2+C^2}}{Z}^{\′})+D=0$

上式化简可得：

$\sqrt{A^2+B^2+C^2}{Z}^{\′}+D=0,---\left(49\right)$

满足公式(43)的形式，进一步化简可得：

$Z^{\′}=-\frac{D}{\sqrt{A^2+B^2+C^2}},---\left(50\right)$

由于公式(45)和(46)中的矩阵R是两个正交矩阵的乘积，因此R为正交矩阵，所以，

R^-1＝R^T。                (51)

由公式(47)可得：

$\left[\begin{array}{c}{X}^{\′}\\ Y^{\′}\\ Z^{\′}\end{array}\right]={\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]}^{-1}\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]={\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\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]$

$={\left[\begin{array}{ccc}\frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}& 0& \frac{A}{\sqrt{A^2+B^2+C^2}}\\ -\frac{\mathrm{AB}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}& \frac{C}{\sqrt{B^2+C^2}}& \frac{B}{\sqrt{A^2+B^2+C^2}}\\ -\frac{\mathrm{AC}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}& -\frac{B}{\sqrt{B^2+C^2}}& \frac{C}{\sqrt{A^2+B^2+C^2}}\end{array}\right]}^T\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]$

$=\left[\begin{array}{ccc}\frac{\sqrt{B^2+C^2}}{\sqrt{A^2+B^2+C^2}}& -\frac{\mathrm{AB}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}& -\frac{\mathrm{AC}}{\sqrt{A^2+B^2+C^2}\sqrt{B^2+C^2}}\\ 0& \frac{C}{\sqrt{B^2+C^2}}& -\frac{B}{\sqrt{B^2+C^2}}\\ \frac{A}{\sqrt{A^2+B^2+C^2}}& \frac{B}{\sqrt{A^2+B^2+C^2}}& \frac{C}{\sqrt{A^2+B^2+C^2}}\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]---\left(52\right)$

$=R^{\′}\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right],$

(4)综上，从世界坐标系到新坐标系的坐标变换方法为：通过把每个点的世界坐标系坐标左乘公式(52)中的旋转矩阵R′，得到每个采样点在最佳重投影坐标系中的坐标(x_w’，y_w’，z_w’)，将步骤(1)中的五元组(x_w，y_w，z_w，u_r，v_r)更新，得到新五元组(x_w’，y_w’，z_w’，u_r，v_r)；并在最佳重投影坐标系下得到简化后的可由公式(50)表示的平面方程以及点到重投影平面的距离；

(5)根据新五元组(x_w’，y_w’，z_w’，u_r，v_r)，利用相机标定完成最佳重投影坐标系和图像坐标系之间所有参数的标定，得到拍摄第一图像时相机的内外参数；

(6)假设照相机遵守针孔模型，那么图像坐标和世界坐标存在如下对应关系：

$Z_c\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{\mathrm{dx}}& 0& u_0\\ 0& \frac{1}{\mathrm{dy}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& 0\\ 0^T& 1\end{array}\right]\left[\begin{array}{cc}1& T\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]$

$=\left[\begin{array}{cccc}{a}_x& 0& u_0& 0\\ 0& a_y& v_0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& 0\\ 0^T& 1\end{array}\right]\left[\begin{array}{cc}1& T\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right],---\left(53\right)$

这儿(x_w，y_w，z_w)是一个物点的世界坐标，(u，v)是该点投影在图像平面的以像素为单位的坐标。a_x、a_y、u₀和v₀所在的矩阵为内参矩阵。(u₀，v₀)是基准点(通常在图像的中心，即图像宽度和高度的一半)，a_x和a_y是以像素为单位的焦距。R和T分别为旋转和平移矩阵也称作外参矩阵，将旋转矩阵R和平移向量T写成分量形式即得：

$Z_c\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{cccc}{a}_x& 0& u_0& 0\\ 0& a_y& v_0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cccc}{R}_{11}& R_{12}& R_{13}& 0\\ R_{21}& R_{22}& R_{23}& 0\\ R_{31}& R_{32}& R_{33}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& t_x\\ 0& 1& 0& t_y\\ 0& 0& 1& t_z\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]$

$=\left[\begin{array}{cccc}{a}_x& 0& u_0& 0\\ 0& a_y& v_0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cccc}{R}_{11}& R_{12}& R_{13}& 0\\ R_{21}& R_{22}& R_{23}& 0\\ R_{31}& R_{32}& R_{33}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_w+t_x\\ y_w+t_y\\ z_w+t_z\\ 1\end{array}\right]$

$=\left[\begin{array}{cccc}{a}_x& 0& u_0& 0\\ 0& a_y& v_0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{R}_{11}(x_w+t_x)+R_{12}(y_w+t_y)+R_{13}(z_w+t_z)\\ R_{21}(x_w+t_x)+R_{22}(y_w+t_y)+R_{23}(z_w+t_z)\\ R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)\\ 1\end{array}\right]---\left(54\right)$

$=\left[\begin{array}{c}{a}_x[R_{11}(x_w+t_x)+R_{12}(y_w+t_y)+R_{13}(z_w+t_z)]+u_0[R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)]\\ a_y[R_{21}(x_w+t_x)+R_{22}(y_w+t_y)+R_{23}(z_w+t_z)]+v_0[R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)]\\ R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)\end{array}\right],$

即：

$\left\{\begin{array}{c}u=\frac{a_x[R_{11}(x_w+t_x)+R_{12}(y_w+t_y)+R_{13}(z_w+t_z)]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)}+u_0\\ v=\frac{a_y[R_{21}(x_w+t_x)+R_{22}(y_w+t_y)+R_{23}(z_w+t_z)]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)}+v_0\end{array}\right.,---\left(55\right)$

如果被摄物体是几何上的平面物体，则z_w为0，那么理想的图像坐标(u_i，v_i)可以写成如下的形式：

$\left\{\begin{array}{c}{u}_i=\frac{a_x[R_{11}(x_w+t_x)+R_{12}(y_w+t_y)+R_{13}{t}_z]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}{t}_z}+u_0\\ v_i=\frac{a_y[R_{21}(x_w+t_x)+R_{22}(y_w+t_y)+R_{23}{t}_z]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}{t}_z}+v_0\end{array}\right.,---\left(56\right)$

此时相邻图像间满足仿射变换关系。但是，实际上，由于表面凹凸z_w的存在，图像上像素点的实际图像坐标(u_r，v_r)为：

$\left\{\begin{array}{c}{u}_r=\frac{a_x[R_{11}(x_w+t_x)+R_{12}(y_w+t_y)+R_{13}(z_w+t_z)]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)}+u_0\\ v_r=\frac{a_y[R_{21}(x_w+t_x)+R_{22}(y_w+t_y)+R_{23}(z_w+t_z)]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)}+v_0\end{array}\right.,---\left(57\right)$

两者之间的差别就会导致拼接误差，该误差为：

$\left\{\begin{array}{c}\mathrm{\Δu}=u_r-u_i=\frac{a_x[R_{11}(x_w+t_x)+R_{12}(y_w+t_y)+R_{13}(z_w+t_z)]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)}-\frac{a_x[R_{11}(x_w+t_x)+R_{12}(y_w+t_y)+R_{13}{t}_z]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}{t}_z}\\ \mathrm{\Δv}=v_r-v_i=\frac{a_y[R_{21}(x_w+t_x)+R_{22}(y_w+t_y)+R_{23}(z_w+t_z)]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}(z_w+t_z)}-\frac{a_y[R_{21}(x_w+t_x)+R_{22}(y_w+t_y)+R_{23}{t}_z]}{R_{31}(x_w+t_x)+R_{32}(y_w+t_y)+R_{33}{t}_z}\end{array}\right.---\left(58\right)$

理论上，只要修正公式(58)中的误差就可以消除由于表面凹凸z_w所导致的拼接误差。

根据公式(58)，利用新五元组(x_w’，y_w’，z_w’，u_r，v_r)和拍摄第一图像时相机的内外参数，得到采样点的理想图像坐标和实际图像坐标之间的误差Δu和Δv。

由于需要在满足仿射变换的前提下误差校正量尽量小，所以不是在原来的世界坐标系而是在最佳重投影平面为xoy平面的新的世界坐标系下完成校正工作。

(7)对步骤(6)的结果进行插值，得出所述的重叠区域中所有像素的理想图像坐标和实际图像坐标之间的误差，进行插值时可以采用现有技术中的算法；

(8)利用每个像素理想图像坐标和实际图像坐标之间的误差对重叠区域中第一图像所有像素进行误差修正，得到消除表面凹凸透视投影误差的第一图像；

同理，也可以得到消除表面凹凸透视投影误差的第二图像，区别仅在于(u_r，v_r)为采样点在第二图像中的坐标；而计算采样点的理想图像坐标和实际图像坐标之间的误差时，利用的是拍摄第二图像时相机的内外参数(根据步骤(5)标定)；

实际上，步骤(7)和步骤(8)可以通过纹理映射完成，即根据(x_w’，y_w’，z_w’，u_r，v_r)可以算出误差，进而得到消除误差后的(u_i，v_i)，将(u_r，v_r)和(u_i，v_i)分别作为输入图像和输出图像中的关键点坐标进行纹理映射完成误差纠正。

(9)将消除表面凹凸透视投影误差的第一图像和第二图像进行拼接。

Claims (2)

1.一种表面凹凸物体照片的高精度拼接方法，其特征在于，包括如下步骤：

(1)根据需要拼接的第一图像和第二图像中的重叠区域，利用结构光扫描仪进行三维扫描得到所述表面凹凸物体上若干采样点上的五元组(x_w，y_w，z_w，u_r，v_r)；所述的五元组中，(x_w，y_w，z_w)为采样点在世界坐标系中的坐标，(u_r，v_r)为采样点在第一图像中的坐标；

所述的第一图像和第二图像通过拍摄所述表面凹凸物体得到，所述的三维扫描至少覆盖所述表面凹凸物体上与所述的重叠区域相对应的部分；

(2)根据所有采样点在世界坐标系中的坐标，用最小二乘法拟合得到最佳重投影平面，以最佳重投影平面作为一个新的坐标系的xoy平面，建立最佳重投影坐标系；

(3)计算世界坐标系和最佳重投影坐标系之间的旋转矩阵R′；

(4)把世界坐标系中的每个采样点的坐标左乘旋转矩阵R′，得到每个采样点在最佳重投影坐标系中的坐标(x_w’，y_w’，z_w’)，将步骤(1)中的五元组(x_w，y_w，z_w，u_r，v_r)更新，得到与第一图像对应的新五元组(x_w’，y_w’，z_w’，u_r，v_r)；

(5)根据与第一图像对应的新五元组(x_w’，y_w’，z_w’，u_r，v_r)，利用相机标定完成最佳重投影坐标系和图像坐标系之间所有参数的标定，得到拍摄第一图像时相机的内外参数；

(6)根据步骤(5)的结果，针对第一图像，计算采样点的理想图像坐标和实际图像坐标之间的误差；

(7)针对第一图像，对步骤(6)的结果进行插值，得出所述的重叠区域中所有像素的理想图像坐标和实际图像坐标之间的误差；

(8)利用重叠区域中所有像素的理想图像坐标和实际图像坐标之间的误差，针对第一图像进行修正，得到消除表面凹凸透视投影误差的第一图像；

(9)对第二图像执行步骤(1)～(8)，得到消除表面凹凸透视投影误差的第二图像；将消除表面凹凸透视投影误差的第一图像和第二图像进行拼接。

2.如权利要求1所述的表面凹凸物体照片的高精度拼接方法，其特征在于，所述的步骤(2)中，最小二乘法拟合得到最佳重投影平面的过程如下：假设最佳重投影平面已经存在，计算任意一点到所述最佳重投影平面的距离后利用最小二乘法，得到所有点到所述最佳重投影平面的距离的平方和；对所述的所有点到所述最佳重投影平面的距离的平方和取最小值，得到最佳重投影平面。

Images