【特征追踪】传统光流方法汇总与介绍（附代码）

最新推荐文章于 2026-07-22 12:29:37 发布

原创最新推荐文章于 2026-07-22 12:29:37 发布 · 1.2k 阅读

18 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

收录于

特征追踪

1 前言

作为连续图像或视频特征追踪的一种经典方法，光流是通过分析图像中像素随时间的位移来捕获场景中物体的运动信息。通过光学流估计获取图像序列中物体的运动轨迹和速度等信息，对于理解视觉场景（前景检测）、运动分析和目标跟踪等任务至关重要。¹

本文将从基本原理、约束条件等方面介绍光流估计的传统算法，并根据算法实现途径的不同，介绍基于梯度、能量、匹配（特征、区域）以及相位等的光流估计的经典算法。同时，由于光流法分为稀疏光流和稠密光流，其区别便是跟踪像素点的多少，一般而言，跟踪几乎所有像素的光流方法才被称为稠密光流，因此该区分在原理部分不再做详细介绍，针对每个方法我将一一列出。接下来，让我们开始吧。

2 传统光流基本原理及约束条件

2.1 基本原理

光流是一种描述像素随时间在图像之间运动的方法，如图所示，随着时间的流逝，同一个像素会在图像中运动，而我们希望追踪它的运动过程。²

像素运动

在传统光流算法中，目标是估计一个视频帧内所有点(稠密光流是所有点，稀疏光流往往是特征点)的位移或速度，这个估计是对所有的点联合执行的。为实现所有像素点的联合估计，光流法必须依赖于某些特定的假设。
传统光流算法的运动只在一个无穷小的距离上被预测，因此，对于大幅度、快速运动的物体光流法存在一定的局限性。其次，虽然可以应用传统光流算法通过时间的推移跟踪目标的移动，但这种方法会以漂移的形式产生累积误差。由于传统光流算法的特定假设，对于图像中的遮挡、运动模糊和无纹理的表面可能表现效果并不出色，因此，相较于近年来基于学习的算法，传统光流算法略显逊色，但其高效的计算效率、较强的可解释性、不同场景的广泛适应能力仍在今天存在着一定的价值。

2.2 光流算法的输出结果

2.2.1 光流场

对于连续的两帧图像，算法输出的光流场描述了图像中每个像素在相邻帧之间的位移。具体来说，它是从一帧图像到另一帧图像之间每个像素的运动向量的集合。

光流场的存储形式是一个二维矩阵，每个元素是一个包含两个分量 $(u, v)$ 的向量，表示水平和垂直方向的位移。我们也可将光流场可视化为彩色图像，其中颜色表示运动方向，亮度表示运动速度。

2.2.2 运动轨迹

对于连续的多帧图像，我们可以通过连接表示同一实际物理意义的像素在不同帧之间的位置从而形成其特定的运动轨迹。
运动轨迹示意图
运动轨迹的保存可以为列表或数组，可将每个特征点的运动轨迹表示在一个列表或数组里，其中每个元素是一个二维坐标 $(x, y)$ ，表示该特征点在某帧中的位置。

2.3 约束条件

根据不同的算法，传统光流的约束条件也各有差别，总得来说，包含以下几点：亮度恒定假设（灰度不变假设）、小运动假设、局部区域一致性假设、全局平滑假设和连续性假设。

2.3.1 亮度恒定假设（灰度不变假设）⭐

亮度恒定假设是指同一场景中的点的灰度值在相邻帧中保持不变，对于 $t$ 时刻位于 $(x, y)$ 处的像素，我们假设它 $t+\Delta t$ 时刻它运动到 $(x + u, y + v)$ 处，由于灰度不变，应有：
$t+\Delta t)$
注意灰度不变假设是一个很强的假设，实际中很可能不成立，事实上，由于物体的材质不同，像素会出现高光和阴影部分；有时，相机会自动调整曝光参数，使得图像整体变亮或变暗。这些时候灰度不变假设都是不成立的，因此光流的结果也不一定可靠。

2.3.2 小运动假设⭐

小运动假设认为在相邻帧之间，像素的位移是很小的。这一假设简化了光流估计的数学模型，使得问题更容易求解。

2.3.3 局部区域一致性假设

