真实的三维空间中，描述物体运动状态的物理概念是运动场。三维空间中的每一个点，经过某段时间的运动之后会到达一个新的位置，而这个位移过程可以用运动场来描述。

而在计算机视觉的空间中，计算机所接收到的信号往往是二维图片信息。由于缺少了一个维度的信息，所以其不再适用以运动场描述。光流场（optical flow）就是用于描述三维空间中的运动物体表现到二维图像中，所反映出的像素点的运动向量场。

Lucas-Kanade Optical Flow

在计算机视觉中，卢卡斯-金出方法是一种广泛使用的光流估计的差分方法，这个方法是由Bruce D. Lucas和Takeo Kanade发明的。它假设光流在像素点的邻域是一个常数，然后使用最小二乘法对邻域中的所有像素点求解基本的光流方程。 通过结合几个邻近像素点的信息，卢卡斯-金出方法(简称为L-K方法)通常能够消除光流方程里的多义性。而且，与逐点计算的方法相比，L-K方法对图像噪声不敏感。

1. Brightnees Constancy : 就是同一点随着时间的变化，其亮度不会发生改变。这是基本光流法的假定（所有光流法变种都必须满足），用于得到光流法基本方程.
2. Temporal Persistence or Small Movements : 就是时间的变化不会引起位置的剧烈变化，这样灰度才能对位置求偏导（换句话说，小运动情况下我们才能用前后帧之间单位位置变化引起的灰度变化去近似灰度对位置的偏导数），这也是光流法不可或缺的假定.
3. Spatial Coherence : 在一个小的领域内，光流是恒定的。

Brightness Constancy & Temporal Persistence or Small Movements

\[ I(x, y, t) = I(x + \Delta x, y + \Delta y, t + \Delta t) \]

Assuming the movement to be small, the image constraint at\(I(x, y, t)\)with Taylor series can be developed to get:

\[ I(x + dx, y + dy, t + dt) = I(x, y, t) + \frac{\partial I}{\partial x} dx + \frac{\partial I}{\partial y} dy + \frac{\partial I}{\partial t} dt + H.O.T \]

From these equations it follows that:

\[ \frac{\partial I}{\partial x} dx + \frac{\partial I}{\partial y} dy + \frac{\partial I}{\partial t} dt = 0 \]

or:

\[ \frac{\partial I}{\partial x} \frac{dx}{dt} + \frac{\partial I}{\partial y} \frac{dy}{dt} + \frac{\partial I}{\partial t} = 0 \]

which results in:

\[ \frac{\partial I}{\partial x} v_x+ \frac{\partial I}{\partial y} v_y+ \frac{\partial I}{\partial t} = 0 \]

Spatial Coherence

L-K方法假设两个相邻帧的图像内容位移很小，且位移在所研究点p的邻域内为大致为常数。所以，可以假设光流方程 在以p点为中心的窗口内对所有的像素都成立。

\[ \begin{cases} I_x(q_1)v_x + I_y(q_1)v_y = -I_t(q_1) \\ I_x(q_2)v_x + I_y(q_2)v_y = -I_t(q_2) \\ \vdots \\ I_x(q_n)v_x + I_y(q_n)v_y = -I_t(q_n) \\ \end{cases} \]

其中\(q_1, q_2, ... q_n\)为窗口的像素，\(I_x(q_i), I_y(q_i), I_t(q_i)\)是图像在点\(q_i\)和当前时间对位置x, y和时间t的偏导。

\[ A={\begin{bmatrix}I_{x}(q_{1})&I_{y}(q_{1})\\[10pt]I_{x}(q_{2})&I_{y}(q_{2})\\[10pt]\vdots &\vdots \\[10pt]I_{x}(q_{n})&I_{y}(q_{n})\end{bmatrix}},\quad \quad v={\begin{bmatrix}V_{x}\\[10pt]V_{y}\end{bmatrix}},\quad {\mbox{and}}\quad b={\begin{bmatrix}-I_{t}(q_{1})\\[10pt]-I_{t}(q_{2})\\[10pt]\vdots \\[10pt]-I_{t}(q_{n})\end{bmatrix}} \]

\[ Av = b \]

此方程组的等式个数多于未知数个数，所以它通常是超定的。L-K方法使用最小二乘法获得一个近似解.

\[ A^TAv = A^Tb \]

\[ v = (A^TA)^{-1}A^Tb \]

即：

