简介：这份凸优化大作业代码包面向学习运动估计与计算机视觉课程的学生，围绕MV子视角的帧间预测任务，给出十一种可运行的实现思路：运动估计部分包含穷尽块匹配、三步搜索与分层块匹配，另设光流估计和单应性矩阵两类方法，并支持三大类方法相互组合，便于对照不同策略的预测效果与参数设置。压缩包共300个文件、约214.99MB，其中216个png多为实验序列与结果图，60个m文件承载MATLAB主脚本和函数，另有mat数据、py辅助脚本及txt说明，目录可支撑逐项复现与二次修改。已有301人学习，适合课程大作业参考、算法对比实验和帧间预测入门实践，也能帮助读者理解运动估计、光流与单应性变换在视频预测中的组合方式。

1. 帧间预测的真相：十一种方法最后都在解同一个能量泛函

把同一对相邻帧输入十一种帧间预测实现，最容易产生的错觉是“这是一组十一种彼此独立的算法”。可把输出摊开看：块匹配给出的是一张整数运动矢量图，Lucas-Kanade 给出的是稀疏角点位移，Horn-Schunck 给出的是稠密流场，单应性配准只给一个 3×3 矩阵，相位相关给的是一个亚像素平移量。形态差别极大，但真正决定精度的关键因素翻来覆去只有两个：运动模型假设了什么形状，正则项假设了什么先验。这份凸优化大作业里的十一份 MATLAB 代码，价值不在于“多”，而在于它把“模型复杂度”这一维展开了，让你能亲眼看到：同一段序列，换运动模型带来的误差跳变，往往比换求解器大一个数量级。搞懂这一维，比背十一个函数名有用得多。

2. 运动模型的分层：从平移矢量到八自由度单应性，多一个参数意味着什么

2.1 块匹配的代价函数：SAD 与 SSD 背后是两种不同的噪声假设

块匹配是帧间预测里最“工程”的一类方法，但它并不粗糙。假设当前帧的一个块在参考帧里只是整体平移了(dy, dx)，那么匹配就是在搜索窗内找一个位移让两块“最像”。而“最像”的定义，直接决定了这个估计在什么噪声下是极大似然解：用绝对差之和（SAD）时，你隐含假设残差服从拉普拉斯分布，尾部厚、对离群点不敏感；用平方差之和（SSD）时，隐含假设残差服从高斯分布，对中等误差敏感，但对一个坏点会付出平方级代价。

function [mv, cost] = block_match_sad(ref, cur, blk, search) % ref/cur: 单通道 double，已归一化到 [0,1] % blk: 块边长（奇数）；search: 搜索半径 [H, W] = size(cur); r = (blk - 1) / 2; mv = zeros(H, W, 2); cost = inf(H, W); pad = padarray(ref, [search + r, search + r], 'replicate'); for y = 1 + r : H - r for x = 1 + r : W - r tpl = cur(y - r : y + r, x - r : x + r); % 当前帧块作模板 best = inf; bv = [0 0]; for dy = -search : search for dx = -search : search cand = pad(y + search + dy : y + search + dy + blk - 1, ... x + search + dx : x + search + dx + blk - 1); c = sum(abs(cand(:) - tpl(:))); % L1 数据项，即 SAD if c < best, best = c; bv = [dy dx]; end end end mv(y, x, :) = bv; cost(y, x) = best; end end end

这段全搜索的复杂度是O(W·H·(2S+1)²·blk²)，S=16、blk=16时单帧就要上亿次运算，而这正是各种快速搜索算法出现的动机：三步搜索、菱形搜索、六边形搜索，本质都是在代价面上做贪心下降。这里有一个必须说清楚的取舍：全搜索在整数格点上做的是穷举，得到的是全局最小，但那个“最小”只是离散格点上的最小；快速搜索连这个保证都放弃了。代价面在纹理丰富区域通常是单峰的，贪心还能收敛到附近；一旦进入大面积平坦区或周期纹理（栅栏、条纹、键盘），代价面会出现多个几乎等高的极小值，贪心搜索会随机落入其中一个，这就是后面避坑章节要专门处理的问题。