局部区域一致性假设是指相邻像素点的运动是相似的。具体来说，将像素在 $x$ 轴和 $y$ 轴上的位移量分别记作 $u, v$ ，在一个小窗口内，我们假定该窗口内的所有像素 $(x, y)$ 都有：
$u(x,y)≈u_0 \\ v(x,y)≈v_0$
这一假设有助于提高光流估计的准确性和稳定性，特别是在具有纹理和边缘的区域。通过考虑局部相似性，可以降低光流场中的噪声和不稳定性。

2.3.4 全局平滑假设

区别于局部区域一致性假设，全局平滑假设表示整个光流场在空间上是平滑的，不会出现突变。全局平滑假设通常通过一个全局的平滑项来实现，平滑项定义如下，是一种在整个图像上的积分：
$E_{smooth}=\alpha \iint (|\nabla u|^2+|\nabla v|^2) {\rm d}x{\rm d}y$
其中， $\nabla u$ 和 $\nabla v$ 代表的是光流场的梯度（空间意义上），由于平滑性，二者均应近似等于 $0$ 。理想的光流场，应在保证预测到全局变化的情况下，使该平滑项最小。
要注意到这是一种全局范围内的假设，意味着相邻像素之间的运动差应该受到限制，光流场不应该有太剧烈的变化。通过引入梯度平方项等正则化项，可以将光流场光滑性的先验知识纳入到优化问题中，有助于避免过拟合和不稳定性。

2.3.5 连续性假设（空间、时间）

连续性假设分为空间连续性假设和时间连续性假设，空间连续性假设的定义等同于全局平滑假设，不过要区分的是，一些算法的空间连续性假设或平滑假设仅限定在窗口内，即只强调局部平滑。
时间连续性假设认为在短时间内，物体的运动是连续的，即光流场在时间上（相邻帧之间）是平滑的，不会出现突变：
$u(x,y,t+\Delta t)≈u(x,y,t)+\Delta u \\ v(x,y,t+\Delta t)≈v(x,y,t)+\Delta v$
其中 $\Delta u$ 和 $\Delta v$ 是小的变化量。
时间连续性假设有助于减少噪声影响，提高光流估计的鲁棒性和准确性。

2.4 图像金字塔——处理大尺度运动

由于小运动假设的存在，在单一尺度上，如果物体的运动较大，传统光流算法将失去鲁棒性，导致光流估计不准确。为解决这一问题，可通过构建⭐图像金字塔⭐，在低分辨率层次上先估计大尺度的运动，然后再逐步细化到高分辨率层次。这样可以逐步累积小的位移，最终得到准确的大尺度运动估计。

应用图像金字塔还可减少图像中的噪声影响，低分辨率层次上的图像噪声较少，可以在这些层次上进行初步的光流估计，然后再逐步细化到高分辨率层次，这样可以减少噪声对最终结果的影响。同时，低分辨率层次上，即使某些区域被遮挡或无纹理，也可以通过其他区域的运动信息来进行补偿，这样可增强算法的鲁棒性。
这一段参考其他博主的介绍³，图像金字塔的可描述如下图：低尺度下找到光流矢量 $d_0$ ，将其扔到下一层去指引我们前行。扔下去后首先要放大两倍，此时你开始抱怨：“上一层找到的 $d_0$ 根本不靠谱，差距好大”（图2中蓝方块与右下角黑色方块的差距）但是你也要庆幸，正是“不靠谱”的 $d_0$ 让你到达了蓝色方块的位置，让你离真实区域前进了很多，否则你还在左上角苦于“小运动”寸步难行呢。在蓝色位置你就满足“小运动”了，继续计算光流，得到 $d_1$ ，然后把 $2d_0+d_1$ 扔到下一层指导我们继续前行，如此往复。

金字塔计算
介绍完了传统光流法的基本原理和假设条件后，我们正式步入具体算法的学习和使用。

3 基于梯度: Lucas-Kanade光流

3.1 原理介绍

