关于遥感图像的控制点片匹配算法的一点想法（一）

    在遥感图像的几何纠正过程中，可能会用到控制点库的点片自动匹配。关于控制点片匹配的算法有很多种，从最基本的分类可以分为基于像元的点片匹配和基于特征的点片匹配，由于基于特征的算法难度较大，一般使用的都是基于像元的。

首先简单说明一下，点片匹配在数字图像处理中交模板匹配（Match Template），模板匹配是数字图像处理的重要组成部分之一。把不同传感器或同一传感器在不同时间、不同成像条件下对同一景物获取的两幅或多幅图像在空间上对准,或根据已知模式到另一幅图中寻找相应模式的处理方法就叫做模板匹配。简单而言，模板就是一幅已知的小图像。模板匹配就是在一幅大图像中搜寻目标，已知该图中有要找的目标，且该目标同模板有相同的尺寸、方向和图像，通过一定的算法可以在图中找到目标，确定其坐标位置。【参考1】

    以下图为例简单说明一下，左边淡蓝色区域表示待搜索的图像，右上角的三角形表示目标，意思就是要在左边的图中找出右边的正三角形在左边图中位置。

    一般的算法思想：将搜索模板 T( m×n 个像素)叠放在被搜索图S( W×H 个像素)上平移，模板覆盖被搜索图的那块区域叫子图Sij。i，j为子图左上角在被搜索图S 上的坐标。搜索范围是：
       1≤ i ≤ W–m
       1≤ j ≤ H–n
    通过比较T 和Sij 的相似性，完成模板匹配过程。 衡量模板T和子图Si,j的匹配程度,可用下列两种测度:

或者

    将上面第一个式子展开得到下面的：

    从展开的公式中可以看出，中间第二项是一个常数，也就是只跟模板有关系，而第一个和第三个是和原图有关的，随着模板在原图上的移动，这两个值也在发生变化。当D(i,j)的值最小时，说明找到了目标，通过这个关系，可以将展开的公式变形得到其他的公式，但是原理都是根据这个来的。在此不再一一列举。

    这个算法的效率是非常的地下的，每移动一个像元就要进行一次计算，计算量是相当的大，曾经试验过一个1000*1000的图，模板大小100*100，大概需要十分钟左右（机器性能也一般，但是这种真的是太慢了）。于是就有人提出下面的优化算法。主要有：

1 粗精匹配结合
观察实际模板匹配运算结果可以发现,匹配点附近的匹配误差迅速下降,明显区别于其它位置。针对这一特点,采用粗精匹配结合的算法迅速锁定匹配点大致区域,可大大降低整体匹配次数。具体实现方法:先跳动着隔几个点进行一次粗匹配,大致框定匹配区域,然后在附近区域逐一检索获得最佳匹配点。运算量可减少到三分之一以下，且目标提取效果相当好。

2 限制最大匹配误差
因为只需找到最小匹配误差的位置,不必完整计算每一位置的绝对匹配误差,而以已经计算的最小匹配误差作为最大允许误差。若计算误差大于该最大允许误差,就肯定不是最佳匹配点,可以提前结束计算,进入下一匹配位置的计算;如果匹配完成后仍小于最大允许误差,就用当前误差替换最大允许误差,并把该点作为潜在的匹配位置记录下来。匹配点和非匹配点的误差常常相差2~3个数量级。经过这种处理后,匹配点后剩余的计算量可以大大降低。

3 乱序匹配
目标出现在匹配区域中的位置不确定。不固定顺序算法可以更快地检索到匹配区域,迅速降低最大匹配误差,减少剩余非匹配点的计算量,降低整体运算量。针对光电探测设备的实际工作情况,在跟踪状态下,目标位移角速度和角加速度有限,导致目标常处于匹配区域中心附近。选择由中心向周围辐射匹配的方式效果最理想。

   程序就不写了，很简单的，可以参考【参考1】的链接。此外在OpenCV库中提供了模板匹配的算法（MatchTemplate），详细可以参考OpenCV的帮助文档，此处就其函数进行一个简单的说明：