另一个常被忽略的点是亚像素。整数格点最优只说明真实位移落在某个半像素范围内，直接用整数矢量做补偿，残差里会残留系统性偏差。标准的修补方式是在最优整数点周围对代价面做一维或二维抛物线拟合：

% 在最优整数位移 (dy, dx) 处做可分离抛物线插值 c0 = cost(dy, dx); cx1 = cost(dy, dx + 1); cx0 = cost(dy, dx - 1); sub_x = 0.5 * (cx0 - cx1) / (cx0 - 2*c0 + cx1); % 分母加 eps 防零除

sub_x的物理含义是“真实极值点相对整数格点的偏移量”，取值应在[-0.5, 0.5]。如果算出来超出这个范围，说明代价面在该邻域不是抛物线形状，通常意味着落进了多峰区或者块内出现了遮挡，此时更稳妥的做法是把亚像素结果截断回 ±0.5，而不是硬信这个拟合值。

2.2 单应性矩阵与 DLT：为什么必须先做坐标归一化

块匹配只能表达平移，仿射变换能表达平移、旋转、缩放、错切共六个自由度，单应性（Homography）则把平面到平面的透视关系也包含进来，共八个自由度。这套代码里出现单应性矩阵，通常是为了处理相机存在旋转或平面场景透视变化的情况。求解它的经典办法是直接线性变换（DLT）：把x' = H x写成齐次坐标下的叉乘形式x' × Hx = 0，每对匹配点贡献两个线性独立的方程，四对点凑够八个方程，解一个 9 维齐次线性方程组，取系数矩阵最小奇异值对应的右奇异向量。

function H = dlt_homography(p1, p2) % p1, p2: 3xN 齐次坐标点，第 1、2 行是 x、y [T1, n1] = normalize_pts(p1); [T2, n2] = normalize_pts(p2); A = zeros(2 * size(n1, 2), 9); for i = 1 : size(n1, 2) x = n1(1,i); y = n1(2,i); X = n2(1,i); Y = n2(2,i); A(2*i-1, :) = [0 0 0 -x -y -1 Y*x Y*y Y]; % 叉乘展开的第一行 A(2*i, :) = [x y 1 0 0 0 -X*x -X*y -X]; % 第二行 end [~, ~, V] = svd(A); Hn = reshape(V(:, end), 3, 3).'; % 注意列优先，必须转置 H = T2 \ Hn * T1; % 反归一化回原坐标系 H = H / H(3, 3); % 固定尺度 end function [T, pn] = normalize_pts(p) % 各向同性归一化：平移到质心，再缩放到平均距离 sqrt(2) c = mean(p(1:2, :), 2); d = mean(sqrt(sum((p(1:2, :) - c).^2, 1))); s = sqrt(2) / d; T = [s 0 -s*c(1); 0 s -s*c(2); 0 0 1]; pn = T * p; end

normalize_pts这十几行是整个 DLT 里最容易被删掉、又最不能删的部分。原因在数值上：图像坐标动辄几百上千，九个未知数里有常数项、有一次项、有二次项，量级相差三四个数量级，A的条件数会大到双精度都撑不住，SVD 给出的最小奇异向量会被舍入误差主导，解出来的H看上去能用，实际重投影误差大得离谱。Hartley 归一化把两组点分别搬到质心在原点、平均半径 √2 的位置，让A各列量级接近，条件数下降几个数量级，这一步是血泪经验换来的，不是可选项。

还有一处细节值得单独提：reshape(V(:,end), 3, 3).'里的转置。MATLAB 是列优先，reshape按列填充，而H我们习惯按行理解，少写这个转置，解出来的矩阵就是转置的H，表现出来的现象是“图像被沿对角线方向拉歪”，很多人会误以为是点对顺序错了。