Lucas-Kanade光流法，也称为‌LK光流法‌，是一种常用的光流估计算法，由Bruce D. Lucas和Takeo Kanade提出。该算法基于三个基本假设：亮度恒定、小运动和空间一致性。‌
根据亮度恒定假设，有：
$t+\Delta t)$
根据小运动假设，可将上式右侧进行泰勒级数展开，并保留一阶项：
$t+\Delta t)≈I(x,y,t)+\frac{\delta I}{\delta x} u+\frac{\delta I}{\delta y} v+\frac{\delta I}{\delta t} \Delta t$
由于灰度不变，则有：
$\frac{\delta I}{\delta x} u+\frac{\delta I}{\delta y} v+\frac{\delta I}{\delta t} \Delta t = 0$
两边同除以 $\Delta t$ ，得到：
$\frac{\delta I}{\delta x} \frac{u}{\Delta t }+\frac{\delta I}{\delta y} \frac{v}{\Delta t }=-\frac{\delta I}{\delta t}$
式中 $\frac{u}{\Delta t }$ 为像素在 $x$ 轴上的速度， $\frac{v}{\Delta t }$ 为像素在 $y$ 轴上的速度。同时， $\frac{\delta I}{\delta x}$ 为图像在该点处 $x$ 方向的梯度， $\frac{\delta I}{\delta y}$ 则是在 $y$ 方向的梯度， $\frac{\delta I}{\delta t}$ 则是在时间方向上的梯度，分别记作 $I_x, I_y,I_t$ 。因此，可将上式写成如下形式：
$I_x\frac{u}{\Delta t }+I_y\frac{v}{\Delta t }+I_t=0$

到这一步，很多细心的读者可能发现问题了，在很多书籍、文章中，亮度恒定和小位移推导出来的公式明明是 $I_xu+I_yv+I_t=0$ ，而这里为什么是上述的形式，并且，一些书籍还把 $u, v$ 定义为像素在 $x, y$ 方向的速度，而通常我们通常理解的 $u, v$ 是光流场中的两个分量，即 $x, y$ 方向上的位移。这到底是什么情况？🤷‍♀️🤷‍♀️🤷‍♀️
这是因为，通常情况下，为了简化计算，进一步假设了帧间时间 $\Delta t=1$ ，即做了时间尺度的归一化处理。
接下来，公式将以 $I_xu+I_yv+I_t=0$ 的形式进行推导。

将上式写为矩阵形式：
$\left [ \begin{matrix} I_x & I_y \end{matrix} \right ] \left [ \begin{matrix} u \\ v \end{matrix} \right ] = -I_t$
推广到一个 $n \times n$ 的窗口，由于空间一致性假设，则它们共享位移 $u, v$ ，因此则有：
$\left [ \begin{matrix} I_{x1} & I_{y1} \\ I_{x2} & I_{y2} \\ ... \\ I_{xn^{2}} & I_{yn^{2}} \end{matrix} \right ] \left [ \begin{matrix} u \\ v \end{matrix} \right ] = \left [ \begin{matrix} -I_{t1} \\ -I_{t2} \\ ... \\ -I_{tn^{2}} \end{matrix} \right ]$
可作如下简写：
$\left [ \begin{matrix} u \\ v \end{matrix} \right ] = -b$
这是一个关于 $u, v$ 的超定线性方程，因此可用最小二乘解求解，因此有：
$\left [ \begin{matrix} u \\ v \end{matrix} \right ]^* = -(A^TA)^{-1}A^Tb$
由此，我们便求解出了该点处的光流，进而推导到全局。
LK光流通常结合角点检测的算法来选择特征点进行可靠的检测和跟踪，常见的角点检测方法包括Harris角点检测、Shi-Tomasi角点检测等，因此，它是一种稀疏光流。
通过结合几个邻近像素点的信息，LK光流法通常能够消除光流方程里的多义性。而且与逐点计算的方法相比，LK方法对图像噪声不敏感。针对于大幅度运动，LK光流可结合图像金字塔进行逐层求解和优化。

3.2 代码演示

以下是一个使用OpenCV进行LK光流计算的例子：

import cv2
import numpy as np

# 读取视频
cap = cv2.VideoCapture('video.mp4')

# 读取第一帧
ret, prev_frame = cap.read()
prev_gray = cv2.cvtColor(prev_frame, cv2.COLOR_BGR2GRAY)

# 检测特征点(Shi-Tomasi角点检测)
feature_params = dict(maxCorners=100, qualityLevel=0.3, minDistance=7, blockSize=7)
p0 = cv2.goodFeaturesToTrack(prev_gray, mask=None, **feature_params)

# 创建掩码用于绘制轨迹
mask = np.zeros_like(prev_frame)