[cpp] view plaincopyprint?
CV_EXPORTS_W void cv::matchTemplate (const Mat & image,  
                                     const Mat & templ,  
                                     /*CV_OUT*/ Mat & result,  
                                     int  method);  
 

    四个参数的意思很好懂，第一个是搜索图像，第二个是模板图像，第三个是搜索的结果图像，第四个参数是匹配方式，可以取值有下面六种，分别对应不同的算法。   

    CV_TM_SQDIFF 平方差匹配法：该方法采用平方差来进行匹配；最好的匹配值为0；匹配越差，匹配值越大。
    CV_TM_CCORR 相关匹配法：该方法采用乘法操作；数值越大表明匹配程度越好。
    CV_TM_CCOEFF 相关系数匹配法：1表示完美的匹配；-1表示最差的匹配。
    CV_TM_SQDIFF_NORMED 归一化平方差匹配法
    CV_TM_CCORR_NORMED 归一化相关匹配法
    CV_TM_CCOEFF_NORMED 归一化相关系数匹配法

   关于OpenCV的匹配效果如下图所示：

（图片来源：http://www.cnblogs.com/xrwang/archive/2010/02/05/MatchTemplate.html ）

    关于如何使用OpenCV的模板匹配可以参考网址http://www.cnblogs.com/xrwang/archive/2010/02/05/MatchTemplate.html 。OpenCV的目标匹配表面上使用的是前面分析的方式，其实实质上在内部做了很多的处理，比如建立金字塔从顶层粗略计算大概范围，然后再详细计算，以及使用傅里叶变换来进行比较等等（如果不进行这些处理，计算速度超级慢）。接下来抽空会将OpenCV的模板匹配的原理写出来。

上文说到使用OpenCV进行模板匹配的函数matchTemplate，下面就matchTemplate函数的内部处理过程做一个简单的说明。matchTemplate函数的源代码在OpenCV的源代码目录下的 modules/imgproc/src/templmatch.cpp 文件中。其核心函数代码如下（其中的注释是我添加的）：

