基于LK光流金字塔算法与TPS变换的连续时间序列图像配准

LK光流金字塔算法是一种经典的稀疏光流算法，前文我们已经详细讲过其原理与实现：

LK光流金字塔算法原理及C++实现

TPS薄板样条变换是一种经典的非刚性形变模型，前文我们也已经详细讲过其原理、特点与实现：

TPS薄板样条变换计算原理及C++实现

我们知道，基于匹配点对的图像配准，关键在于找到基准图像与浮动图像的多组匹配点对：(p1,p1')、(p2,p2')、...、(pn,pn')。然后使用多组匹配点对计算坐标变换模型（比如仿射变换、TPS变换）。如果使用特征点（点+描述子）匹配的方法寻找匹配点对，需要计算点的描述子，然后根据不同点的描述子的匹配程度来寻找匹配点对，这是比较耗时的。那么，有没有不需要计算描述子就能寻找匹配点对的方法，答案是有的，LK金字塔光流算法就能实现。基于LK金字塔光流算法与TPS变换的图像配准主要步骤如下：

1. 选择一张图像作为固定的参考图像

对其进行直方图均衡化的预处理。直方图均衡化的原理我们前文有讲过，作用是增加图像的对比度和亮度，利于LK金字塔光流算法的跟踪。

数字图像处理之直方图均衡化

直方图均衡化的效果如下图所示：

2. 检测参考图像的角点

首先在参考图像上面检测一些适合光流跟踪的点，要求这些点具有明显的角点特征，且分布均匀，否则不利于光流跟踪或者计算的坐标变换模型不稳定。在数字图像处理中，Harris角点是一种很常用的角点，具有平移和旋转的鲁棒性。后来有研究人员对Harris角点算法进行改进，提出了一种改进后的适合用于光流跟踪的Shi-Tomasi角点算法，在很多情况下其效果比Harris角点算法更好，因此我们选择Shi-Tomasi角点作为跟踪点。

通过使用一个尺寸不变的窗口在图像上滑动，如果在任意方向的滑动都有较大的灰度值变化，则认为该窗口中存在角点，这就是Harris角点检测的基本依据。如图下图所示：

移动窗口时，可将窗口内的灰度变化表示为下式：

其中，[u,v]为窗口偏移量，(x,y)为窗口内的像素坐标，w(x,y)为窗口函数，通常直接取固定值1，有时候也取以窗口中心为原点的二维正态分布。

将上式中的I(x+u,y+v)泰勒展开，由于u，v很小，故把余项忽略掉之后得到：

  

于是有：

上式中M矩阵为：

通常取w(x,y)为固定值1，则有：

                                               

定义角点响应函数：

上式中detM与traceM分别为矩阵M的行列式与迹，k为常数，通常取0.04~0.06。行列式与迹的计算如下式，其中λ₁和λ₂是矩阵M的两个特征值：

判断角点时，首先设定一个阈值，当角点响应函数的值R大于阈值则认为该点为角点。

Shi-Tomasi角点算法的检测过程与Harris角点大同小异，主要区别在于最后的角点判断与剔除，其判断矩阵M的两个特征值λ₁和λ₂中较小的那个是否大于设定的阈值，如果大于则认为该点为角点：

检测到若干Shi-Tomasi角点之后，再对一些不符合要求的角点进行剔除，以确保得到的都是高质量角点，利于后面的光流跟踪：

(1) M矩阵特征值低于qualityLevel.λ_m的特征点被剔除。其中qualityLevel为设定的0~1之间的阈值，该阈值越大被剔除的角点越多，留下的角点质量越好。其中λ_m为检测到的所有角点中最大的M矩阵特征值；

(2) 根据设定的角点之间的最小距离，剔除距离小于该最小距离的点，确保各角点均匀分布，利于后面的图像配准。

对直方图均衡化之后的参考图像进行Shi-Tomasi角点检测，得到p1、p2、...、pn这n个角点检测结果如下图所示：

3. 对参考图像上的角点进行LK光流跟踪获取匹配点对

参考图像是固定的，其它需要配准的多帧连续图像都是浮动图像，使用LK光流算法跟踪参考图像的角点，找到这些角点在浮动图像上的对应位置，从而得到一系列匹配点对，如下图所示：

对于每一个角点的追踪结果，包含一个标志，该标志不为0则表示追踪成功，否则表示追踪失败，将该点及其对应的追踪点剔除。同时，使用RANSAC算法进一步剔除异常匹配点，RANSAC算法的原理我们在前文也有讲过：

图像配准系列之Sift特征点检测与匹配(2)

剔除个别异常匹配点对之后，留下一系列准确的匹配点对：

连续的时间序列图像，是由相机以一定的帧率连续拍摄得到的，比如一个视频文件就是由多帧连续时间序列图像构成的。如果在连续拍摄期间拍摄对象具有运动特征，那么拍摄得到的相邻帧图像的运动具有一定的关联性，比如当前帧的运动是在上一帧运动的基础上进一步运动的结果。对于每一个参考图像上待跟踪的角点，对其进行LK光流跟踪时首先需要确定一个初始跟踪点，在初始跟踪点的基础上进行多轮迭代，从而得到最终跟踪结果。通常初始跟踪点取待跟踪角点的位置，不过我们为了使跟踪更加准确稳定，根据相邻帧的运动关联性我们可以充分利用上一帧的跟踪结果作为当前帧的初始跟踪点。比如对于参考图像上的角点p(x,y)，其在上一帧的跟踪点为p1(x1,y1)，那么当前帧的初始跟踪点p2(x2,y2)我们可以这样取：

x2=x1 + (x1 - x)

y2=y1 + (y1 - y)

4. 计算坐标映射关系

把上一步得到的多组准确匹配点对输入TPS变换模型，得到参考图像与浮动图像的坐标映射关系。

5. 像素重采样

根据上一步得到的坐标映射关系，对浮动图像进行像素重采样，从而使其位置与参考图像对齐。有关像素重采样以及插值问题，前文我们也已经讲过：

常见图像插值算法的原理与C++实现

本文中，我们使用C++和Opencv来实现基于LK光流金字塔算法与TPS变换的图像配准。

我们调用Opencv的函数实现相关功能：

1. 调用equalizeHist函数实现直方图均衡化。

2. 调用goodFeaturesToTrack函数检测Shi-Tomasi角点。

3. 调用calcOpticalFlowPyrLK实现LK光流跟踪。

4. 调用findFundamentalMat实现RANSAC算法剔除异常匹配点。

5. 调用remap函数实现像素重采样。

同时调用前文实现的Tps_TxTy函数实现TPS变换。

TPS薄板样条变换计算原理及C++实现

代码如下：

void lk_test(void)
{
  Mat img0 = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);  //读取参考图像


  Mat img0_e, imgi_e;
  equalizeHist(img0, img0_e);  //对参考图像进行直方图均衡化
  
  vector<Point2f> p1;
  //检测Shi-Tomasi角点
  goodFeaturesToTrack(img0_e, p1, 200, 0.015, 20.0, Mat(), 3, false, 0.04);    //提取shi-tomasi角点
  vector<Point2f> p2 = p1;
  vector<uchar> status;
  vector<float> err;         //跟踪过程中的错误
  TermCriteria termcrit(CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 35, 0.01);
  Size winSize(51, 51);


  for (int i = 1; i <= 100; i++)  //分别配准100张浮动图像
  {
    printf("i=%d\n", i);


    char str[100] = { 0 };
    sprintf(str, "lena/%d.tif", i);
    Mat imgi = imread(str, CV_LOAD_IMAGE_GRAYSCALE);
    equalizeHist(imgi, imgi_e);  //对浮动图像进行直方图均衡化
    
    //跟踪参考图像上角点在浮动图像上的对应位置
    calcOpticalFlowPyrLK(img0_e, imgi_e, p1, p2, status, err, winSize, 3, termcrit, OPTFLOW_USE_INITIAL_FLOW, 0.001);


    vector<uchar> RANSACStatus;//用以标记每一个匹配点的状态，等于0则为外点，等于1则为内点。
    findFundamentalMat(p1, p2, RANSACStatus, CV_FM_RANSAC, 3);//p1 p2必须为float型
    vector<Point2f> f1_features_ok;
    vector<Point2f> f2_features_ok;
    for (int k = 0; k < p1.size(); k++)   //剔除跟踪失败点
    {
      if (status[k] && RANSACStatus[k])  //如果追踪成功、且满足RANSAC条件，则认为该匹配点对使准确的，否则将其剔除
      {
        f1_features_ok.push_back(p1[k]);       //基准图特征点
        f2_features_ok.push_back(p2[k]);     //流动图特征点
      }
    }


    
    Mat out;
    Mat Tx, Ty;
    //将准确的匹配点对输入TPS变换，获取坐标映射
    Tps_TxTy(f1_features_ok, f2_features_ok, imgi.rows, imgi.cols, Tx, Ty);
    remap(imgi, out, Tx, Ty, INTER_CUBIC);  //像素重采样
    
    //使用当前帧的跟踪信息，计算下一帧图像的初始跟踪点
    for (int k = 0; k < p1.size(); k++)
    {
      p2[k] = p2[k] + p2[k] - p1[k];
    }


    Mat diff0 = abs(imgi - img0);
    Mat diff1 = abs(out - img0);


    imshow("imgi", imgi);
    imshow("out", out);
    imshow("imgi-img0", diff0);
    imshow("out-img0", diff1);
    waitKey();
  }
}

运行上述代码，以一张无形变的Lena图像作为参考图像，对100张具有形变的Lena图像进行配准，结果如下所示。

参考图像

待配准的多张Lena图像（浮动图像）

配准之后的多张Lena图像（配准图像）

浮动图像与参考图像的差值图

配准图像与参考图像的差值图

由上述结果可知，配准之后，形变的Lena图像得到了一定的矫正，使其与参考图像位置对齐。然而TPS变换是比较耗时的，实际使用时我们可以使用GPU对其进行并行加速，使其满足实时应用场景。

欢迎扫码关注以下微信公众号，接下来会不定时更新更加精彩的内容噢～