图像拼接是计算机视觉中的重要分支，它是将两幅以上的具有部分重叠的图像进行无缝拼接从而得到较高分辨率或宽视角的图像。图像拼接在虚拟现实、全景图制作、医学影像等领域起到重要作用。图像拼接必须完成关键的图像配准和图像融合两个过程。通过图像配准运算可以确定图像的旋转和平移等配准参数，通过图像融合运算可以使合成图像实现自然过渡，最终实现无缝拼接。但是，在实际应用中，由于待拼接数字图像存在曝光程度的差异、图像配准的误差、被摄物本身的形态变化以及光学成像过程产生的畸变等原因都会出现明显的缝合接缝现象，即，在视觉效果上能够明显看得出存在“接缝线”，能够清晰地分辨出处于“接缝线”两侧存在着明显色差或亮度的原始两幅图像，所谓“无缝拼接”就是拼接处具有平滑过渡的、富有良好视觉效果的、分辨不出存在任何接缝线的一幅完整图像，所以理想的图像拼接算法必须具备能够自动校正待拼接图像的畸变、消除待拼接图像的曝光差异和形变影响的技术性能，就当前的技术现状来看，现有的全景图拼接方法大多针对圆柱形和球形映射，能够有效应用于大型宽幅平面图像拼接的方法较为少见，功能也不完善。
提高数字图像自动拼接质量的技术关键在于如何克服图像形变和畸变所造成的影响，就当前的多幅图像拼接技术来看，通常对图片拍摄的焦距做出了限制，一般来说景物通过长焦距成像出现畸变的概率较小，通过广角即短焦距成像出现畸变的概率最大，因此至今为止，在该技术领域对广角图像进行自动拼接过程希望能够得到最终高质量的宽幅无缝拼接图像存在着较大的技术难度。
经对现有技术文献的检索发现，中国发明专利号为200510058832.X的专利，专利名称为：一种自动球形全景图的生成方法，该技术公开了一种自动球形全景图的生成方法。它通过图像对齐、投影变换、图像拉伸等操作，将一系列在某一固定视点拍摄的图像合成为一幅无缝的球形全景图。该技术不使用局部和全局优化技术来提高图像对齐精度，而是近似求解图像的旋转角度，并且允许对齐结果存在一定的误差积累。在生成球形全景图时，通过对图像进行分块拉伸的方法弥补旋转角度的误差，以达到无缝的球形拼接。但该技术的不足之处在于：1、不具备对图像的形变和畸变的自动纠正功能；2、图像经“竖置分块”和“水平分块”后映射至球面坐标，再通过竖直和水平方向的“拉升”，较难实现精确融合和无缝拼接；3、图像拼接的基本方法不适用于大型宽幅图像的自动拼接。