void matchTemplate( const Mat& _img, const Mat& _templ, Mat& result, int method )
{
    CV_Assert( CV_TM_SQDIFF <= method && method <= CV_TM_CCOEFF_NORMED );

    //numType用来表示模板匹配的方式，0表示相关匹配法，1表示相关系数匹配法，2表示平方差匹配法
    //isNormed表示是否进行归一化处理，true表示进行归一化，false表示不进行归一化处理

    int numType = method == CV_TM_CCORR || method == CV_TM_CCORR_NORMED ? 0 :
                  method == CV_TM_CCOEFF || method == CV_TM_CCOEFF_NORMED ? 1 : 2;
    bool isNormed = method == CV_TM_CCORR_NORMED ||
                    method == CV_TM_SQDIFF_NORMED ||
                    method == CV_TM_CCOEFF_NORMED;

    //判断两幅图像的大小关系，如果输入的原始图像比匹配图像要小，则将原始图像作为模板，原来的模板图像作为搜索图
    Mat img = _img, templ = _templ;
    if( img.rows < templ.rows || img.cols < templ.cols )
        std::swap(img, templ);
    
    CV_Assert( (img.depth() == CV_8U || img.depth() == CV_32F) &&
               img.type() == templ.type() );

   //crossCorr函数是将输入图像做了一次DFT变换（离散傅里叶变换），将空间域的图像转换到频率域中来进行处理，并将处理的结果存放在result中
    int cn = img.channels();
    crossCorr( img, templ, result,
               Size(img.cols - templ.cols + 1, img.rows - templ.rows + 1),
               CV_32F, Point(0,0), 0, 0);
    //如果是相关匹配方法，此处已经计算完毕，返回
 if( method == CV_TM_CCORR )
        return;
    //将模板看作单位1，计算每一个像元所占的百分比（也可以理解为整个模板面积为1，计算每个像元的面积）
    double invArea = 1./((double)templ.rows * templ.cols);

    Mat sum, sqsum;
    Scalar templMean, templSdv;
    double *q0 = 0, *q1 = 0, *q2 = 0, *q3 = 0;
    double templNorm = 0, templSum2 = 0;
    //相关系数匹配算法
    if( method == CV_TM_CCOEFF )
    {
        integral(img, sum, CV_64F);//对原始图像进行求和 templMean = mean(templ);//计算模板图像的均值向量
    }
    else//其他匹配算法 {
        integral(img, sum, sqsum, CV_64F);//计算原始图像的和以及平方和 meanStdDev( templ, templMean, templSdv );//计算模板图像的均值向量和方差向量 templNorm = CV_SQR(templSdv[0]) + CV_SQR(templSdv[1]) +
                    CV_SQR(templSdv[2]) + CV_SQR(templSdv[3]);//计算所有通道的方差和 if( templNorm < DBL_EPSILON && method == CV_TM_CCOEFF_NORMED )
        {//如果所有通道的方差的和等于0，并且使用的方法是归一化相关系数匹配方法，则返回 result = Scalar::all(1);
            return;
        }
        
        templSum2 = templNorm +
                     CV_SQR(templMean[0]) + CV_SQR(templMean[1]) +
                     CV_SQR(templMean[2]) + CV_SQR(templMean[3]);//计算所有通道的均值的平方和 if( numType != 1 )//匹配方式不是相关系数，对模板均值向量和templNorm重新赋值
        {
            templMean = Scalar::all(0);
            templNorm = templSum2;
        }
        
        templSum2 /= invArea;
        templNorm = sqrt(templNorm);
        templNorm /= sqrt(invArea); // care of accuracy here

        q0 = (double*)sqsum.data;
        q1 = q0 + templ.cols*cn;
        q2 = (double*)(sqsum.data + templ.rows*sqsum.step);
        q3 = q2 + templ.cols*cn;
    }
    //下面就是在结果图像中进行查找匹配的结果位置，代码略去，具体可参考OpenCV源代码

    在测试过程中，发现直接用OpenCV打开超过1G（可能比这个大小还要小）的图像就会导致错误，查看OpenCV代码，发现其一次将所有的图像数据全部载入内存，导致内存不够出错，而且OpenCV只支持普通的常用的图像数据格式，并不能直接打开Erdas的img格式或者PCI的pix格式，针对这两个问题，我提出的解决思路是实用GDAL来作为图像数据的读取库，然后分块读取一部分数据（指定一个大小，我指定的是每次读64M的数据），如果图像数据超过该大小，则分块处理，每次处理一块，OpenCV可以建立一个内存数据，就是将图像的RGB值按照特定的顺序读取到内存中，然后将该指针交给OpenCV，OpenCV可以通过该内存数据来构建一个图像，然后再调用模板匹配算法，计算出的结果在还原到原来大的图像中的位置。

    在实际的执行过程中，上述方式是可行的，执行效率大概在一个1G的img图像查找一个128*128的模板需要1分钟左右，相比按照前一篇博客中的方式速度提升了很多倍（1G的数据大概计算了4个小时还没有算完）。对于1min中查找出结果应该是可以接受的，但是感觉还是有点慢，于是想到首先将模板和原始图像同时按照一定的采样比例建立金字塔，然后从金字塔的顶层开始进行匹配，根据顶层找到的位置，再到下一层对应的区域中查找，依次知道查找到原始图像级别。这样就可以减少很多的运算。经过测试，使用这种方式，还是上面的1G的数据，查找只要10S钟即可完成，而且查找的结果也是很准确的。

    在编写程序的过程中，可以使用GDAL的RasteIO函数，直接使用最后面的三个参数，可以直接将图像中的数据读取成OpenCV中的图像存储方式 BGRBGRBGR…

对于OpenCV在进行第一步匹配的时候，是在频率域中进行处理的，相关的函数也在OpenCV的源代码目录下的 modules/imgproc/src/templmatch.cpp 文件中，关于DFT变换（离散傅里叶变换）比较复杂，理论知识尤其是需要很好的数学功底才能够理解，感兴趣的童鞋可以google或者参考相关的数字图像处理的书籍，一般的数字图像处理的书籍中都会有的，因为DFT和FFT在数字图像处理以及信号处理方面使用的非常广泛。

Technorati 标签: OpenCV,模板匹配