while True:
    ret, frame = cap.read()
    if not ret:
        break
    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    # 计算光流
    p1, st, err = cv2.calcOpticalFlowPyrLK(prev_gray, gray, p0, None)
    # 选择好的点
    good_new = p1[st==1]
    good_old = p0[st==1]
    # 绘制轨迹
    for i, (new, old) in enumerate(zip(good_new, good_old)):
        a, b = new.ravel()
        c, d = old.ravel()
        a, b = int(a), int(b)
        c, d = int(c), int(d)
        mask = cv2.line(mask, (a, b), (c, d), (0, 255, 0), 2)
        frame = cv2.circle(frame, (a, b), 5, (0, 0, 255), -1)
    img = cv2.add(frame, mask)
    cv2.imshow('Feature Trajectories', img)
    if cv2.waitKey(30) & 0xFF == ord('q'):
        break
    # 更新前一帧和前一点
    prev_gray = gray.copy()
    p0 = good_new.reshape(-1, 1, 2)

cap.release()
cv2.destroyAllWindows()

若要采用金字塔光流，仅需设置cv2.calcOpticalFlowPyrLK()的maxLevel属性即可，此参数为图像金字塔的层数，该参数默认为3，即采用三层金字塔进行计算。

3.3 其他基于梯度法的光流

在基于梯度的光流算法之上，衍生除了许多光流变种，以下是两种常见的光流计算方法，感兴趣的读者可以自行了解。
(1) Gunnar Farneback 方法
Gunnar Farneback 方法是一种高效的稠密光流估计方法，通过多项式展开和频域技术来估计光流场。这种方法首先在局部窗口内对像素值进行多项式拟合，然后利用拟合结果来估计光流向量。Farneback 方法的一个优点是可以快速计算光流，因此在实时应用中非常受欢迎。其具体原理入下图所示。
Farneback方法

(2) Variational Methods 变分方法或能量方法
变分方法通过定义一个目标函数并寻找使该函数达到最小值的光流场来求解光流。这种方法可以灵活地包含各种先验知识，例如光流场的平滑性、亮度恒定假设等。变分方法可以处理复杂的场景，但通常计算成本较高。

读者可能好奇，变分方法已经提及能量函数，为何还划归到梯度方法内？因为基于能量的方法在优化过程中确实使用了梯度信息，例如，欧拉-拉格朗日方程通常涉及对能量函数的一阶导数（梯度）和二阶导数（Hessian）。因此，某种意义上，大家公认的基于能量的光流方法分类部分算法也属于基于梯度的方法，只是在实现时有所不同。接下来，针对基于能量的光流法，我将重点介绍经典的Horn-Schunck光流。

4 基于能量: Horn-Schunck光流

4.1 算法原理

HS光流法是一种全局方法估算光流的算法，属于稠密光流估计方法，它除了要满足LK光流前两个假设（亮度恒定、小运动），还需要满足全局平滑假设⁴。首先，同LK光流法的定义， $u, v$ 分别是像素在 $x, y$ 轴上的位移， $I_x, I_y$ 分别是图像在 $x, y$ 方向的梯度或空间导数，定义一个能量函数，如下：

$E(u,v)=\iint [(I_xu+I_yv+I_t)^2+\alpha ^2(||\nabla u||^2+||\nabla v||^2)] {\rm d}x{\rm d}y$
其中，该能量函数的前半部分 $I_xu+I_yv+I_t)^2$ 的积分结果是灰度变化因子 $E_{data}$ ，后半部分 $\alpha ^2(||\nabla u||^2+||\nabla v||^2)$ 的积分结果为平滑因子 $E_{smooth}$ 。理想的光流场，应使这两项的值最小，即灰度变化小（亮度恒定）并且速度变化小（小运动）。为解决上述泛函极值问题，需引入欧拉-拉格朗日方程进行求解，具体推导见下图。
HS光流法推导

4.2 代码演示

由于OpenCV没有直接提供Horn-Schunck算法的实现，在此，可参考其他博主的代码实现，其中，C++版本来自作者诗眼天涯⁵及原作者Eric Yuan。

#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <math.h>
#include <fstream>
#include <iostream>
 
using namespace cv;
using namespace std;
 
#define ATD at<double>
#define elif else if
 
