前言

又是编程作业。

工程文件请见:数字摄影测量

影像匹配

相关系数+最小二乘匹配

Mat相关知识

  • 行列与坐标的对应关系,经常容易混

opencv影像矩阵Matrows对应y方向,cols对应x方向。

构造Mat时,rows在前,cols在后:

Mat(int rows, int cols, int type);

使用Mat::at<_Tp>时,rows在前,cols在后:

intereMat.at<int>(row, col);//也就是(y,x)

使用坐标时,顺序是正常的:

Point2f point=Point2f(x,y);
  • 数据类型

    变量类型大小范围Mat中对应的数据类型
    uchar8bits0~255CV_8U
    char8bits-128~127CV_8S
    ushort16bits0~65535CV_16U
    short int16bits-32768~32767CV_16S
    int32bits-2147483648~2147483647CV_32S
    float32bits1.18e-38~3.40e38CV_32F
    double64bits2.23e-308~1.79e308CV_64F

    表格不全,但基本上可以看出类型的规律,数据类型以CV_+位数+有无符号(+波段数)组成。使用Mat::at<_Tp>时变量类型要与Mat的数据类型对应。

    在数据类型后加上C1C2C3,可以指明Mat存储数据的波段数,此时变量类型使用对应的向量类型,源码如下。

    typedef Vec<uchar, 2> Vec2b;
    typedef Vec<uchar, 3> Vec3b;
    typedef Vec<uchar, 4> Vec4b;

    typedef Vec<short, 2> Vec2s;
    typedef Vec<short, 3> Vec3s;
    typedef Vec<short, 4> Vec4s;

    typedef Vec<ushort, 2> Vec2w;
    typedef Vec<ushort, 3> Vec3w;
    typedef Vec<ushort, 4> Vec4w;

    typedef Vec<int, 2> Vec2i;
    typedef Vec<int, 3> Vec3i;
    typedef Vec<int, 4> Vec4i;
    typedef Vec<int, 6> Vec6i;
    typedef Vec<int, 8> Vec8i;

    typedef Vec<float, 2> Vec2f;
    typedef Vec<float, 3> Vec3f;
    typedef Vec<float, 4> Vec4f;
    typedef Vec<float, 6> Vec6f;

    typedef Vec<double, 2> Vec2d;
    typedef Vec<double, 3> Vec3d;
    typedef Vec<double, 4> Vec4d;
    typedef Vec<double, 6> Vec6d;

特征提取

采用的是Moravec特征提取算法,为了使经验阈值能不受窗口大小变化的影响,算法中计算兴趣值时进行了归一化(除以了窗口大小)。

相关系数匹配

从幻灯片里复制下来的左右影像对,其大小并不完全相同。因而先对影像进行了一次重采样,避免后续处理中的问题。

  • 用目标窗口直接对右影像进行整影像扫描

    对于第一组影像来说,由于左右影像基本上可以视为只有平移变换的影像,因而直接用形状相同的正方形窗口进行相关系数匹配,结果也很好。但是对于第二组影像来说,由于影像存在倾斜,匹配效果明显下降。此外,在第一步特征提取时也可以发现,Moravec算法提取的第二组的特征点兴趣值,也比第一组低,这可能也是匹配效果不佳的一个原因。

    此外,这种遍历整个影像的匹配算法,计算量大且浪费时间,匹配特征点的数量一旦变多,耗费时间就成倍增加,并不是一种可行的方案。

  • 优化算法,减少重复计算

    第一次扫描时,每次移动窗口,都把目标窗口(左影像上的窗口)重新计算了一遍,但是其实每次扫描中目标窗口都是不变的,所以进行了大量的无意义计算。同时,窗口移动时,右影像窗口的计算有一部分是重叠的,也可以进行优化。

    计算速度虽然变快了,但能计算的点的数量依然有限,耗费时间依然太长。

  • 优化算法,加入对上下视差和左右视差的预测

    注意到第一组影像基本上没有倾斜与旋转的影像,因而考虑可以由前若干次的扫描结果预测一个大概的上下视差和左右视差,然后在此基础上,预测右影像待匹配的窗口。

    匹配速度较上一次改进有很大的提升,并且匹配的效果也很好。

    但是第二对影像的匹配效果依然差得离谱。

    从理论上,对倾斜影像,如果依然采用正方形窗口(矩形窗口)而不根据倾斜的特点调整形状重采样的话,匹配效果是不可能有所改善的。因为相关系数本质上就是灰度分布相似的区域,而倾斜影像产生的像点位移会破坏掉原有的灰度分布(或者说是形状),导致相关系数失效。

  • 对于透视变换的重采样

    数字摄影测量课上所讲的核线重采样需要内方位元素,而题目所给数据并没有任何方位元素数据。虽然通过5对以上同名像点可以进行相对定向,但是缺少内方位元素不能直接运用。根据射影几何的知识,左右影像之间应该存在着射影变换关系。nn维射影变换可用n+2n+2对点确定,因此我们可以通过4对同名像点,确定左右影像像点之间的关系。将右影像按射影变换重采样到左影像上,然后相关系数匹配就可以使用相同大小的窗口进行匹配了。

    主要用到的函数:

    getPerspectiveTransform:由四对点计算透视变换矩阵。

    warpPerspective:根据透视变换矩阵,对图像进行透视变换。

    经过变换之后我们可以看到,这对影像基本上与第一对影像类似了。因此,我们可以使用对第一幅影像的处理方法进行处理。

    结果虽然有小部分错点,但效果比未变换时要好很多很多。

    参考:Opencv学习(15)—warpPerspective()_华山令狐冲、的博客-CSDN博客_warpperspective

