- dataset: 三个测试数据集
- test_data：20张图片
- HeJiaDong：145张图片
- yongzhou：60张图片
- parfor_progress：在matlab terminal中创建进度条的函数
- vlfeat-0.9.14：VLFeat工具包，用于高速的SIFT检测和匹配
- adjustmentSparse.fm: 全局匹配单应矩阵求解，利用稀疏矩阵减少内存开销
- applyTransform.m: 对点进行单应矩阵变换
- findMaxConnected.m: 寻找匹配图最大连通分量
- GetImageMatches.m：获得两个图像的SIFT匹配
- imageWarping.m：图像扭曲
- lapMultiBlend.m：多个图像的拉普拉斯融合
- lapBlend.m：内存优化分批次的拉普拉斯融合
- linearBlend.m: Alpha融合
- main.m：主程序
- myimpyramid.m：改变的图像金字塔生成函数
- rebuildMask：生成融合时的mask
- resampleByOverlap.m：根据重叠面积大小选择图像
- 实现了基于SIFT特征算子的大范围图像拼接
- 对拼接速度和内存开销进行了优化
- 没有使用无人机相机的姿态和位置信息作为额外参数进行相片位置的确定
本实验首先利用图像的特征信息对无人机遥感图像进行特征点的提取与匹配。特征提取使用了尺度不变特征算子SIFT进行图像特征提取并进行特征点的匹配。然后通过RANSAC算法提取内点，删除错误的匹配，这也是图像配准过程中的一般性操作步骤和算法，接下来主要对这一过程进行详细讲解。
尺度不变特征转换即SIFT (Scale-invariant feature transform)是一种计算机视觉的算法。它用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量。其步骤简述如下：
- 尺度空间的极值检测：搜索所有尺度空间上的图像，通过高斯微分函数（DOG）来识别潜在的对尺度和旋转不变的兴趣点。
- 特征点定位：在每个候选的位置上，通过一个拟合精细模型来确定位置尺度，关键点的选取依据他们的稳定程度。
- 特征方向赋值： 基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向，后续的所有操作都是对于关键点的方向、尺度和位置进行变换，从而提供这些特征的不变性。
- 特征点描述： 在每个特征点周围的邻域内，在选定的尺度上测量图像的局部梯度，这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变换。
RANSAC是“RANdom SAmple Consensus（随机抽样一致）”的缩写。它可以从一组包含“局外点”的观测数据集中，通过迭代方式估计数学模型的参数。它是一种不确定的算法——它有一定的概率得出一个合理的结果；为了提高概率必须提高迭代次数。
RANSAC算法的思想简单而巧妙：首先随机地选择两个点，这两个点确定了一条直线，并且称在这条直线的一定范围内的点为这条直线的支撑。这样的随机选择重复数次，然后，具有最大支撑集的直线被确认为是样本点集的拟合。在拟合的误差距离范围内的点被认为是内点，它们构成一致集，反之则为外点。根据算法描述，可以很快判断，如果只有少量外点，那么随机选取的包含外点的初始点集确定的直线不会获得很大的支撑，值得注意的是，过大比例的外点将导致RANSAC算法失败。在直线拟合的例子中，由点集确定直线至少需要两个点；而对于透视变换，这样的最小集合需要有4个点。
图中蓝色点属于内点（正确点），而第红色点属于外点（偏移点）。此时用最小二乘法拟合这组数据，实际的最佳拟合直线是那条穿越了最多点数的蓝色实线。
在这一步中，我们需要对所有的图像进行特征的提取并对图像之间两两进行特征的匹配。这也是本实验中计算耗时最长的步骤。为了对这一步进行加速，我使用了VLFeat开源工具包，这也是公认的计算速度比较快的SIFT匹配工具包，同时我还利用了MATLAB的并行工具集对该过程进行了并行化的处理。
为了后面的对于图像单应变换矩阵的全局优化计算有效，我们需要先找出所有匹配中连通在一起的部分。因此我们需要构造图像的匹配图并寻找其中的最大连通分量
首先我们根据前一步中计算得到的图像特征匹配信息构造匹配图。尽管我们使用了RANSAC算法剔除错误匹配，但是在所有图像均两两匹配的情况，仍然可能存在毫不相关的图像由于特征点存在一定程度的线性关系而被视为已经匹配，因此我们需要设定一个最小匹配数量来对匹配图进行优化，在本实验中我们规定当特征点匹配对数大于全局所有图像配对中匹配对数的最大值的5%才算有效匹配，即：$N_{i,j}>thresh*\max{N} $。匹配图中每一个图像为一个顶点，如果两幅图像有效匹配，就将它们代表的顶点进行连接。
在得到匹配图后，我们从匹配图中寻找最大连通分量，并将它们对应的图像作为接下来进行全局优化计算的输入。
对于大范围的多幅拼接，如果逐一进行图像的配对和变换矩阵的计算，很有可能会因为累计误差而造成最后的大范围拼接结果十分不理想。为了规避这一问题，本实验基于单一的参考图像，构造目标函数，进行全局优化计算。其中参考图像的选取可以是任意的。
假设在$M$张无人机遥感影像中总共可以找到$N$组SIFT特征匹配。第i幅对于选定的参考图像变换矩阵为$H_i$,其对应的列向量为$X_i$，定义$\textbf{X} = (X_1^T,X_2^T, \cdots X_M^T)$。则我们需要优化的目标函数$E(\textbf{X})$为：
其中$\textbf{e}i = X_mp{i,m} - X_np_{i,n},,r_i = X_{ref}p_{i,ref}-p_{i,ref}\cdot(p_{i,m}, p_{i,n})$，$p_{i,m}$为第i组匹配中图像m中的特征点坐标，$ref$为参考图像的编号。这个目标函数实际上可以看做是方程$\textbf{A}\textbf{X}^T = \textbf{B}$，其中$\textbf{A}$为非参考图像匹配点构成的稀疏矩阵，$\textbf{B}$ 为参考图像匹配点坐标构成的稀疏矩阵。在实际求解的过程中由于$N$的值会很大，为了提高运算性能，节省内存，我们使用稀疏矩阵对$\textbf{A}$进行存储。
在具体的求解过程中，由于单应变换最终会构成也一个非线性的方程，稀疏矩阵的非线性方程在求解难度上比较大，因此我们使用正下视假定，也就是假设每一幅图像的视角都是固定垂直向下的，这样可以将单应变换变成仿射变换，进而解出线性方程。
基于我们前面得到的单应变换矩阵，我们可以将所有图像都旋转到和参考图像相对应的平面上，并将它们进行融合。融合主要分为两步：融合预处理和拉普拉斯融合
因为拉普拉斯融合需要构建图像金字塔，对于内存的消耗会比较大，同时无人机影像之间的重叠度可能比较大，过多的图像同时做拼接可能会导致鬼影较多的问题。为了提高拼接的准确度，我们在这一步还需要对图像做进一步的筛选。我们对于变换后的图像之间的重叠面积进行计算，删除与其他图像重叠部分面积较大的图像，得到最终用于拼接的图像数据集。
除此之外，我们还对拉普拉斯融合的掩膜进行了重新计算，让不同图像的mask之间尽量没有重叠，便于我们得到最终的融合结果。首先我们可以获得每一幅图像在经过单应变换后的图像范围，以及最终融合图像的总大小；然后我们构建一个与融合图像大小相同的矩阵，对于融合图像内的每一个像元做如下操作
- 如果该像元没有位于任何一个变换后图像的范围，矩阵的对应位置记为0；
- 如果该像元位于某些变换后图像的范围内，则矩阵对应位置的值为该像元距离图像边界最小值的最大值所对应的图像
拉普拉斯融合分为以下几步：
- 建立图像的拉普拉斯金字塔和高斯金字塔
- 求掩膜图像的高斯金字塔；
- 对图像的拉普拉斯金字塔进行拼接：在每一层上将图像透过掩膜的部分拼接得到结果的拉普拉斯金字塔
- 计算结果的图像金字塔顶层：将每张图像的高斯金字塔顶层透过对应掩膜金字塔顶层的部分拼接
- 重建图像：从结果的图像金字塔顶层开始，将结果的拉普拉斯金字塔的每一层相加然后上采样，最后得到融合结果。
但是拉普拉斯融合由于需要构建多级金字塔，因此存在很大的内存开销，在图片数目较多时，我机器的内存无法容纳全部的图像同时进行拉普拉斯融合，构建金字塔。为了解决这一问题，我采用了一种分批次的拉普拉斯融合，使得融合可以勉强进行。
内存优化的分批次融合分为以下几步：
- 根据用户输入的分批次数对图像进行分批次处理，如果用户输入的为n，那么会产生2n个拉普拉斯融合的图像批次，每个批次的大小为img_num/n，第i个批次的图像范围是i*img_num/2n~(i+1)img_num/2n,这样可以最大程度保证图像之间的相互覆盖，减少不同批次图像之间的由于上下采样过程中的误差而产生的接缝。
- 对所有批次进行拉普拉斯融合
- 对融合结果进行alpha融合
[1] Xu Y , Ou J , Hu H , et al. Mosaicking of Unmanned Aerial Vehicle Imagery in the Absence of Camera Poses[J]. Remote Sensing, 2016, 8(3):204.