计算机视觉——OpenCV实现Lucas-Kanade 光流追踪

1.光流

光流法是计算机视觉中用于估计图像序列中物体运动的关键技术。它类似于观察夜空中的彗星，通过其在天空中的运动轨迹来追踪它的路径。在图像处理中，光流帮助我们理解像素点如何在连续的帧之间移动。

1.1 稀疏光流法

稀疏光流法关注于图像中的关键点（通常是角点或显著的特征点），并计算这些点在连续帧中的运动。Lucas-Kanade算法是这种方法的一个经典例子，它通过比较特征点在连续两帧中的灰度值变化来估计这些点的运动。Lucas-Kanade方法适用于跟踪图像序列中的局部运动，尤其是当特征点清晰且显著时。

1.2 稠密光流法

与稀疏光流法不同，稠密光流法计算图像中每个像素的运动，生成一个速度场，其中每个像素都有一个对应的运动向量。Horn-Schunck算法是稠密光流法的一个代表，它通过平滑约束来优化光流场，假设图像亮度在物体运动的方向上变化不大。

1.3 Lucas-Kanade光流

在Lucas-Kanade光流中，图像的灰度值被视为位置和时间的函数。对于一个固定空间点，尽管其在世界坐标系中的位置是固定的，但在图像平面上的像素坐标会随着相机的运动会发生变化。Lucas-Kanade算法通过最小化重投影误差来估计这些像素坐标的变化，即：

$$I(x,y,t)=I(x+\Delta x,y+\Delta y,t+\Delta t)$$

其中，$I(x,y,t)$是在时间$t$的图像中的灰度值，$(dx,dy)$是像素在图像平面上的位移，而$dt$是时间间隔。通过建立一个关于位移$(dx,dy)$的方程组，Lucas-Kanade算法可以估计出特征点的运动。

2. 光流基本假设推导过程

2.1. 光流法的基本假设：
光流法的基本假设是同一个空间点的像素灰度值，在各个图像中的是固定不变的。公式描述为：
$$I(x,y,t)=I(x+\Delta x,y+\Delta y,t+\Delta t)$$

2.2. 泰勒展开：
对上式右侧进行泰勒展开，得到：
$$I(x+\Delta x,y+\Delta y,t+\Delta t)\approx I(x,y,t)+\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t$$

2.3. 光流约束方程：
因为假设了灰度不变，即 $I(x,y,t)$不随 $\Delta x,\Delta y,\Delta t$变化，因此有：
$$\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t=0$$

2.4. 速度和梯度：
将像素在 x 轴上的速度记为 $u$，在 y 轴上的速度记为 $v$，上述方程可以写为：
$$\frac{\partial I}{\partial x}u\Delta t+\frac{\partial I}{\partial y}v\Delta t=-\frac{\partial I}{\partial t}$$