\[ {\begin{bmatrix}V_{x}\\[10pt]V_{y}\end{bmatrix}}={\begin{bmatrix}\sum _{i}I_{x}(q_{i})^{2}&\sum _{i}I_{x}(q_{i})I_{y}(q_{i})\\[10pt]\sum _{i}I_{y}(q_{i})I_{x}(q_{i})&\sum _{i}I_{y}(q_{i})^{2}\end{bmatrix}}^{-1}{\begin{bmatrix}-\sum _{i}I_{x}(q_{i})I_{t}(q_{i})\\[10pt]-\sum _{i}I_{y}(q_{i})I_{t}(q_{i})\end{bmatrix}} \]

经过研究发现，通常靠近图像中心像素的点的贡献要小于较远的点，因此对于上式加上权重。

\[ A^TWAv = A^TWb \]

\[ v = (A^TWA)^{-1}A^TWb \]

\[ {\begin{bmatrix}V_{x}\\[10pt]V_{y}\end{bmatrix}}={\begin{bmatrix}\sum _{i}w_iI_{x}(q_{i})^{2}&\sum _{i}w_iI_{x}(q_{i})I_{y}(q_{i})\\[10pt]\sum _{i}w_iI_{y}(q_{i})I_{x}(q_{i})&\sum _{i}w_iI_{y}(q_{i})^{2}\end{bmatrix}}^{-1}{\begin{bmatrix}-\sum _{i}w_iI_{x}(q_{i})I_{t}(q_{i})\\[10pt]-\sum _{i}w_iI_{y}(q_{i})I_{t}(q_{i})\end{bmatrix}} \]

对于LKT\(hessian = A^TA\)的特征值\(\lambda_1 \, \, \lambda_2\)都很大的时候才能保证空间一致性。

KLT

对于运动速度比较大的情况下，LK光流容易失效，打破了光流的条件约束。因此提出一种金字塔的算法。

金字塔特征跟踪算法描述如下：首先，光流和仿射变换矩阵在最高一层的图像上计算；将上一层的计算结果作为初始值传递给下一层图像，这一层的图像在这个初始值的基础上，计算这一层的光流和仿射变化矩阵；再将这一层的光流和仿射矩阵作为初始值传递给下一层图像，直到传递给最后一层，即原始图像层，这一层计算出来的光流和仿射变换矩阵作为最后的光流和仿射变换矩阵的结果。

算法流程如下：

Demo

/**
******************************************************************************
* @file
* @author  maky <chengwei920412@outlook.com>
* @version
* @date    2019-03-17 19::40::38
* @brief
******************************************************************************
* @attention
*
*
******************************************************************************
*/

#include <iostream>
#include <algorithm>
#include <Eigen/Eigen>
#include <opencv2/opencv.hpp>

class IKanadeLucasTomasi
{
public:
    IKanadeLucasTomasi(){}
    virtual ~IKanadeLucasTomasi(){}
public:
    virtual bool operator()(cv::Mat query, cv::Mat train, const std::vector<cv::Point2f> &query_pts, std::vector<cv::Point2f> &train_pts, std::vector<bool> &error) = 0;
};

class KanadeLucasTomasi
{
public:
    KanadeLucasTomasi(int pyramid_level_total, int window_half_size);
    virtual ~KanadeLucasTomasi();
public:
    bool operator()(cv::Mat query, cv::Mat train, const std::vector<cv::Point2f> &query_pts, std::vector<cv::Point2f> &train_pts, std::vector<bool> &error);
protected:
    IKanadeLucasTomasi *impl_;
};


class KanadeLucasTomasiBase : public IKanadeLucasTomasi
{
public:
    KanadeLucasTomasiBase(){}
    virtual ~KanadeLucasTomasiBase(){}
protected:
    int pyramid_level_total_;
    int window_half_size_;
    std::vector<cv::Mat> query_pyramid_;
    std::vector<cv::Mat> train_pyramid_;
};

class KanadeLucasTomasiCPU : public KanadeLucasTomasiBase
{
public:
    KanadeLucasTomasiCPU(int pyramid_level_total, int window_half_size);
    virtual ~KanadeLucasTomasiCPU();
public:
    virtual bool operator()(cv::Mat query, cv::Mat train, const std::vector<cv::Point2f> &query_pts, std::vector<cv::Point2f> &train_pts, std::vector<bool> &error);
private:
    bool lkOpticalFlow(const cv::Mat &query, const cv::Mat &train, const std::vector<cv::Point2f> &query_pts, const std::vector<cv::Point2f> &guess, std::vector<cv::Point2f> &optical_flows, std::vector<bool> &error);
    bool lkOpticalFlow(const cv::Mat &query, const cv::Mat &train, const cv::Point2f &query_pt, const cv::Point2f &guess, cv::Point2f &optical_flow);
protected:
};