2.3 模型自由度与可用数据量的匹配：单应性不是越强越好

把三种模型放在一起看自由度：平移 2、相似 4、仿射 6、单应性 8。每增加一个自由度，就需要更多的独立观测来约束它，而且这些观测必须在几何上“张满”参数空间。四个点解单应性，前提是任意三点不共线；如果四个点几乎落在一条直线上，A矩阵接近秩亏，H的解会在某个方向上剧烈放大，重投影看起来还行，但只要场景稍微偏离这四个点的凸包，误差就会爆炸。

更实际的问题是尺度。在小块区域上估计单应性，等于用很少的像素去拟合一个能表达透视的模型，过拟合几乎必然发生：H会把匹配噪声也一并拟合进去，导致变换后的图像在块内残差极小，块外完全崩掉。所以在这套代码里，单应性通常只用在“整帧全局配准”这条路径上，局部运动仍然交给块匹配或光流。判断标准其实很朴素：当你手里的匹配点对数量只是自由度的两三倍时，优先降模型，而不是加正则。仿射在多数场景下已经能吃掉旋转和缩放，剩下的透视畸变，除非相机明显倾斜，否则收益经常小于它引入的方差。

3. 光流法的凸优化内核：亮度恒常假设如何变成一个可解的最小二乘

3.1 光流约束方程的泰勒展开与孔径问题

光流法的出发点只有一个假设：同一个物理点在相邻两帧里亮度不变。写成式子就是I(x, y, t) = I(x+u, y+v, t+1)，其中(u, v)是这一点的位移。对右端在(x, y)处做一阶泰勒展开并丢掉高阶项，就得到光流约束方程：

I_x · u + I_y · v + I_t = 0

I_x、I_y是空间梯度，I_t是时间差分，全部可以从两帧图像直接算出来，唯独u和v是未知的。一个方程、两个未知数，解空间是一条直线而非一个点——这就是经典的孔径问题。几何上它的含义很直白：单个像素只能观测到沿梯度方向的变化，垂直于等亮度线的那一个分量（法向流）是可测的，沿等亮度线方向的分量完全不可测。

这一步的推导还藏着第二个坑：泰勒展开成立的前提是位移足够小，小到一阶近似能覆盖真实的亮度变化。经验阈值是 1 到 2 个像素，位移超过这个量级，I_t里混入的就不再是线性项而是高阶项，“解出来的流”会系统性地小于真实位移。这个偏差不是靠调参数能修掉的，它必须靠 3.3 节的分层策略解决。

3.2 Lucas-Kanade 与 Horn-Schunck：局部常数与全局平滑的对偶

补上缺失约束的两条路，正好构成了这套代码里两个标志性方法的分野。Lucas-Kanade 的思路是“在一个小窗口内，所有像素共享同一个(u, v)”，把窗口内几十上百个方程堆起来做加权最小二乘，正规方程是：

[Σ w I_x² Σ w I_x I_y] [u] [Σ w I_x I_t] [Σ w I_x I_y Σ w I_y²] [v] = -[Σ w I_y I_t]

左边这个 2×2 矩阵就是结构张量，它的两个特征值决定了这个窗口“能不能定住”位移。两个特征值都大，说明窗口内有角点，两个方向都有梯度，可解；一个大一个小，说明是边缘，只能定住法向；两个都小，说明是平坦区，矩阵接近奇异，解出来的值本质上是噪声放大。

function [u, v, conf] = lk_flow(I0, I1, win) % I0/I1: 已高斯平滑的 double 灰度图；win: 窗口边长 [Ix, Iy] = gradient(I0); It = I1 - I0; % 时间差分，注意符号方向 k = ones(win) / (win^2); Sxx = conv2(Ix.*Ix, k, 'same'); Sxy = conv2(Ix.*Iy, k, 'same'); Syy = conv2(Iy.*Iy, k, 'same'); bx = conv2(Ix.*It, k, 'same'); by = conv2(Iy.*It, k, 'same'); det = Sxx.*Syy - Sxy.^2 + eps; % eps 保证平坦区不产生 Inf u = -( Syy.*bx - Sxy.*by) ./ det; v = -(-Sxy.*bx + Sxx.*by) ./ det; lam_min = 0.5*(Sxx + Syy) - sqrt(0.25*(Sxx - Syy).^2 + Sxy.^2); conf = lam_min > 1e-4; % 最小特征值门限，标记可信像素 u(~conf) = 0; v(~conf) = 0; end