#ifndef bool
#define bool int
#define false ((bool)0)
#define true  ((bool)1)
#endif
 
 
Mat get_fx(Mat &src1, Mat &src2){
	Mat fx;
	Mat kernel = Mat::ones(2, 2, CV_64FC1);
	kernel.ATD(0, 0) = -1.0;
	kernel.ATD(1, 0) = -1.0;
 
	Mat dst1, dst2;
	filter2D(src1, dst1, -1, kernel);
	filter2D(src2, dst2, -1, kernel);
 
	fx = dst1 + dst2;
	return fx;
}
 
Mat get_fy(Mat &src1, Mat &src2){
	Mat fy;
	Mat kernel = Mat::ones(2, 2, CV_64FC1);
	kernel.ATD(0, 0) = -1.0;
	kernel.ATD(0, 1) = -1.0;
 
	Mat dst1, dst2;
	filter2D(src1, dst1, -1, kernel);
	filter2D(src2, dst2, -1, kernel);
 
	fy = dst1 + dst2;
	return fy;
}
 
Mat get_ft(Mat &src1, Mat &src2){
	Mat ft;
	Mat kernel = Mat::ones(2, 2, CV_64FC1);
	kernel = kernel.mul(-1);
 
	Mat dst1, dst2;
	filter2D(src1, dst1, -1, kernel);
	kernel = kernel.mul(-1);
	filter2D(src2, dst2, -1, kernel);
 
	ft = dst1 + dst2;
	return ft;
}
 
bool isInsideImage(int y, int x, Mat &m){
	int width = m.cols;
	int height = m.rows;
	if (x >= 0 && x < width && y >= 0 && y < height) return true;
	else return false;
}
 
double get_Average4(Mat &m, int y, int x){
	if (x < 0 || x >= m.cols) return 0;
	if (y < 0 || y >= m.rows) return 0;
 
	double val = 0.0;
	int tmp = 0;
	if (isInsideImage(y - 1, x, m)){
		++tmp;
		val += m.ATD(y - 1, x);
	}
	if (isInsideImage(y + 1, x, m)){
		++tmp;
		val += m.ATD(y + 1, x);
	}
	if (isInsideImage(y, x - 1, m)){
		++tmp;
		val += m.ATD(y, x - 1);
	}
	if (isInsideImage(y, x + 1, m)){
		++tmp;
		val += m.ATD(y, x + 1);
	}
	return val / tmp;
}
 
Mat get_Average4_Mat(Mat &m){
	Mat res = Mat::zeros(m.rows, m.cols, CV_64FC1);
	for (int i = 0; i < m.rows; i++){
		for (int j = 0; j < m.cols; j++){
			res.ATD(i, j) = get_Average4(m, i, j);
		}
	}
	return res;
}
 
void saveMat(Mat &M, string s){
	s += ".txt";
	FILE *pOut = fopen(s.c_str(), "w+");
	for (int i = 0; i<M.rows; i++){
		for (int j = 0; j<M.cols; j++){
			fprintf(pOut, "%lf", M.ATD(i, j));
			if (j == M.cols - 1) fprintf(pOut, "\n");
			else fprintf(pOut, " ");
		}
	}
	fclose(pOut);
}
 
void getHornSchunckOpticalFlow(Mat img1, Mat img2){
	double lambda = 0.05;
	Mat u = Mat::zeros(img1.rows, img1.cols, CV_64FC1);
	Mat v = Mat::zeros(img1.rows, img1.cols, CV_64FC1);
 
	Mat fx = get_fx(img1, img2);
	Mat fy = get_fy(img1, img2);
	Mat ft = get_ft(img1, img2);
 
	int i = 0;
	double last = 0.0;
	while (1){
		Mat Uav = get_Average4_Mat(u);
		Mat Vav = get_Average4_Mat(v);
		Mat P = fx.mul(Uav) + fy.mul(Vav) + ft;
		Mat D = fx.mul(fx) + fy.mul(fy) + lambda;
		Mat tmp;
		divide(P, D, tmp);
		Mat utmp, vtmp;
		utmp = Uav - fx.mul(tmp);
		vtmp = Vav - fy.mul(tmp);
		Mat eq = fx.mul(utmp) + fy.mul(vtmp) + ft;
		double thistime = mean(eq)[0];
		cout << "i = " << i << ", mean = " << thistime << endl;
		if (i != 0 && fabs(last) <= fabs(thistime)) break;
		i++;
		last = thistime;
		u = utmp;
		v = vtmp;
	}
	saveMat(u, "U");
	saveMat(v, "V");
 
	imshow("U", u); 
	imshow("v", v);
	waitKey(20000);
}
 
 
 