KanadeLucasTomasiCPU::KanadeLucasTomasiCPU(int pyramid_level_total, int window_half_size)
{
    pyramid_level_total_ = pyramid_level_total;
    window_half_size_ = window_half_size;
}
KanadeLucasTomasiCPU::~KanadeLucasTomasiCPU()
{}

bool KanadeLucasTomasiCPU::operator()(cv::Mat query, cv::Mat train, const std::vector<cv::Point2f> &query_pts, std::vector<cv::Point2f> &train_pts, std::vector<bool> &error)
{
    if (query.cols != train.cols || query.rows != train.rows) {
        return false;
    }
    if (query_pts.empty()) {
        return false;
    }
    if (query.empty() || train.empty() || query.size() != train.size()) {
        return false;
    }

    auto width = query.cols;
    auto height = query.rows;
    auto size = query_pts.size();
    train_pts.resize(size, cv::Point2f(0., 0.));
    error.resize(size, false);

    query_pyramid_.resize(pyramid_level_total_);
    train_pyramid_.resize(pyramid_level_total_);
    query_pyramid_[0] = query;
    train_pyramid_[0] = train;
    for (int pos = 1; pos < pyramid_level_total_; pos++) {
        cv::pyrDown(query_pyramid_[pos - 1], query_pyramid_[pos]);
        cv::pyrDown(train_pyramid_[pos - 1], train_pyramid_[pos]);
    }

    auto query_points = query_pts;
    std::vector<cv::Point2f> guess(query_pts.size(), cv::Point2f());
    std::vector<cv::Point2f> optical_flows(query_pts.size(), cv::Point2f());

    for (int pos = pyramid_level_total_ - 1; pos >= 0; pos --) {
        //std::printf("optical flow level: %d.\n", pos);
        //project query_pts to top level of pyramid
        for (int iter = 0; iter < query_points.size(); iter++) {
            query_points[iter] = query_pts[iter] / (1 << (pos));
        }
        if (!lkOpticalFlow(query_pyramid_[pos], train_pyramid_[pos], query_points, guess, optical_flows, error)) {
            return false;
        }
        for (int iter = 0; iter < query_points.size(); iter++) {
            guess[iter] = optical_flows[iter] * 2;
        }
#if 0
        cv::Mat output;
        cv::hconcat(query_pyramid_[pos], train_pyramid_[pos], output);
        cv::cvtColor(output, output, CV_GRAY2RGB);
        cv::Mat query_debug = output.colRange(0, output.cols / 2);
        cv::Mat train_debug = output.colRange(output.cols / 2, output.cols);
        for (auto iter = 0; iter < size; iter++) {
            cv::circle(query_debug, query_points[iter], pyramid_level_total_ - pos, cv::Scalar(255, 0, 0));
            cv::circle(train_debug, cv::Point2f(optical_flows[iter].x + query_points[iter].x, optical_flows[iter].y + query_points[iter].y), pyramid_level_total_ - pos, cv::Scalar(0, 255, 0));
        }
        cv::imshow(cv::format("%d", cv::getTickCount()), output);
        cv::waitKey(1);
#endif
    }
    for (int iter = 0; iter < query_points.size(); iter++) {
        train_pts[iter] = query_pts[iter] + optical_flows[iter];
    }
#if 0
    cv::waitKey(0);
#endif
    return true;
}