conf这个掩膜是 LK 能不能用的关键，很多人把这一步省了，直接拿整张流场去算端点误差，结果被平坦区的天文数字污染，指标难看还以为算法不行。Shi-Tomasi 用最小特征值做角点判据，用的就是同一个量。

Horn-Schunck 走的是另一条路：不假设局部常数，而是假设整个流场是光滑的，把问题写成一个全局能量泛函：

E(u, v) = ∫∫ (I_x u + I_y v + I_t)² dxdy + α ∫∫ (|∇u|² + |∇v|²) dxdy

前一项是数据项，惩罚违反亮度恒常；后一项是正则项，惩罚流场剧烈变化。对它变分求极值会得到一对耦合的欧拉-拉格朗日方程，离散后可以用 Jacobi 或 Gauss-Seidel 迭代求解，每次迭代用邻域均值作为平滑估计。把两个方法对照着看，会发现它们其实是同一个泛函的两种极端正则：LK 的正则是“窗口内等于常数”的硬约束，HS 的正则是“梯度平方要小”的软约束。这个视角很有用——后面要讲的 TV-L1，无非是把 HS 的 L2 正则换成 L1 正则、把 L2 数据项换成 L1 数据项。

3.3 金字塔分层：绕开非凸性的标准做法

回到 3.1 节那个“位移必须小于 1 到 2 像素”的硬约束。真实视频里几十像素的位移随处可见，于是所有能用在实际场景里的光流实现，都会在时间维度之外再加一个尺度维度：把两帧都建成图像金字塔，在最顶层（最小分辨率）先估一个粗糙的流，上采样后作为下一层的初值，把图像按当前流场 warp 回去，再在残差上估增量。位移在每降一层就减半，顶层哪怕只有 40 像素的位移，到最底层也变成 10 像素以内，泰勒展开就能站得住。

从优化角度看，这一步的作用是绕开非凸性，而不是加速。原始能量泛函在位移较大时是非凸的，直接迭代很容易停在错误的局部极小；金字塔提供了一条“由粗到细”的同伦路径，粗尺度上的代价面被平滑过，极小值少而宽，先落进正确的盆地，再逐层细化。层数的经验取法是log2(max_displacement)，四层能覆盖大约 16 像素，五层覆盖 32 像素。

warp 的实现细节值得单独说：必须用双线性或双三次插值，且 warp 的方向要和流场的定义方向一致。interp2(I1, X + v, Y + u)表示“从I1的(x+v, y+u)处采样放到(x, y)”，如果写成X - v而流场又按前向定义，误差会被放大一倍，表现出来的现象是流场幅值系统性偏小、边缘处明显错位。这个符号错误在初学者代码里出现频率高得惊人。

4. 正则项与求解器：TV-L1 为什么必须靠 ADMM 或原始-对偶拆开

4.1 二次正则的过度平滑与全变差的边缘保持

HS 用的 L2 正则有一个几何后果：它对应各向同性扩散，扩散系数是常数，于是无论边缘两侧亮度差多大，平滑力度都一样。结果是运动边界（前景物体和背景之间）被抹成一条渐变带，物体轮廓附近的流场既不像前景也不像背景。数据项用 L2 也有类似问题：一个被遮挡的像素会产生巨大的亮度差，平方之后在总能量里占压倒性权重，求解器会不惜扭曲整个邻域的流场去迁就它。