int main(){
	Mat img1 = imread("table1.jpg", 0);
	Mat img2 = imread("table2.jpg", 0);
	
	img1.convertTo(img1, CV_64FC1, 1.0 / 255, 0);
	img2.convertTo(img2, CV_64FC1, 1.0 / 255, 0);
 
	getHornSchunckOpticalFlow(img1, img2);
	//    waitKey(0);
	return 0;
}

以下是Python版本的HS算法实现：

import numpy as np
import cv2
import matplotlib.pyplot as plt

def horn_schunck(image1, image2, alpha=0.001, iterations=100):
    """
    实现Horn-Schunck光流算法
    :param image1: 第一帧图像 (灰度)
    :param image2: 第二帧图像 (灰度)
    :param alpha: 平滑项权重
    :param iterations: 迭代次数
    :return: 光流场 (u, v)
    """
    # 图像尺寸
    height, width = image1.shape

    # 初始化光流场
    u = np.zeros((height, width))
    v = np.zeros((height, width))

    # 计算梯度
    Ix = cv2.Sobel(image1, cv2.CV_64F, 1, 0, ksize=3)
    Iy = cv2.Sobel(image1, cv2.CV_64F, 0, 1, ksize=3)
    It = image2.astype(np.float64) - image1.astype(np.float64)

    # 迭代求解
    for _ in range(iterations):
        u_avg = (u[1:-1, :-2] + u[1:-1, 2:] + u[:-2, 1:-1] + u[2:, 1:-1]) / 4
        v_avg = (v[1:-1, :-2] + v[1:-1, 2:] + v[:-2, 1:-1] + v[2:, 1:-1]) / 4

        common_term = (Ix[1:-1, 1:-1] * u_avg + Iy[1:-1, 1:-1] * v_avg + It[1:-1, 1:-1]) / \
                      (alpha**2 + Ix[1:-1, 1:-1]**2 + Iy[1:-1, 1:-1]**2)

        u[1:-1, 1:-1] = u_avg - Ix[1:-1, 1:-1] * common_term
        v[1:-1, 1:-1] = v_avg - Iy[1:-1, 1:-1] * common_term

    return u, v

def draw_flow(image, flow, step=16):
    """在图像上绘制光流矢量"""
    h, w = image.shape[:2]
    y, x = np.mgrid[step//2:h:step, step//2:w:step].reshape(2, -1).astype(int)
    fx, fy = flow[y, x].T
    lines = np.vstack([x, y, x + fx, y + fy]).T.reshape(-1, 2, 2)
    lines = np.int32(lines + 0.5)
    vis = cv2.cvtColor(image, cv2.COLOR_GRAY2BGR)
    cv2.polylines(vis, lines, isClosed=False, color=(0, 255, 0), thickness=1)
    for (x1, y1), (_x2, _y2) in lines:
        cv2.circle(vis, (x1, y1), 1, (0, 128, 255), -1)
    return vis

def main():
    # 读取两帧图像
    image1 = cv2.imread('image1.jpg', 0)  # 灰度图像
    image2 = cv2.imread('image2.jpg', 0)  # 灰度图像

    # 计算Horn-Schunck光流
    u, v = horn_schunck(image1, image2, alpha=0.001, iterations=100)

    # 合并光流场
    flow = np.stack((u, v), axis=-1)

    # 绘制光流
    flow_image = draw_flow(image1, flow)

    # 显示结果
    plt.imshow(cv2.cvtColor(flow_image, cv2.COLOR_BGR2RGB))
    plt.title('Horn-Schunck Optical Flow')
    plt.show()

if __name__ == "__main__":
    main()

5 基于匹配的光流

5.1 算法原理

5.1.1 PatchMatch算法

块匹配算法(Patch Match)最早由Barnes于2009年提出，虽然论文中指出的应用场景是立体匹配、纹理合成、图像修复等，但该匹配策略对相关光流算法的产生奠定了理论基础。在块匹配出现之前，为了匹配图像往往使用最原始的暴力匹配，后来的加速策略如KD-tree搜索、PCA降维、局部敏感哈希（LSH)等，虽然提速了很多但仍存在瓶颈。块匹配极大地提升了图像匹配的速率，从而被广泛应用。
要理解基于匹配类型的光流，必须先行理解PatchMatch算法，其中的重中之重更是其⭐匹配传递(propagation)思想⭐，它利用了图像的局部相关性实现了快速的块匹配，通过阅读原始论文及代码⁶和知乎博文⁷，我将PatchMatch简要概述如下。