bool KanadeLucasTomasiCPU::lkOpticalFlow(const cv::Mat &query, const cv::Mat &train, const std::vector<cv::Point2f> &query_pts, const std::vector<cv::Point2f> &guess, std::vector<cv::Point2f> &optical_flows, std::vector<bool> &error)
{
    if (query_pts.size() != guess.size() || query_pts.size() != optical_flows.size() || query_pts.size() != error.size()) {
        return false;
    }
    auto size = query_pts.size();
    for (auto pos = 0; pos < size; pos++) {
        if (error[pos]) {
            continue;
        }
        error[pos] = !lkOpticalFlow(query, train, query_pts[pos], guess[pos], optical_flows[pos]);
    }
}
bool KanadeLucasTomasiCPU::lkOpticalFlow(const cv::Mat &query, const cv::Mat &train, const cv::Point2f &query_pt, const cv::Point2f &guess, cv::Point2f &optical_flow)
{
    int windows_size = 2 * window_half_size_ + 1;
    Eigen::MatrixXd a = Eigen::MatrixXd::Zero(windows_size * windows_size, 2);
    Eigen::MatrixXd w = Eigen::MatrixXd::Identity(windows_size * windows_size, windows_size * windows_size);
    cv::Mat kernel_x = (cv::Mat_<float>(3, 3) << -3,   0, +3, -10,  0, +10, -3,   0, +3);
    cv::Mat kernel_y = (cv::Mat_<float>(3, 3) << -3, -10, -3,   0,  0,   0, +3, +10, +3);
    //cv::Mat kernel_x = (cv::Mat_<float>(3, 3) << -1, 0, +1, -2, 0, +2, -1, 0, +1);
    //cv::Mat kernel_y = (cv::Mat_<float>(3, 3) << -1, -2, -1, 0, 0, 0, +1, +2, +1);
    int pos = -1;
    for (int y = -window_half_size_; y <= window_half_size_; y++) {
        for (int x = -window_half_size_; x <= window_half_size_; x++) {
            pos++;
            cv::Point point(std::floor(query_pt.x) + x, std::floor(query_pt.y) + y);
            if (point.x + 1 >= query.cols || point.x + 1 < 0 || point.y >= query.rows || point.y < 0) {
                w(pos, pos) = 0;
                continue;
            }
            if (point.x - 1 >= query.cols || point.x - 1 < 0 || point.y >= query.rows || point.y < 0) {
                w(pos, pos) = 0;
                continue;
            }
            if (point.x >= query.cols || point.x < 0 || point.y + 1 >= query.rows || point.y + 1 < 0) {
                w(pos, pos) = 0;
                continue;
            }
            if (point.x >= query.cols || point.x < 0 || point.y - 1 >= query.rows || point.y - 1 < 0) {
                w(pos, pos) = 0;
                continue;
            }
            // scharr
            cv::Mat domain;
            query.rowRange(point.y - 1, point.y + 1 + 1).colRange(point.x - 1, point.x + 1 + 1).convertTo(domain, CV_32F);
            cv::Mat ix = kernel_x.mul(domain);
            cv::Mat iy = kernel_y.mul(domain);
            a(pos, 0) = cv::sum(ix)[0];
            a(pos, 1) = cv::sum(iy)[0];

            //float ix = 0.5 * (query.at<uchar>(point.y, point.x + 1) - query.at<uchar>(point.y, point.x - 1));
            //float iy = 0.5 * (query.at<uchar>(point.y + 1, point.x) - query.at<uchar>(point.y - 1, point.x));
            //a(pos, 0) = ix;
            //a(pos, 1) = iy;
        }
    }
    auto a_t = a.transpose();
    Eigen::Vector2d v = Eigen::Vector2d::Zero();
    for (auto iter = 0; iter < 100; iter++) {
        Eigen::VectorXd b = Eigen::VectorXd::Zero(windows_size * windows_size);
        pos = -1;
        for (int y = -window_half_size_; y <= window_half_size_; y++) {
            for (int x = -window_half_size_; x <= window_half_size_; x++) {
                pos++;
                if (w(pos, pos) < 1e-3) {
                    continue;
                }
                cv::Point point(std::floor(query_pt.x) + x, std::floor(query_pt.y) + y);
                cv::Point prediction(std::floor(query_pt.x + guess.x + v.x()) + x, std::floor(query_pt.y + guess.y + v.y()) + y);
                if (point.x >= query.cols || point.x < 0 || point.y >= query.rows || point.y < 0) {
                    w(pos, pos) = 0;
                    continue;
                }
                if (prediction.x >= train.cols || prediction.x < 0 || prediction.y >= train.rows || prediction.y < 0) {
                    w(pos, pos) = 0;
                    continue;
                }
                float it = query.at<uchar>(point.y, point.x) - train.at<uchar>(prediction.y, prediction.x);
                b[pos] = it;
            }
        }
        auto hessian = (a_t * w * a);
        auto values = hessian.eigenvalues();

        Eigen::Vector2d dv = (hessian).inverse() * a_t * w * b;
        //std::printf("(%8.5f, %8.5f) ", dv.x(), dv.y());
        v = v + dv;
        if (dv.norm() < 1e-4) {
            break;
        }
        if (v.x() > window_half_size_ || v.y() > window_half_size_) {
            return false;
        }
    }
    optical_flow = cv::Point2f(guess.x + v.x(), guess.y + v.y());
    //std::printf("query: (%8.5f, %8.5f), guess: (%8.5f, %8.5f), train: (%8.5f, %8.5f), v: (%8.5f, %8.5f)\n", query_pt.x, query_pt.y, query_pt.x + guess.x, query_pt.x + guess.x, query_pt.x + optical_flow.x, query_pt.y + optical_flow.y, v.x(), v.y());
    return true;
}


KanadeLucasTomasi::KanadeLucasTomasi(int pyramid_level_total, int window_size)
#ifdef KLT_GPU
    : impl_(new KanadeLucasTomasiGPU(pyramid_level_total, window_size))