最小二乘匹配

最小二乘匹配的前提是要有近似的匹配点作为初值。因此,需要先进行相关系数匹配之后,才能开始最小二乘匹配。

  • 几点需要注意的地方

    1. g2x2,g2y2\cfrac{\partial g_2}{\partial x_2},\cfrac{\partial g_2}{\partial y_2} 要怎么求?

      虽然是对x2,y2x_2,y_2求偏导,但实际上要理清楚,求导是在g2g2的影像上直接通过g2(x+1,y)g2(x1,y)2\cfrac{g_2(x+1,y)-g_2(x-1,y)}{2}计算,与x2,y2x_2,y_2是无关的,更不要用g2g_2几何纠正之后的图计算偏导数。

    2. 插值公式

      课件上的插值公式图中的xx轴和yy轴与opencv上的是相反的。

    3. 未知数更新的式子很奇怪

      和正常情况下的线性化然后更新的过程不一样,但是实际上确实收敛了。

  • 代码太乱,不放上来了

核线重采样

使用"水平"影像法制作核线影像,卫星摄影测量中讲得很仔细。

水平影像纠正关系

在解析摄影测量中,我们学习过倾斜影像与水平影像的关系,这种关系是基于共线条件方程建立的

[XXS YYS ZZS]=λR[x y f]\begin{equation} \begin{bmatrix}X-X_S \\\ Y-Y_S \\\ Z-Z_S\end{bmatrix} =\lambda R \begin{bmatrix}x \\\ y \\\ -f\end{bmatrix} \end{equation}

若要建立与水平影像的关系,就应该找到两者之间的联系——倾斜影像和水平影像对应的像点具有相同的地面坐标。

需要注意的是水平影像与倾斜影像的外方位线元素是相同的(纠正前后投影中心位置是不变的)

[XXS YYS ZZS]=λnRn[xn yn fn]\begin{equation} \begin{bmatrix}X-X_S \\\ Y-Y_S \\\ Z-Z_S\end{bmatrix} =\lambda_n R_n \begin{bmatrix}x_n \\\ y_n \\\ -f_n\end{bmatrix} \end{equation}

两式联立,求得

λnRn[xn yn fn]=λR[x y f]\begin{equation} \lambda_nR_n\begin{bmatrix}x_n \\\ y_n \\\ -f_n\end{bmatrix} = \lambda R\begin{bmatrix}x \\\ y \\\ -f\end{bmatrix} \end{equation}

对于水平影像,其RnR_n是单位矩阵。

xn=fna1x+a2ya3fc1x+c2yc3f yn=fnb1x+b2yb3fc1x+c2yc3f\begin{equation} \begin{aligned} x_n&=-f_n\cfrac{a_1x+a_2y-a_3f}{c_1x+c_2y-c_3f}\\\ y_n&=-f_n\cfrac{b_1x+b_2y-b_3f}{c_1x+c_2y-c_3f} \end{aligned} \end{equation}

平行于基线的影像纠正关系

课本中所谓的“水平”影像实际上指平行于基线的影像。关系式依然是式(3),但是RnR_n不能简单地用单位矩阵表示。而是应该根据基线计算出一个旋转矩阵,使得经这个旋转矩阵旋转后,影像恰好能平行于基线。

xn=fnm11x+m12ym13fm31x+m32ym33f yn=fnm21x+m22ym23fm31x+m32ym33f\begin{equation} \begin{aligned} x_n&=-f_n\cfrac{m_{11}x+m_{12}y-m_{13}f}{m_{31}x+m_{32}y-m_{33}f}\\\ y_n&=-f_n\cfrac{m_{21}x+m_{22}y-m_{23}f}{m_{31}x+m_{32}y-m_{33}f} \end{aligned} \end{equation}

卫星摄影测量课程中这种核线重采样方法讲得很详细,但是具体如何计算旋转矩阵,还是要靠自己理解,否则计算的角度可能正负会反。

建立了关系之后,利用间接法,对原始影像进行内插,灰度赋值,完成核线重采样工作。