下图最左侧的两幅图A和B是一对存在横向偏移的图片，对于A中的每个patch都要在B中找到其最匹配的patch，匹配到的patch构成了从A中patch到B中patch的映射。现定义一个最近邻场（NNF: Nearest-Neighbor Field） $f$ 的概念，用于记录所有从A到B的patch的映射，这种映射记录的是目标块（图像B中的patch）相对于参考块（图像A中的patch）的坐标偏移量（offset）。如：某对patch在图像A中的坐标为 $a$ ，在图像B中的坐标为 $b$ ，则 $f (a) = b - a$ 。显然，这个最近邻场 $f$ 携带的映射信息是 $A→\mathbb{R}^2$ ，如果图像A,B完全一致，则这个场处处都是零向量。到这相信不少同学已经觉得这个最近邻场 $f$ 很接近于光流场啦！别急，求解在下面。
PatchMatch
基于最近邻场的快速块匹配算法分为三步，分别对应上图右侧的a、b、c三幅图像，请重点关注图像中的蓝块，并假定它是我们当前的重点匹配块，假设蓝块的位置是 $(x, y)$ ，其左邻居红块的位置是(x-1,y)，上邻居绿块的坐标是 $(x, y - 1)$ ，定义三个匹配块的NNF最近邻场的值（offset值）为 $f (x, y), f (x - 1, y)$ 和 $f (x, y - 1)$ 。算法的三个步骤分别是如下。
(1)初始化——对应图a
初始化是为NNF最近邻场赋值，可以完全赋予随机的初值，也可以加入一些先验信息（如正确匹配关系）作为向导。
(2)迭代——对应图b和c
初始化完成后便开始迭代，每次迭代都是一个全图扫描过程：奇数次迭代是从上到下逐行扫描，每一行从左到右扫描；偶数此迭代反过来，从下到上，从右到左。每个patch被扫描到时，都要先后地经过两个子过程：图b对应的匹配传递（propagation）和图c对应的随机搜索（random search）。
PatchMatching迭代过程
匹配传递的思想就是试图借用邻居的匹配关系，来得到更好的匹配。对于图中的蓝块（假定现在是迭代过程中的当前块，且迭代过程为奇数次迭代），它会在原图(图A)基础上对比其已经完成迭代的左邻居红块和上邻居绿块的最新最近邻场值(offset值)作为自己的最近邻场值，看看哪对匹配最好。换句话说就是，拿左邻居的offset值作为自己的offset值进行匹配对比得到匹配度a，再拿上邻居的offset值作为自己的offset值进行匹配对比得到匹配度b，最后和自己原始offset值的匹配度结果c，三个数值进行对比，取a,b,c中匹配度最大的值对应的offset作为自己的offset。
这种情况下，一旦有patch的匹配效果良好，且处于一致性较好的区域里，那么它会在接下来的扫描中带动其右侧和下侧的邻居们，使它们都得到较好的匹配。因此，在上上图中，假定(a)中红块得到了较好的匹配，在一次迭代中蓝块便可更新到更适合的位置，进而得到了图(b)。
虽然匹配传递可以快速实现匹配，但很容易陷入局部最优问题，为了让其跳出局部最优，需对其施加一个随机搜索策略。
随机搜索策略是以当前块 $(x, y)$ 的匹配块 $(x, y) + f (x, y)$ 为中心，在不断指数衰减的半径区域里随机匹配若干次，直到半径缩小到 $1$ 以下才停止。数学描述如下：
$u_i=v_0+w\alpha ^i \mathrm{R}_i$
其中， $u_i$ 是第i次随机搜索的最近邻场值（offset值）, $v_0$ 是匹配传递后随即搜索前的offset值， $w$ 是最大搜索半径， $\alpha$ 是取值在 $(0, 1)$ 之间的衰减因子， $\mathrm{R}_i$ 是位于 $[- 1, 1] \times [- 1, 1]$ 中服从均匀分布的二维随机数， $i$ 是搜索的次数，直到 $w\alpha ^i$ 小于 $1$ 为止。这个随机搜索的过程可以早停，如刚算到一半就已经比原来更差了，便可提前终止随机搜索。
原文中，作者五次迭代得到的NNF场矩阵见下图。
NNF场