#else
    : impl_(new KanadeLucasTomasiCPU(pyramid_level_total, window_size))
#endif
{}
KanadeLucasTomasi::~KanadeLucasTomasi()
{}

bool KanadeLucasTomasi::operator()(cv::Mat query, cv::Mat train, const std::vector<cv::Point2f> &query_pts, std::vector<cv::Point2f> &train_pts, std::vector<bool> &error)
{
    return impl_->operator()(query, train, query_pts, train_pts, error);
}

int main(int argc, const char **argv)
{
    cv::Mat left = cv::imread("./res/1403715283362142976.png");
    cv::Mat right = cv::imread("./res/1403715283412143104.png");
    cv::Mat output;
    cv::hconcat(left, right, output);
    cv::Mat left_debug = output.colRange(0, output.cols / 2);
    cv::Mat right_debug = output.colRange(output.cols / 2, output.cols);

    cv::cvtColor(left, left, CV_RGB2GRAY);
    cv::cvtColor(right, right, CV_RGB2GRAY);

    std::vector<cv::Point2f> left_pts, right_pts, feedback_pts;
    std::vector<uchar> status;
    cv::goodFeaturesToTrack(left, left_pts, 500, 0.01, 20, cv::Mat(), 10);
    cv::cornerSubPix(left, left_pts, cv::Size(5, 5), cv::Size(-1, -1), cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 30, 0.001));
#if 0
    {
        std::vector<float> error;
        cv::calcOpticalFlowPyrLK(left, right, left_pts, right_pts, status, error, cv::Size(31, 31), 3, cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 30, 0.001), 0, 0.001);
        cv::calcOpticalFlowPyrLK(right, left, right_pts, feedback_pts, status, error, cv::Size(31, 31), 3, cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 30, 0.001), 0, 0.001);
    }
#else
    {
        std::vector<bool> error;
        KanadeLucasTomasi tracker(3, 15);
        tracker(left, right, left_pts, right_pts, error);
        tracker(right, left, right_pts, feedback_pts, error);
        status.resize(left_pts.size());
        for (auto pos = 0; pos < status.size(); pos ++) {
            if (error.at(pos)) {
                status.at(pos) = 0;
            }
            else {
                status.at(pos) = 0x1;
            }
        }
    }
#endif

    if (left_pts.size() != feedback_pts.size()) {
        return -1;
    }
    auto size = left_pts.size();
    for (auto pos = 0; pos < size; pos++) {
        if (!status.at(pos)) {
            continue;
        }
        auto &point = feedback_pts.at(pos);
        auto &left_pt = left_pts.at(pos);
        if (cv::norm(point - left_pt) > 3.) {
            //continue;
        }

        cv::circle(left_debug, left_pts.at(pos), 2, cv::Scalar(0, 255, 0), -1);
        cv::circle(right_debug, right_pts.at(pos), 2, cv::Scalar(0, 255, 0), -1);

        cv::line(left_debug, cv::Point(point.x - 4, point.y), cv::Point(point.x + 4, point.y), cv::Scalar(255, 0, 0), 1);
        cv::line(left_debug, cv::Point(point.x, point.y - 4), cv::Point(point.x, point.y + 4), cv::Scalar(255, 0, 0), 1);

        cv::circle(right_debug, left_pts.at(pos), 1, cv::Scalar(0, 0, 255), -1);
        cv::line(right_debug, left_pts.at(pos), right_pts.at(pos), cv::Scalar(0, 0, 155), 1);
        std::printf("(%8.5f, %8.5f) --> (%8.5f, %8.5f)\n", left_pts.at(pos).x, left_pts.at(pos).y, right_pts.at(pos).x, right_pts.at(pos).y);
    }
    cv::imshow("klt tracker", output);
    cv::waitKey();
    return 0;
}

---

Reference

Pyramidal Implementation of the Lucas Kanade Feature Tracker Description of the algorithm

https://en.wikipedia.org/wiki/Lucas%E2%80%93Kanade_method

http://cecas.clemson.edu/~stb/klt/

https://homes.cs.washington.edu/~shapiro/EE596/notes/Optical_Flow.pdf

http://www.mpig.com.cn/

http://www.cs.ucf.edu/~gvaca/REU2013/p4_opticalFlow.pdf

http://eric-yuan.me/coarse-to-fine-optical-flow/

https://zh.m.wikipedia.org/wiki/%E7%B4%A2%E8%B2%9D%E7%88%BE%E7%AE%97%E5%AD%90

http://www.cnblogs.com/lancidie/archive/2011/07/17/2108885.html

Lucas-Kanade 20 Years On: A Unifying Framework