全变差（TV）正则把|∇u|²换成|∇u|，效果是“代价随梯度线性增长”。在平坦区梯度小，TV 和 L2 差不多；在边缘处梯度大，L2 的代价按平方爆炸，TV 只是线性增长，于是优化器宁愿保留一个跳变也不愿意把它抹平。这就是 TV 能保边的机理，代价是它会把缓慢的亮度斜坡变成一段段平台，出现所谓的“阶梯效应”。数据项换成 L1 则是为了鲁棒：|ρ|对大的残差只线性惩罚，遮挡区、光照突变区的像素不再主导整个解。

需要提前打个预防针：TV-L1 的能量泛函是不可微的，|·|在零点没有导数，|∇u|在梯度为零处也不可微。这意味着梯度下降那一套直接失效，必须换求解器。

4.2 ADMM 的增广拉格朗日与软阈值闭式解

ADMM 的思路是把一个难解的问题拆成两个各自好解的子问题。对min f(x) + g(z) s.t. x = z这种形式，交替做三件事：固定z优化x（通常是一个二次问题或梯度步）、固定x优化z（通常有闭式解，叫 prox 算子）、更新对偶变量。对于 TV-L1，把梯度算子∇u单独拆成一个变量，TV 项就变成“约束该变量落在 L1 球内”，其 prox 就是逐点的软阈值或球面投影，一次算术就完事。