5.1.2 相关光流算法

让我们回到光流本身，许多块匹配的光流算法类似于Patch Matching，它们都有一个共同的流程：

1.划分图像：将原图像和目标图像划分为大小相等的块，如 $N \times N$ 的像素块；
2.选择相似度度量：常用的相似度度量包括均方误差MSE、绝对误差和SAD、归一化互相关NCC等；
3.匹配块初始化及迭代求解匹配块；
4.结果整合：将所有块的运动向量整合成光流场，还可通过插值方法来平滑光流场，使其连续且光滑；
5.后处理：包含滤波、一致性检查等，避免孤立的错误匹配。

其中最经典的算法属BMA(Block Matching Algorithm)，它在搜索方法上做出了改进，如三步法、四步法、钻石搜索法等，为了提高运动估计的准确度，还提出了分层搜索、可变块大小搜索等⁸。除此之外，还有其他一些改进的算法，论文有《DIP: Deep Inverse Patchmatch for High-Resolution Optical Flow》、《PatchMatch Stereo - Stereo Matching with Slanted Support Windows》、《Multiple View Stereo with quadtree-guided priors》等，感兴趣的读者可自行了解查看。

5.2 代码演示

以下是一个简单的代码演示：

import cv2
import numpy as np

def block_matching_optical_flow(frame1, frame2, block_size=16, search_area=15):
    """
    使用块匹配算法计算光流。

    参数:
    frame1 (numpy.ndarray): 第一帧图像 (灰度图)
    frame2 (numpy.ndarray): 第二帧图像 (灰度图)
    block_size (int): 块的大小
    search_area (int): 搜索区域的大小

    返回:
    optical_flow (numpy.ndarray): 光流场 (形状为 (height, width, 2))
    """
    height, width = frame1.shape
    optical_flow = np.zeros((height, width, 2), dtype=np.float32)

    for y in range(0, height - block_size + 1, block_size):
        for x in range(0, width - block_size + 1, block_size):
            # 获取当前块
            block1 = frame1[y:y + block_size, x:x + block_size]
            # 初始化最小误差和最佳匹配位置
            min_error = float('inf')
            best_match = (0, 0)
            # 在搜索区域内寻找最佳匹配块
            for dy in range(-search_area, search_area + 1, 1):
                for dx in range(-search_area, search_area + 1, 1):
                    # 计算搜索块的边界
                    y2_start = max(0, y + dy)
                    y2_end = min(height, y + dy + block_size)
                    x2_start = max(0, x + dx)
                    x2_end = min(width, x + dx + block_size)
                    if y2_end - y2_start != block_size or x2_end - x2_start != block_size:
                        continue
                    # 获取搜索块
                    block2 = frame2[y2_start:y2_end, x2_start:x2_end]
                    # 计算相似度（这里使用SAD）
                    error = np.sum(np.abs(block1 - block2))
                    # 更新最小误差和最佳匹配位置
                    if error < min_error:
                        min_error = error
                        best_match = (dx, dy)
            # 计算运动向量
            motion_vector = np.array(best_match, dtype=np.float32)
            # 将运动向量赋值到光流场
            for i in range(y, y + block_size):
                for j in range(x, x + block_size):
                    optical_flow[i, j] = motion_vector
    
    return optical_flow

def main(f1, f2):
    # 读取两帧图像
    frame1 = cv2.imread(f1, cv2.IMREAD_GRAYSCALE)
    frame2 = cv2.imread(f2, cv2.IMREAD_GRAYSCALE)
    # 计算光流
    optical_flow = block_matching_optical_flow(frame1, frame2, block_size=16, search_area=15)
    # 可视化光流
    flow_image = np.zeros_like(frame1)
    for y in range(0, frame1.shape[0], 16):
        for x in range(0, frame1.shape[1], 16):
            flow_vector = optical_flow[y, x]
            end_point = (int(x + flow_vector[0]), int(y + flow_vector[1]))
            cv2.arrowedLine(flow_image, (x, y), end_point, (255, 255, 255), 1)
    cv2.imshow('Optical Flow', flow_image)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main('1.jpg', '2.jpg')