本发明的目的在于克服上述现有技术中的不足，提供一种广角图像自动拼接方法，使广角摄像机在移动过程中，能够对广角摄像机定时所拍摄的近距景物系列图像实现无缝拼接。
本发明是通过以下技术方案实现的，本发明包括如下步骤：
步骤一、设定广角摄像机移动速度V与摄像机图像采集周期T；
所述设定摄像机图像采集周期T，图像采集周期即相邻两幅图像采集的间隔时间，为了确保相邻图像具有足够的重叠部分，图像采集周期T的取值范围限定为：其中，u为每幅图像像素宽度，即每行像素数；Vu为对应摄像机移动速度V的像素移动速度。随着摄像机移动速度提高，采集周期T就要适当减小，亦即提高采样频率，才能确保相邻两幅图像具有足够的重叠部分；反之，采集周期T就要适当放大，亦即降低采样频率，才能确保相邻两幅图像不必要具备过大的重叠部分。
所述广角摄像机，其光学视场角达到120°以上。
为了可以降低景物图像形变发生率和拼接运算工作量，取V为常数，即摄像机采取匀速移动的方式，此时图像采集周期T为常数。
步骤二、根据设定的速度V和图像采集周期T拍摄图像，得到两幅相邻图像重叠宽度，具体为：Δu＝umax-TVu，单位为像素数pels，umax为采样图像最大宽度下的行像素数。
步骤三、建立图像非线性畸变矫正数学模型，并根据该数学模型中的参数构成一个畸变向量；
所述图像非线性畸变矫正数学模型，具体为：
$\left\{\begin{array}{c}\tilde{x}=x(1+k_1{r}^2+k_2{r}^4)+l_1(3x^2+y^2)+2l_2\mathrm{xy}+m_1(x^2+y^2)\\ \tilde{y}=y(1+k_1{r}^2+k_2{r}^4)+2l_1\mathrm{xy}+l_2(x^2+3y^2)+m_2(x^2+y^2)\end{array}\right.$
该方程考虑了光学镜头成像过程受多种畸变因素影响，忽略了高次项影响。其中，为图像平面坐标系中受畸变影响后的点坐标；r为图像平面坐标系中真实点坐标[x y]T到原点的距离，r2＝x2+y2；k1、k2为径向畸变参数，l1、l2为切向畸变参数，m1、m2为薄棱镜畸变参数；
所述畸变向量，由径向畸变参数、切向畸变参数和薄棱镜畸变参数构成，具体为：[k1 k2 l1 l2 m1 m2]T，该畸变向量建立在归一化虚平面上，不因摄像机的位置改变而改变，同时也与摄像机电荷耦合图像传感器(CCD)参数无关，当图像的分辨率改变时，内参数将发生变化，但是畸变向量参数不会变化；
所述归一化虚平面坐标系，是指：沿光轴在与摄像机光学成像平面形成对称的平面称为虚平面，在其上所建立的坐标系称为虚平面坐标系，将图像像素特征量进行归一化后在虚平面坐标系上的表达称为图像的归一化虚平面坐标系表达，使用归一化表达的虚平面坐标系所处平面，又称为归一化虚平面。
步骤四、实现图像数据坐标转换，包括：从世界坐标系到摄像机坐标系的转换、从摄像机坐标系到归一化虚平面坐标系的转换、从归一化虚平面坐标系到计算机图像像素坐标系的转换；
所述建立从世界坐标系到摄像机坐标系的转换，具体为：
$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_1\\ r_{21}& r_{22}& r_{23}& t_2\\ r_{31}& r_{32}& r_{33}& t_3\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]\Rightarrow \left[\begin{array}{c}{P}^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}P\\ 1\end{array}\right]$
其中，P＝[xw yw zw]T与P/＝[x′y′z′]T分别为世界坐标系与摄像机坐标系中点p和p′的坐标，称为三维旋转矩阵，t＝[t1 t2 t3]T为平移向量，由R与t构成的矩阵称为摄像机的外参数矩阵。
所述世界坐标系，即景物点的实际立体位置所处的地球空间三维坐标系，这是一个唯一的全局坐标系，坐标轴X、Y、Z的方向根据现场实际情况和讨论问题的出发点来决定，而景物点的实际位置用世界坐标系中的坐标P＝[x y z]T来表示，坐标单位以物理距离m为单位；
所述摄像机坐标系，即摄像机成像过程所使用的坐标系，其坐标轴X′、Y′、Z′的取向与摄像机密切相关，坐标系的原点在摄像机的光心处，其中，Z′轴取光轴的方向，X′和Y′建立在镜头平面上，方向与归一化虚平面坐标系x、y及其计算机图像像素坐标系u、v一致；
所述建立从摄像机坐标系到归一化虚平面坐标系的转换，具体为：
$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{z^{\prime }}\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]$
其中，[x y]T即归一化虚平面坐标系中未受畸变影响的真实图像点坐标，f为摄像机有效焦距。
所述建立从归一化虚平面坐标系到计算机图像像素坐标系的转换，具体为：
$\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}\frac{k_x}{d_x}x\\ \frac{1}{d_y}y\end{array}\right]+\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right]$
其中，kx为比例因子，用来适应在归一化虚平面坐标系横轴方向上取样带来的种种不确定因素，dx为计算机图像像素坐标系横轴方向上相邻两像素之间的有效距离，dy为计算机图像像素坐标系纵轴方向上相邻两象素之间的有效距离，[u v]T为计算机图像像素坐标系像素点坐标，[u0 v0]T为计算机图像像素坐标系中光学中心坐标。
u0、v0、kx、f、k1、k2、l1、l2、m1、m2均属于摄像机的内参数。
步骤五、根据图像数据坐标变换关系，通过迭代计算来求取畸变向量[k1 k2 l1 l2 m1 m2]T参数和摄像机的全部外参数值，具体为：
第一步：对计算机图像像素坐标系中光学中心坐标u0和v0，以及比例因子kx取初始预标定值，一般由产品说明书给出；
第二步：对其余内参数f、畸变向量[k1 k2 l1 l2 m1 m2]T参数、外参数采用迭代法求解，具体如下：
①将公式
$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_1\\ r_{21}& r_{22}& r_{23}& t_2\\ r_{31}& r_{32}& r_{33}& t_3\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]\Rightarrow \left[\begin{array}{c}{P}^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}P\\ 1\end{array}\right]$
代入公式
$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{z^{\prime }}\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]$
得到方程组
$\left\{\begin{array}{c}{x}_w\frac{{\mathrm{fr}}_{11}}{t_3}+y_w\frac{{\mathrm{fr}}_{12}}{t_3}+z_w\frac{{\mathrm{fr}}_{13}}{t_3}+\frac{{\mathrm{ft}}_1}{t_3}-x_{w}x\frac{r_{31}}{t_3}-y_{w}x\frac{r_{32}}{t_3}-z_{w}x\frac{r_{33}}{t_3}=x\\ x_w\frac{{\mathrm{fr}}_{21}}{t_3}+y_w\frac{{\mathrm{fr}}_{22}}{t_3}+z_w\frac{{\mathrm{fr}}_{23}}{t_3}+\frac{{\mathrm{ft}}_2}{t_3}-x_{w}y\frac{r_{31}}{t_3}-y_{w}y\frac{r_{32}}{t_3}-z_{w}y\frac{r_{33}}{t_3}=y\end{array}\right.$
令$R_1=\frac{{\mathrm{fr}}_{11}}{t_3},$$R_2=\frac{{\mathrm{fr}}_{12}}{t_3},$$R_3=\frac{{\mathrm{fr}}_{13}}{t_3},$$R_4=\frac{{\mathrm{fr}}_{21}}{t_3},$$R_5=\frac{{\mathrm{fr}}_{22}}{t_3},$$R_6=\frac{{\mathrm{fr}}_{23}}{t_3},$$R_7=\frac{r_{31}}{t_3},$$R_8=\frac{r_{32}}{t_3},$$R_9=\frac{r_{33}}{t_3},$$R_{10}=\frac{t_1}{t_3},$$R_{11}=\frac{t_2}{t_3},$并将其代入上式后表达为
$\left\{\begin{array}{c}x=(x_w{R}_1+y_w{R}_2+z_w{R}_3+R_{10})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\\ y=(x_w{R}_4+y_w{R}_5+z_w{R}_6+R_{11})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\end{array}\right.$
再采用最小二乘法求得中间参数R1～R11的最优解；
②将中间参数R1～R11代入
$\left\{\begin{array}{c}x=(x_w{R}_1+y_w{R}_2+z_w{R}_3+R_{10})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\\ y=(x_w{R}_4+y_w{R}_5+z_w{R}_6+R_{11})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\end{array}\right.$
求得x、y代入
$\left\{\begin{array}{c}\tilde{x}=x(1+k_1{r}^2+k_2{r}^4)+l_1(3x^2+y^2)+2l_2\mathrm{xy}+m_1(x^2+y^2)\\ \tilde{y}=y(1+k_1{r}^2+k_2{r}^4)+2l_1\mathrm{xy}+l_2(x^2+3y^2)+m_2(x^2+y^2)\end{array}\right.$
得到
$\left\{\begin{array}{c}x(k_1{r}^2+k_2{r}^4)+l_1(3x^2+y^2)+2l_2\mathrm{xy}+m_1(x^2+y^2)=\tilde{x}-x\\ y(k_1{r}^2+k_2{r}^4)+2l_1\mathrm{xy}+l_2(x^2+3y^2)+m_2(x^2+y^2)=\tilde{y}-y\end{array}\right.$
可见，上式中除了畸变向量[k1 k2 l1 l2 m1 m2]T外，其余均为已知，再使用最小二乘法就能够在迭代中获得畸变参数值；
③再转到第①步，重复上述迭代步骤，当迭代过程收敛到事先设定的标准时停止迭代，最终获得畸变参数的最优值；
④利用迭代第①～第③步的计算结果，由公式：t1＝t3R10/f、t2＝t3R11/f、r11＝t3R1/f、r12＝t3R2/f、r13＝t3R3/f、r21＝t3R4/f、r22＝t3R5/f、r23＝t3R6/f、r31＝t3R7、r32＝t3R8、r33＝t3R9同时计算获得全部外参数。
步骤六、将步骤五所确立的畸变向量[k1 k2 l1 l2 m1 m2]T参数值代入步骤三建立的图像非线性畸变矫正数学模型，通过该数学模型对受畸变影响后的图像点坐标进行求逆运算获得矫正非线性畸变后的真实点坐标[x y]T，进而获得矫正畸变后的整幅真实图像，因此实现对被采集景物图像的非线性畸变的自动消除，习惯上称这种矫正图像非线性畸变过程为“去畸变”。
步骤七、针对经过矫正畸变之后的图像，在重叠宽度为Δu＝umax-TVu的区域内快速寻找相邻两幅图像上的公共特征点，具体如下：
①利用基于图像灰度的角点检测方法在图像的重叠区域寻找4～5个角点作为特征点；
②将第一幅图像的重叠区域若干角点作为控制点，以每个控制点为中心建立4～5匹配模板，再分别以每个匹配模板对第二幅图像的重叠区域特征点进行搜寻，当被搜寻区域出现与匹配模板一致或者两者相似度大于等于95％，而且当前被搜寻区域内仅有一个角点时，即可确认两者的角点为公共特征点；重复如此搜寻与匹配过程，最终确定两幅图像上的4～5个公共特征点。
步骤八、将两幅图像进行透视变换恢复到同一个透视关系上，获得两幅图像形状一致的公共部分，并实现两幅图像在非重叠部分具有良好的形状平滑连续性，实现两幅图像的配准。
所述将两幅图像进行透视变换恢复到同一个透视关系上，具体如下：
①建立两幅图像的透视变换方程：
$w\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)=\left(\begin{array}{ccc}{t}_{11}& t_{12}& t_{13}\\ t_{21}& t_{22}& t_{23}\\ t_{31}& t_{32}& t_{33}\end{array}\right)\left(\begin{array}{c}\hat{x}\\ \hat{y}\\ 1\end{array}\right)$
其中，是透视变换矩阵，[x y]T分别为按照上式进行透视变换前后图像平面点坐标，即像素点坐标，w为比例因子；
矩阵T中的元素为透视变换前像素点坐标相对于透视变换后像素点坐标[x y]T的旋转参数，为前者相对后者的平移参数，[t31 t32]为透视失真参数，t33为整个图像的比例因子，规定为1；
②步骤七获得的公共特征点通过透视变换方程，将后一幅图像的特征点逐一与前一幅图像的对应特征点重合；
③将透视变换扩展至整个后一幅图像，对后一幅图像，从左上角的像素开始直至右下角的像素为止，除了公共特征点以外的所有像素点逐一转换成透视变换后的图像像素坐标图像。
步骤九、将已经配准好的两幅图像进行拼接，具体如下：
①分别以经过最左边和最右边公共特征点的两根竖直线作为两幅图像待拼接区域的左、右边线，即待拼接区域的左、右列像素，ul与ur代表对应左、右列的行坐标；
②对待拼接区域左、右边线宽度(ul-ur)进行N等分获得N+1根拼接线和N个列条块：令第一幅图像的待拼接图像拼接区域左、右边线宽度(ul-ur)图像为E，第二幅图像的待拼接图像拼接区域左、右边线宽度(ul-ur)图像为F，则被分块后的E依次得到列条块图像：E(1)、E(2)、…、E(N)，被分块后的F依次得到列条块图像：F(1)、F(2)、…、F(N)，此时，待拼接图像拼接区域左、右边线宽度(ul-ur)被N等分后，每列条块所包含的列像素条数
③在上述分块基础上，第一待拼接图像和第二待拼接图像的拼接区域左边线为L(0)，拼接区域每个分块的右边线依次为L(1)、L(2)、…、L(N-1)，拼接区域的右边线为L(N)，L(1)、L(2)、…、L(N-1)被分别称为第一、第二、…第N-1拼接线；
④第一待拼接图像和第二待拼接图像拼接线处的像素进行合成，具体如下：
$I\left(i\right)=\frac{1}{N}[(N-i)I_E\left(i\right)+i{I}_F\left(i\right)],i=\mathrm{1,2},...,N-1$
其中，I(i)表示E和F经过上述公式合成后第i条拼接线上的列像素值，像素值包括：黑白图像灰度值和彩色图像红色分量R、绿色分量G和蓝色分量B的亮度值，IE(i)、IF(i)分别代表待拼接区域宽度图像E和F在第i条拼接线上的原列像素值。
步骤十、采用基于小波双三次插值搜索方法对图像的拼接处进行融合，包括：拼接区域左边线与第一拼接线之间、相邻拼接线之间以及第N-1拼接线与拼接区域右边线之间的融合；
所述小波双三次插值搜索方法，具体为：原图像被小波分解成沿水平、垂直和对角三个方向的一系列高频信号带图像与一个低频图像，其中同一方向各高频信号带图像之间具有相似性，小波图像插值就是利用这种相似性，进行高频外推，然后利用小波反变换重构出比原始图像分辨力更高的插值图像，该方法在高频外推上采用双三次曲面拟合方法代替双线性方法做小波分解后高频外推的相似变换，由于双三次重建图像点的值和导数均连续，比较有效地抑制了双线性处理的高频失真，重建图像的细节丰富，畸变小，更适合人眼观察，因此小波双三次插值搜索算法是一种提高光学遥感图像空间分辨力的有效方法。具体实现过程如下：
①对第一待拼接图像拼接区域原图E进行分解
第一次对E分解：[EIL1，EIH1，EIV1，EID1]＝DWT(E)
第二次对EIL1分解：[EIL2，EIH2，EIV2，EID2]＝DWT(EIL1)
其中，EIH1和EIH2、EIV1和EIV2、EID1和EID2分别为水平、垂直、对角方向高频细节，EIL1和EIL2为低分辨率图像，DWT为离散小波变换；
②对EIL2、EIH2、EIV2、EID2进行小波反变换，重构高分辨率图像EIL1，
EIL1＝DWT-1(EIL2，EIH2，EIV2，EID2)
③再对EIL1、EIH1、EIV1、EID1进行小波反变换，重构高分辨率图像E，
E＝DWT-1(EIL1，EIH1，EIV1，EID1)
小波分解后水平、垂直、对角三个方向的细节图反映了原始图像在这三个方向上的边沿特征，由于边沿特征在同一个方向上的不同频带中的高频细节具有相似性；
④重复上述①～③的过程，对第二待拼接图像拼接区域原图F进行分解，先后获得
第一次对F分解：[FIL1，FIH1，FIV1，FID1]＝DWT(F)
第二次对FIL1分解：[FIL2，FIH2，FIV2，FID2]＝DWT(FIL1)
重构高分辨率图像FIL1：FIL1＝DWT-1(FIL2，FIH2，FIV2，FID2)
重构高分辨率图像F：F＝DWT-1(FIL1，FIH1，FIV1，FID1)
其中，FIH1和FIH2、FIV1和FIV2、FID1和FID2分别为水平、垂直、对角方向高频细节，FIL1和FIL2为低分辨率图像；
⑤根据E和F的插值结果对列像素逐列进行合成，合成基本公式如下：
$\left\{\begin{array}{c}I\left(i\right)=\frac{1}{N}[(N-i)I_E\left(i\right)+i{I}_F\left(i\right)],i=\mathrm{1,2},...,N-1\\ I_E\left(i\right)=I_E[{\mathrm{DWT}}^{-1}{(\mathrm{EIL}1,\mathrm{EIH}1,\mathrm{EIV}1,\mathrm{EID}1)]}_i\\ I_F\left(i\right)=I_F{\left[{\mathrm{DWT}}^{-1}(\mathrm{FIL}1,\mathrm{FIH}1,\mathrm{FIV}1,\mathrm{FID}1)\right]}_i\end{array}\right.$
其中，IE(i)、IF(i)分别代表E、F被等分后依次隶属第i条拼接线的经过小波双三次插值搜索运算重构的高分辨率列条块像素值；
对列条块图像E(1)、E(2)、…、E(N)和F(1)、F(2)、…、F(N)，按j＝1，2，...，SN进行再分解，将上述合成基本公式转换为
$\left\{\begin{array}{c}I\left(j\right)=\frac{1}{S_N}[(S_N-j)I_E\left(j\right)+j{I}_F\left(j\right)],j=\mathrm{1,2},...,S_N\\ I_E{\left(i\right)}^j=I_E{[{\mathrm{DWT}}^{-1}(\mathrm{EIL}1,\mathrm{EIH}1,\mathrm{EIV}1,\mathrm{EID}1{)}_i]}^{\left(j\right)}\\ I_F{\left(i\right)}^j=I_F{\left[{\mathrm{DWT}}^{-1}{(\mathrm{FIL}1,\mathrm{FIH}1,\mathrm{FIV}1,\mathrm{FID}1)}_i\right]}^{\left(j\right)}\end{array}\right.$
其中，IE(i)j、IF(i)j、分别代表E、F在第i列条块中第j列像素值，I(j)为合成后的第j列像素值，SN为每列条块所包含的列像素条数。
通过上述拼接处融合过程，能够克服图像畸变以及两幅图像因各自形变与曝光差异造成拼接部分与原图其余部分的明显区别，最终实现相邻两幅图像的无缝拼接。
当完成第一与第二两幅图像得无缝拼接后，继续上述步骤六至步骤十的过程完成第二幅与第三幅、第三幅与第四幅及其后续的所有相邻两幅图像的无缝拼接，最终实现一幅完整的、接近实际视觉效果的宽幅无缝拼接图像。
与现有技术相比，本发明具有如下有益效果：1、自动消除广角镜头引起的图像非线性畸变；2、对待拼接图像重叠部分自动进行图像精确配准；3、自动消除形变与曝光差异引起的两幅待拼接图像重叠部分在亮度与对比度上的差异；4、自动实现相邻图像的连续融合和无缝拼接。
附图说明
图1为本发明实施例采用的坐标体系图。

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
如图1所示，X、Y、Z为世界坐标系，对应物点的实际位置坐标为P＝[x y z]T，坐标单位以物理距离m为单位；X′、Y′、Z′为摄像机坐标系，即摄像机成像过程所使用的坐标系，坐标系的原点在摄像机的光心处，其中，Z′轴取光轴的方向，X′和Y′建立在镜头平面上，方向与归一化虚平面坐标系x、y及其计算机图像像素坐标系u、v一致，计算机图像像素坐标系点坐标为p＝[u v]T，单位为像素数，u′、v′为光学成像平面图像像素坐标系，对应像素点的坐标为p′＝[u′v′]T，单位为像素数。
实施前，设定世界坐标系、摄像机坐标系、光学成像平面图像像素坐标系、归一化虚平面坐标系和计算机图像像素坐标系，及其相互转换关系；安装视场角为120°的广角摄像机的支架沿着与世界坐标系的水平轴Y平行的导轨上移动，对摄像机前方的景物图像进行间隔采样拍摄，被拍摄画面总宽度在导轨允许的距离内随意确定，本实施例所能拍摄的画面总宽度长达700m。
本实施例包括如下具体步骤：
步骤一、采用光学视场角为120°的广角摄像机；设定摄像机沿着直线导轨的行走速度V＝10m/s，匀速移动；图像采集周期T＝3.06s。
步骤二、由V＝10m/s对应Vu＝226pels/s，以及T＝3.06s，得到两幅相邻图像重叠宽度：Δu＝umax-TVu＝76pels。
步骤三、建立图像非线性畸变矫正数学模型，具体为：
$\left\{\begin{array}{c}\tilde{x}=x(1+k_1{r}^2+k_2{r}^4)+l_1(3x^2+y^2)+2l_2\mathrm{xy}+m_1(x^2+y^2)\\ \tilde{y}=y(1+k_1{r}^2+k_2{r}^4)+2l_1\mathrm{xy}+l_2(x^2+3y^2)+m_2(x^2+y^2)\end{array}\right.$
待求畸变向量[k1 k2 l1 l2 m1 m2]T参数值。
步骤四、实现图像数据坐标转换，具体是：
①从世界坐标系到摄像机坐标系的转换为
$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_1\\ r_{21}& r_{22}& r_{23}& t_2\\ r_{31}& r_{32}& r_{33}& t_3\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]\Rightarrow \left[\begin{array}{c}{P}^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}P\\ 1\end{array}\right]$
②从摄像机坐标系到归一化虚平面坐标系的转换为
$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{z^{\prime }}\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]$
③从归一化虚平面坐标系到计算机图像像素坐标系的转换为
$\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}\frac{k_x}{d_x}x\\ \frac{1}{d_y}y\end{array}\right]+\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right]$
步骤五、根据图像数据坐标变换关系，通过迭代计算来求取畸变向量[k1 k2 l1 l2 m1 m2]T参数和全部外参数值，具体为：
第一步：给定u0、v0、kx的初始预标定值；
第二步：对其余内参数f、畸变向量[k1 k2 l1 l2 m1 m2]T参数、外参数
$\left[\begin{array}{cc}R& t\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_1\\ r_{21}& r_{22}& r_{23}& t_2\\ r_{31}& r_{32}& r_{33}& t_3\end{array}\right]$
进行迭代求解：
①：将公式
$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_1\\ r_{21}& r_{22}& r_{23}& t_2\\ r_{31}& r_{32}& r_{33}& t_3\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]\Rightarrow \left[\begin{array}{c}{P}^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}P\\ 1\end{array}\right]$
代入公式
$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{z^{\prime }}\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]$
得到方程组
$\left\{\begin{array}{c}{x}_w\frac{{\mathrm{fr}}_{11}}{t_3}+y_w\frac{{\mathrm{fr}}_{12}}{t_3}+z_w\frac{{\mathrm{fr}}_{13}}{t_3}+\frac{{\mathrm{ft}}_1}{t_3}-x_{w}x\frac{r_{31}}{t_3}-y_{w}x\frac{r_{32}}{t_3}-z_{w}x\frac{r_{33}}{t_3}=x\\ x_w\frac{{\mathrm{fr}}_{21}}{t_3}+y_w\frac{{\mathrm{fr}}_{22}}{t_3}+z_w\frac{{\mathrm{fr}}_{23}}{t_3}+\frac{{\mathrm{ft}}_2}{t_3}-x_{w}y\frac{r_{31}}{t_3}-y_{w}y\frac{r_{32}}{t_3}-z_{w}y\frac{r_{33}}{t_3}=y\end{array}\right.$
令$R_1=\frac{{\mathrm{fr}}_{11}}{t_3},$$R_2=\frac{{\mathrm{fr}}_{12}}{t_3},$$R_3=\frac{{\mathrm{fr}}_{13}}{t_3},$$R_4=\frac{{\mathrm{fr}}_{21}}{t_3},$$R_5=\frac{{\mathrm{fr}}_{22}}{t_3},$$R_6=\frac{{\mathrm{fr}}_{23}}{t_3},$$R_7=\frac{r_{31}}{t_3},$$R_8=\frac{r_{32}}{t_3},$$R_9=\frac{r_{33}}{t_3},$$R_{10}=\frac{t_1}{t_3},$$R_{11}=\frac{t_2}{t_3},$代入上式并将上式表达为
$\left\{\begin{array}{c}x=(x_w{R}_1+y_w{R}_2+z_w{R}_3+R_{10})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\\ y=(x_w{R}_4+y_w{R}_5+z_w{R}_6+R_{11})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\end{array}\right.$
采用最小二乘法求取中间参数R1～R11；
②：将中间参数R1～R11代入
$\left\{\begin{array}{c}x=(x_w{R}_1+y_w{R}_2+z_w{R}_3+R_{10})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\\ y=(x_w{R}_4+y_w{R}_5+z_w{R}_6+R_{11})/(x_w{R}_7+y_w{R}_8+z_w{R}_9+1)\end{array}\right.$
求得x、y后，代入
$\left\{\begin{array}{c}\tilde{x}=x(1+k_1{r}^2+k_2{r}^4)+l_1(3x^2+y^2)+2l_2\mathrm{xy}+m_1(x^2+y^2)\\ \tilde{y}=y(1+k_1{r}^2+k_2{r}^4)+2l_1\mathrm{xy}+l_2(x^2+3y^2)+m_2(x^2+y^2)\end{array}\right.$
得到
$\left\{\begin{array}{c}x(k_1{r}^2+k_2{r}^4)+l_1(3x^2+y^2)+2l_2\mathrm{xy}+m_1(x^2+y^2)=\tilde{x}-x\\ y(k_1{r}^2+k_2{r}^4)+2l_1\mathrm{xy}+l_2(x^2+3y^2)+m_2(x^2+y^2)=\tilde{y}-y\end{array}\right.$
再使用最小二乘法在迭代中获得畸变参数值；
③：再转到第①步，重复上述迭代步骤，当迭代过程收敛到事先设定的标准时停止迭代，最终获得畸变参数的最优值；
利用上述计算结果，最后再由公式：t1＝t3R10/f、t2＝t3R11/f、r11＝t3R1/f、r12＝t3R2/f、r13＝t3R3/f、r21＝t3R4/f、r22＝t3R5/f、r23＝t3R6/f、r31＝t3R7、r32＝t3R8、r33＝t3R9获得其余全部外参数值。
步骤六、根据所建立图像非线性畸变矫正数学模型和所确立的畸变向量[k1 k2 l1 l2 m1 m2]T参数值，对受畸变影响后的图像点坐标通过上述数学模型进行求逆运算获得矫正非线性畸变后的真实点坐标[x y]T，进而获得矫正畸变后的整幅真实图像，因此实现对被采集景物图像的非线性畸变的自动消除。
步骤七、针对经过矫正畸变之后的图像，在重叠宽度为Δu＝umax-TVu＝76pels的区域内快速寻找相邻两幅图像上的公共特征点，具体过程如下：
①利用基于图像灰度的角点检测方法在图像的重叠区域筛选4个角点作为特征点；
②将第一幅图像的重叠区域4个特征点作为控制点，以每个控制点为中心建立4匹配模板，再分别以每个匹配模板对第二幅图像的重叠区域特征点进行搜寻，当被搜寻区域出现与匹配模板一致或者两者相似度≥95％，而且当前被搜寻区域内仅有一个角点时，即可确认两者的角点为公共特征点；重复如此搜寻与匹配过程，最终确定两幅图像上的4个公共特征点。
步骤八、将两幅图像进行透视变换恢复到同一个透视关系上，获得两幅图像形状一致的公共部分，并实现两幅图像在非重叠部分具有良好的形状平滑连续性，实现两幅图像的配准。
具体过程如下：
①建立两幅图像的透视变换方程
$w\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)=\left(\begin{array}{ccc}{t}_{11}& t_{12}& t_{13}\\ t_{21}& t_{22}& t_{23}\\ t_{31}& t_{32}& t_{33}\end{array}\right)\left(\begin{array}{c}\hat{x}\\ \hat{y}\\ 1\end{array}\right)$
②根据所获得的公共特征点通过透视变换方程，将后一幅图像的特征点逐一与前一幅图像的对应特征点重合；
③将透视变换扩展至整个后一幅图像，对后一幅图像，从左上角的像素开始直至右下角的像素为止，除了公共特征点以外的所有像素点逐一转换成透视变换后的图像像素坐标图像。
步骤九、将已经配准好的两幅图像进行拼接，具体如下：
①分别以经过最左边和最右边公共特征点的两根竖直线作为两幅图像待拼接区域的左、右边线，即待拼接区域的左、右列像素，ul与ur代表对应左、右列的行坐标，计算得到(ul-ur)＝70；
②对待拼接区域左、右边线宽度(ul-ur)进行5等分获得6根拼接线和5个列条块，于是，待拼接图像拼接区域左、右边线宽度(ul-ur)＝70图像为E和F，被分别等分为竖直条块图像：E(1)、E(2)、…、E(5)和F(1)、F(2)、…、F(5)，待拼接图像拼接区域左、右边线宽度(ul-ur)＝70被5等分后，每列条块所包含的列像素条数SN＝14；
③根据上述分块结果，得到第一待拼接图像和第二待拼接图像的拼接区域的左边线为L(0)，拼接区域每个分块的右边线依次为L(1)、L(2)、…、L(4)，待拼接区域宽度的右边线为L(5)；
④按公式合成第一待拼接图像和第二待拼接图像拼接线像素值。
步骤十、采用基于小波双三次插值搜索方法对图像的拼接处进行融合，包括：拼接区域左边线L(0)与第一拼接线L(1)之间、相邻拼接线之间以及第4拼接线L(4)与拼接区域右边线L(5)之间的融合；
具体过程如下：
①对第一待拼接图像拼接区域原图E进行分解
第一次对E分解：[EIL1，EIH1，EIV1，EID1]＝DWT(E)
第二次对EIL1分解：[EIL2，EIH2，EIV2，EID2]＋DWT(EIL1)
其中，EIH1和EIH2、EIV1和EIV2、EID1和EID2分别为水平、垂直、对角方向高频细节，EIL1和EIL2为低分辨率图像，DWT为离散小波变换；
②对EIL2、EIH2、EIV2、EID2进行小波反变换，
由EIL1＝DWT-1(EIL2，EIH2，EIV2，EID2)重构高分辨率图像EIL1
③再对EIL1、EIH1、EIV1、EID1进行小波反变换，
由E＝DWT-1(EIL1，EIH1，EIV1，EID1)重构高分辨率图像E
④重复上述①～③的过程，对第二待拼接图像拼接区域原图F进行分解，先后获得
第一次对F分解：[FIL1，FIH1，FIV1，FID1]＝DWT(F)
第二次对FIL1分解：[FIL2，FIH2，FIV2，FID2]＝DWT(FIL1)
重构高分辨率图像FIL1：FIL1＝DWT-1(FIL2，FIH2，FIV2，FID2)
重构高分辨率图像F：F＝DWT-1(FIL1，FIH1，FIV1，FID1)
⑤根据E和F的插值结果通过合成基本公式
$\left\{\begin{array}{c}I\left(i\right)=\frac{1}{5}[(5-i)I_E\left(i\right)+i{I}_F\left(i\right)],i=\mathrm{1,2},...,5\\ I_E\left(i\right)=I_E[{\mathrm{DWT}}^{-1}(\mathrm{EIL}1,\mathrm{EIH}1,\mathrm{EIV}1,\mathrm{EID}1{)}_i]\\ I_F\left(i\right)=I_F\left[{\mathrm{DWT}}^{-1}{(\mathrm{FIL}1,\mathrm{FIH}1,\mathrm{FIV}1,\mathrm{FID}1)}_i\right]\end{array}\right.$
对列条块图像E(1)、E(2)、…、E(N)和F(1)、F(2)、…、F(N)，依次按进行再分解，将上述合成基本公式转换为
$\left\{\begin{array}{c}I\left(j\right)=\frac{1}{14}[(14-j)I_E\left(j\right)+j{I}_F\left(j\right)],j=\mathrm{1,2},...,14\\ I_E{\left(i\right)}^j=I_E{[{\mathrm{DWT}}^{-1}(\mathrm{EIL}1,\mathrm{EIH}1,\mathrm{EIV}1,\mathrm{EID}1{)}_i]}^{\left(j\right)}\\ I_F{\left(i\right)}^j=I_F{\left[{\mathrm{DWT}}^{-1}{(\mathrm{FIL}1,\mathrm{FIH}1,\mathrm{FIV}1,\mathrm{FID}1)}_i\right]}^{\left(j\right)}\end{array}\right.$
进而对逐列像素值进行合成，表示取不超过的整数值；
通过上述拼接处融合过程，克服了图像畸变以及两幅图像因各自形变与曝光差异造成拼接部分与原图其余部分的明显区别，最终实现相邻两幅图像的无缝拼接。
当完成第一与第二两幅图像得无缝拼接后，继续上述步骤六至步骤九的过程完成第二幅与第三幅、第三幅与第四幅及其后续的所有相邻两幅图像的无缝拼接，最终实现一幅完整的、接近实际视觉效果的宽幅无缝拼接图像。
本实施方法能够自动消除广角镜头引起的图像非线性畸变；实现待拼接图像重叠部分的自动精确配准；自动消除形变与曝光差异引起的两幅待拼接图像重叠部分在亮度与对比度上的差异；自动实现相邻图像的连续融合和无缝拼接；最终获得一幅700m宽的视觉效果良好的宽幅料场堆料图像。