2.5. 矩阵形式：
写成矩阵形式有：
$$\left[\begin{array}{cc}{I}_x& I_y\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=-I_t$$

2.6. 超定线性方程组：
在 Lucas-Kanade 光流中，引入了新的假设作为约束，即某一个窗口内的像素具有相同的运动。考虑一个大小为 $w\times h$ 的窗口，其包含 $k$个像素。因为假设该窗口内像素具有相同的运动，因此可以得到 $k$个方程：
$$\left[\begin{array}{cc}{I}_x^1& I_y^1\\ I_x^2& I_y^2\\ ⋮& ⋮\\ I_x^k& I_y^k\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=-\left[\begin{array}{c}{I}_t^1\\ I_t^2\\ ⋮\\ I_t^k\end{array}\right]$$

记 $A=\left[\begin{array}{cc}{I}_x^1& I_y^1\\ ⋮& \\ I_x^k& I_y^k\end{array}\right]$ 和$b=\left[\begin{array}{c}{I}_t^1\\ ⋮\\ I_t^k\end{array}\right]$，则方程组可以写为：
$$A\left[\begin{array}{c}u\\ v\end{array}\right]=b$$

2.7. 最小二乘解：
该方程是关于 $u$ 和 $v$ 的超定线性方程组，可以使用最小二乘法求解：
$${\left[\begin{array}{c}u\\ v\end{array}\right]}^{\ast }=-(A^{T}A{)}^{-1}{A}^{T}b$$

值得注意的是，上述公式中的 $I_x,I_y,I_t$ 分别表示图像在该点的 x 方向梯度、y 方向梯度和时间方向的梯度。$\Delta x,\Delta y$是像素点在图像平面上的位移，而 $\Delta t$ 是时间间隔。$u$和 $v$是像素点在 x 和 y 方向上的速度。

3.OpenCV中calcOpticalFlowPyrLK函数

OpenCV中calcOpticalFlowPyrLK方法使用迭代Lucas-Kanade算法计算稀疏特征点的光流，用来做特征点跟踪，该方法使用了金字塔，因此具有一定的尺度不变性。
函数原型：

void cv::calcOpticalFlowPyrLK(	
    InputArray 	        prevImg,
    InputArray 	        nextImg,
    InputArray 	        prevPts,
    InputOutputArray 	nextPts,
    OutputArray 	    status,
    OutputArray 	    err,
    Size 	            winSize = Size(21, 21),
    int 	            maxLevel = 3,
    TermCriteria 	    criteria = TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01),
    int 	            flags = 0,
    double          	minEigThreshold = 1e-4 
);

• prevImg:上一帧图像
• nextImg:下一帧图像
• prevPts:上一帧图像中关键点
• nextPts:根据光流计算的上一帧关键点在当前帧中的位置
• status:关键点的跟踪状态，vector,1表示OK,0表示LOST
• err:每个特征点的跟踪误差vector<float>有len(status)==len(nextPts)==len(prePts)==len(err)
• winSize:在每层金字塔中，LK算法中用来求解计算像素运动而假设具有相同运动的窗口大小。
• maxLevel:金字塔的层数，层数多，尺度不变性能更好，运算时间更久
• criteria:迭代搜索算法的终止条件,默认值表示在指定的最大迭代次数criteria.maxCount(30)之后或当搜索窗口移动小于criteria.epsilon(0.01)时终止迭代
• flags:设置误差或者初始值参数,可选下面两个值：

  • OPTFLOW_USE_INITIAL_FLOW设置使用nextPts中的值作为迭代的初始值，如果不设置为OPTFLOW_USE_INITIAL_FLOW,初始状态就使用prevPts中的值，直接从prevPts复制到nextPts,OpenCV源码中对OPTFLOW_USE_INITIAL_FLOW的使用方式为:

    if( flags & OPTFLOW_USE_INITIAL_FLOW )
        nextPt = nextPts[ptidx]*(float)(1./(1 << level));
    else
        nextPt = prevPt;
    

    • OPTFLOW_LK_GET_MIN_EIGENVALS,flags设置为这个值时使用光流运动方程2x2的正规矩阵，也即空间梯度矩阵的最小特征值作为误差项。如果不设置成OPTFLOW_LK_GET_MIN_EIGENVALS,将原始点和移动点周围像素的距离除以窗口中的像素作为误差项。
• minEigThreshold:迭代LK算法会计算光流运动方程2x2的正规矩阵，也即空间梯度矩阵的最小特征值，然后再除以运动不变窗口中的像素总数作为一个误差评价标准，当其小于minEigThreshold时，说明这个点已经追踪不到了，会将其从追踪特征点中移除，避免其对应相素运动的计算，可提升性能。

calcOpticalFlowPyrLK通常和goodFeatureToTrack方法一起使用，先使用GFTTDetector提取特征点的位置，再使用calcOpticalFlowPyrLK追踪其在连续视频流中的位置，避免了特征描述子的计算和特征点的匹配，可以极大的提升追踪的性能。

c++实现：

#include <memory> 
#include <vector>
#include <cstdlib>

#include <opencv2/features2d.hpp>
#include <opencv2/opencv.hpp>

// 定义一个名为TestOpticalFlowLK的类，用于处理光流跟踪
class TestOpticalFlowLK {
    public:
        // 使用std::shared_ptr来定义智能指针，便于管理类的实例
        typedef std::shared_ptr<TestOpticalFlowLK> Ptr;

        // 构造函数
        TestOpticalFlowLK();
        // 默认的析构函数

        // track函数用于处理特征点的跟踪
        void track(std::vector<cv::String> &filenames) const;

    private:
        // 使用cv::Ptr来定义一个智能指针，指向GFTTDetector对象
        cv::Ptr<cv::GFTTDetector> gftt_ptr_;
};

// 实现TestOpticalFlowLK类的构造函数
TestOpticalFlowLK::TestOpticalFlowLK()
{
    // 初始化GFTTDetector对象，用于特征点检测
    gftt_ptr_ = cv::GFTTDetector::create(500, 0.2, 50);
}

// 实现track函数，用于处理特征点的跟踪
void TestOpticalFlowLK::track(std::vector<cv::String> &filenames) const
{
    // 确保filenames至少有一个元素
    assert(filenames.size() > 1);

    // 存储检测到的特征点
    std::vector<cv::KeyPoint> kps1;
    // 存储上一帧和当前帧中特征点的位置
    std::vector<cv::Point2f> pts1, pts2;
    // 存储绘制特征点时使用的颜色
    std::vector<cv::Scalar> colors;

    // 读取第一帧图像，并初始化last_img
    cv::Mat last_img = cv::imread(filenames[0], 0), cur_img;
    // 创建一个掩码，用于绘制特征点
    cv::Mat mask(last_img.size(), CV_8UC1, 255);

    // 使用GFTTDetector检测第一帧图像中的特征点
    gftt_ptr_->detect(last_img, kps1, mask);

    // 遍历所有检测到的特征点
    for(auto &kp : kps1) {
        // 为每个特征点生成一个随机颜色
        int r = (int)(255. * rand() / (RAND_MAX + 1.f));
        int g = (int)(255. * rand() / (RAND_MAX + 1.f));
        int b = (int)(255. * rand() / (RAND_MAX + 1.f));
        // 输出随机颜色的值
        std::cout << "r:" << r << "g:" << g << "b:" << b << std::endl;
        // 将颜色添加到colors数组中
        colors.emplace_back(r, g, b);
        // 将特征点的位置添加到pts1和pts2数组中
        pts1.push_back(kp.pt);
        pts2.push_back(kp.pt);
    }

    // 存储特征点的跟踪状态
    std::vector<uchar> status;
    // 存储每个特征点的跟踪误差
    std::vector<float> err;

    // 将掩码转换为BGR格式，以便绘制
    cv::cvtColor(mask, mask, cv::COLOR_GRAY2BGR);

    // 创建一个Mat对象，用于存储绘制特征点后的图像
    cv::Mat frame;

    // 遍历所有提供的图像文件名
    for (auto &filename : filenames)
    {
        // 输出当前处理的文件名
        std::cout << "filename: " << filename << std::endl;

        // 读取当前帧图像
        cur_img = cv::imread(filename, 0);

        // 使用calcOpticalFlowPyrLK计算特征点的光流
        cv::calcOpticalFlowPyrLK(last_img,
                                 cur_img,
                                 pts1,
                                 pts2,
                                 status,
                                 err,
                                 cv::Size(13, 13),
                                 3,
                                 cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 30, 0.01),
                                 cv::OPTFLOW_USE_INITIAL_FLOW
                                 );

        // 将当前帧图像转换为BGR格式
        cv::cvtColor(cur_img, cur_img, cv::COLOR_GRAY2BGR);

        int cnt = 0; // 用于计数成功跟踪的特征点数量
        // 遍历所有特征点
        for(size_t i = 0; i < status.size(); i++) {
            // 输出当前特征点的跟踪误差
            std::cout << " " << err[i];
            // 如果特征点未成功跟踪，则跳过
            if(!status[i]) continue;

            // 如果特征点的移动距离超过80像素，则跳过
            if(abs((pts1[i].x - pts2[i].x)) > 80 || 
               abs((pts1[i].y - pts2[i].y)) > 80) continue;

            // 在掩码上绘制特征点之间的连线
            cv::line(mask, pts1[i], pts2[i], colors[i], 2);
            // 在当前帧图像上绘制特征点
            cv::circle(cur_img, pts2[i], 10, colors[i], 1);

            // 更新pts1中的坐标为pts2中的坐标
            pts1[i].x = pts2[i].x;
            pts1[i].y = pts2[i].y;

            // 增加成功跟踪的特征点计数
            cnt += 1;
        }
        // 输出每帧成功跟踪的特征点数量
        std::cout << std::endl;

        // 将掩码和当前帧图像混合，以便同时显示原始图像和特征点
        cv::addWeighted(mask, 0.5, cur_img, 0.5, -65, frame);

        // 显示混合后的图像
        cv::imshow("frame", frame);

        // 等待用户按键，0表示无限等待
        cv::waitKey(0);

        // 更新last_img为当前帧图像
        cv::cvtColor(cur_img, cur_img, cv::COLOR_BGR2GRAY);
        last_img = cur_img;

        // 保存混合后的图像到文件
        cv::imwrite("frame.png", frame);
    }
}