function [u, v] = tvl1_flow_step(I0, I1, u, v, lambda, tau, sigma, iters) % 单层 TV-L1 光流：原始-对偶交替，结构示意 % 注意 grad/div 必须互为负伴随，否则迭代会缓慢漂移 [Ix, Iy] = gradient(I0); px = zeros(size(I0)); py = zeros(size(I0)); % TV 项的对偶变量 for k = 1 : iters % (1) 按当前流场 warp，得到亮度残差 I1w = interp2(I1, (1:size(I1,2)) + v, (1:size(I1,1))' + u, 'linear', 0); rho = I1w - I0; % (2) 数据项 L1 的软阈值：步长 lambda*tau，逐点 sign u = u - tau * lambda * Ix .* sign(rho); v = v - tau * lambda * Iy .* sign(rho); % (3) TV 项的对偶上升，再逐点投影到单位球 px = px + sigma * grad_x(u); py = py + sigma * grad_y(v); nrm = max(1, sqrt(px.^2 + py.^2)); px = px ./ nrm; py = py ./ nrm; % (4) 对偶变量回写原始变量 u = u + tau * div_x(px); v = v + tau * div_y(py); end end

sign(rho)这一步就是 L1 数据项的全部内容：残差为正就往负梯度方向推lambda*tau，为负就往正方向推，残差恰好为零则不动。这就是著名的软阈值收缩，它天然产生“稀疏残差”——大部分像素的残差被压到零，只有遮挡和真实运动不连续处保留非零值，这正好符合我们对帧间预测残差的预期。TV 项那边的px ./ max(1, |p|)是投影到单位球，保证对偶变量满足 TV 的对偶约束|p| ≤ 1。

这里的tau和sigma不是随便填的，Chambolle-Pock 型算法的收敛条件是tau * sigma * ||K||² < 1，K是梯度算子，其谱范数上界是√8。所以常用的安全取值是tau = sigma = 1/√8 ≈ 0.35，实际代码里常取 0.2 留余量。步长设大了，现象是能量震荡不降，甚至几步之内就溢出成 NaN。

4.3 数值实现细节：warp 方向、步长与收敛判据

外层还要套一层对u、v本身的固定点迭代，因为上面的线性化只在当前流场附近成立。典型写法是每轮用当前(u, v)warp 一次，内层跑 20 到 30 次原始-对偶，再进入下一轮。外层轮数不必太多，5 到 10 轮基本收敛，因为位移在金字塔的每一层已经被压得很小。

收敛判据不能用“能量变化小于阈值”这一条，因为能量里各项量级差很多，L1 项在遮挡区贡献巨大，掩盖了流场本身的变化。更稳的做法是同时看两个量：一是平均端点误差的逐轮变化，二是|u_new - u_old|的最大值。前者反映精度是否还在提升，后者反映解是否已经稳定。任一个指标连续几轮不动，就可以停。

边界处理也容易被忽略。梯度算子在一像素宽的边界上只能做单边差分，量级和内部不一致，流场在边界会往外“甩”。常见做法是把边缘几行几列的结果复制或裁剪掉，或者直接在计算前把图像 edge-pad 一圈再算，算完裁回原尺寸。

5. 避坑与排查：十一种方法跑同一段序列时最容易翻车的五处细节

5.1 光流场整体接近纯黑或纯白

现象最典型：quiver画出来几乎看不到箭头，或者用色轮显示时整张图落在同一颜色上。原因排第一的是数据类型和量级不匹配——图像如果是uint8，梯度值在 0 到 255 之间，It也在同一量级，直接代进 LK 的正规方程，解出来的位移量级只有零点几甚至更小；如果是 HSV 那类范围更窄的输入，问题会更严重。解决办法是在入口处统一转成double并归一化到[0,1]，再做高斯平滑（sigma取 1 到 1.5，太小压不住噪声，太大把细节 blur 掉）。

第二个原因是结构张量奇异，平坦区det ≈ 0，除法把噪声放大到1e6量级，色轮显示时这些值被截断，整张图看上去就是糊的。解决方案就是 3.2 节里那个conf掩膜，用最小特征值门限把不可信像素剔掉，图上只保留可信区域的箭头。第三个原因更隐蔽：It = I1 - I0的符号方向。如果把It写成I0 - I1，解出来的流场整体反向，quiver画出来箭头全指向反方向，却又不至于完全看不见，很容易被当成“算法精度差”而不是“符号写反”。

5.2 单应性矩阵出现镜像或尺度错乱

这个坑的现象很好认：把H作用到参考帧上，输出的图像要么左右颠倒，要么被缩放到只剩一个角。原因集中在三处。最常见的是reshape的转置缺失，MATLAB 列优先，reshape(h, 3, 3)得到的是H的转置，表现出来像是沿对角线做了一个奇怪的错切，而不是简单的镜像，所以很多人第一反应不会怀疑这里。第二处是反归一化写反，H = T2 \ Hn * T1的顺序对应“先归一化输入、再归一化输出、最后反解”，写成T1 \ Hn * T2时在对称的归一化矩阵下有时看不出来，一旦两组点分布差异大就会崩。第三处是点坐标的排列顺序，(x, y)和(y, x)在A矩阵里对应完全不同的列，而代码里p(1,:)和p(2,:)哪个是 x 完全靠约定。

排查手段很直接：拿四对已知的正确匹配点做单元测试，H作用后重投影误差应该在1e-10量级；然后把det(H)打印出来，单应性矩阵的行列式在纯旋转平移下应接近 1，如果出现非常小或非常大的值，基本可以确定尺度出了问题。

5.3 块匹配在平坦区域疯狂跳变

现象是运动矢量场在天空、白墙、纯色背景这些区域出现随机方向、满幅值的大箭头，而同一帧的纹理区明明估得很好。原因就是 2.1 节说的代价面多峰：平坦区所有候选块的 SAD 都接近零，差值只在噪声量级，argmin取到哪个完全取决于数值抖动。这不是实现 bug，是问题本身欠定。

解决方案有三条，可以叠加使用。一是加运动矢量平滑惩罚，把能量从SAD(d)改成SAD(d) + λ·|d - d_pred|，d_pred用左邻和上邻的中值（中值比均值抗离群点），这也是实际编码器里做运动矢量预测的标准做法。二是给零矢量一个偏置，SAD(0) - bias，让平坦区默认落回静止，bias取 1 到 2 个灰度级就够。三是事后用中值滤波对矢量场做一次后处理，medfilt2分别作用在两个通道上，窗口 3×3，能吃掉大部分孤立跳变。

5.4 重构残差出现规律性网格

把预测帧和当前帧相减，残差图上出现整齐的方块边界，这是块效应的典型表现，原因不在估计而在补偿。块匹配是逐块独立决定位移的，相邻块可能选中完全不同的矢量，补偿时每个块整体搬运，边界处就出现不连续。块越大越明显，16×16的块在斜边上能差出好几个像素。

对症的修补是重叠块运动补偿（OBMC）：每个像素的预测值由覆盖它的多个块加权平均，权重在块中心最大、向边界递减，通常用升余弦窗。另一个来源是亚像素插值的核截断，双线性插值只有 2×2 支撑，在放大倍数高时高频损失明显，换成双三次或 6 抽头 Wiener 插值能缓解。还有一个容易漏的点是块大小与搜索窗的奇偶性，块边长为偶数时块中心落在半像素上，坐标取整会引入系统性的半像素偏差，所有块统一偏半个像素，残差上就表现为全图性的轻微网格。

5.5 迭代到几百轮代价不降反升

TV-L1 那类原始-对偶算法，能量曲线应该是单调下降的（对偶间隙单调收窄），一旦出现震荡上升，几乎一定是步长违反了对偶条件。tau * sigma * 8 >= 1时算法理论上就不保证收敛，实际表现是前几十轮正常，之后误差慢慢爬升，最后溢出。检查方法很简单，把tau和sigma各减半再跑，曲线立刻变平滑，说明问题就在这里。

第二个原因是金字塔层数不够。初始位移离真值太远，每层要处理的位移超过泰勒展开的适用范围，线性化误差积累，外层固定点迭代发散。判断依据是看最顶层输出的平均位移幅值，如果它已经超过该层图像尺寸的十分之一，说明顶层没把位移压下来，得加层。第三个原因是 warp 时用了最近邻插值，残差里出现阶梯状的不连续，sign(rho)在跳变处频繁翻转，u、v被反复推来推去。改成双线性是必须的，这不是精度优化，是收敛性的前提。

6. 取舍与验证：什么场景该用哪一种，怎么证明它真的更好

把十一种方法排在一起比较，最容易犯的错是拿 PSNR 当唯一尺子。PSNR 对运动边界的错位不敏感，一张整体偏移一像素的预测帧和一张边界精确但内部有噪声的预测帧，PSNR 可能前者更高，但人的观感和后续编码效率都是后者更好。这也是为什么这类对比实验里通常要同时报三个量：重构 PSNR、SSIM，以及对有真值的光流数据集的端点误差（EPE）。前两个衡量预测质量，后一个衡量运动估计本身准不准，三者不一致的时候，恰恰是最有信息量的地方。

选型的判断链条其实很短：先问场景是不是近似平面且相机有旋转，是就上单应性，不是就往下走；再问需不需要逐像素的稠密流，需要就看纹理和位移量，纹理够、位移小用 LK，位移大就套金字塔，边界重要就上 TV-L1；只关心全局平移且要求快，相位相关是性价比最高的选择；如果只是给编码器做块级预测，块匹配加中值平滑和矢量预测，工程上最稳。

验证环节有一个习惯值得保持：任何新方法进来，先在一对人工合成的图上跑。左边一张，右边把它按已知的平移、旋转、单应性搬过去，这样真值完全可控，估计出来的量能直接和参数对照。我见过太多次在真实序列上调半天参，最后发现是 warp 方向写反了，而这个问题在合成数据上五秒钟就能暴露。这份资源里的十一份 MATLAB 实现，最大的用处也在这里：它们把同一个问题的十一种视角摆在了一起，你可以用同一个合成用例去喂它们，看谁在哪种退化条件下先崩溃，这比读十篇综述更能建立直觉。从那以后我每次碰运动估计的代码，都强制先用合成图跑一遍单位测试，确认位移方向和量级都对得上，再往真实数据上放。希望帮到你。

本文还有配套的精品资源，点击获取