Python 代码：

import cv2
import numpy as np
import random

class TestOpticalFlowLK:
    def __init__(self):
        # 初始化GFTTDetector对象，用于特征点检测
        self.gftt_ptr_ = cv2.GFTTDetector_create(500, 0.2, 50)


    def track(self, filenames):
        assert len(filenames) > 1, "At least one filename is required"

        last_img = cv2.imread(filenames[0], 0)
        kps1 = self.gftt_ptr_.detect(last_img)
        pts1 = np.array([kp.pt for kp in kps1], dtype='float32')
        pts2 = np.zeros((len(pts1), 2), dtype='float32')

        colors = [((random.randint(0, 255), random.randint(0, 255), random.randint(0, 255)) for _ in range(len(kps1)))]

        status = np.zeros((len(pts1),), dtype='uint8')
        err = np.zeros((len(pts1),), dtype='float32')

        for filename in filenames[1:]:
            print("Processing file:", filename)
        cur_img = cv2.imread(filename, 0)
        pts2, status, err = cv2.calcOpticalFlowPyrLK(last_img, cur_img, pts1, pts2, None, winSize=(13, 13), maxLevel=3,
                                                     criteria=(cv2.TERM_CRITERIA_COUNT | cv2.TERM_CRITERIA_EPS, 30, 0.01),
                                                     flags=cv2.OPTFLOW_USE_INITIAL_FLOW)

        img = cv2.cvtColor(cur_img, cv2.COLOR_GRAY2BGR)

        for i, (pt1, pt2, stat) in enumerate(zip(pts1, pts2, status)):
            if stat:
        # 确保坐标是一维数组
                pt1 = pt1.reshape(-1)
                pt2 = pt2.reshape(-1)
                cv2.line(img, (int(pt1[0]), int(pt1[1])), (int(pt2[0]), int(pt2[1])), (255), 2)
                cv2.circle(img, (int(pt2[0]), int(pt2[1])), 10, (127), 1)

                pts1 = pts2.copy()
                last_img = cur_img

        cv2.imshow("frame", img)
        cv2.waitKey(0)

        cv2.destroyAllWindows()


if __name__ == "__main__":
    # 假设我们有一个包含图像文件路径的列表
    filenames = ['data/00001.jpg', 'data/00001.jpg', 'data/00002.jpg','data/00003.jpg']  # 以此类推
    flow_tracker = TestOpticalFlowLK()
    flow_tracker.track(filenames)