简介这份凸优化大作业代码包面向学习运动估计与计算机视觉课程的学生围绕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·(2S1)²·blk²)S16、blk16时单帧就要上亿次运算而这正是各种快速搜索算法出现的动机三步搜索、菱形搜索、六边形搜索本质都是在代价面上做贪心下降。这里有一个必须说清楚的取舍全搜索在整数格点上做的是穷举得到的是全局最小但那个“最小”只是离散格点上的最小快速搜索连这个保证都放弃了。代价面在纹理丰富区域通常是单峰的贪心还能收敛到附近一旦进入大面积平坦区或周期纹理栅栏、条纹、键盘代价面会出现多个几乎等高的极小值贪心搜索会随机落入其中一个这就是后面避坑章节要专门处理的问题。另一个常被忽略的点是亚像素。整数格点最优只说明真实位移落在某个半像素范围内直接用整数矢量做补偿残差里会残留系统性偏差。标准的修补方式是在最优整数点周围对代价面做一维或二维抛物线拟合% 在最优整数位移 (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; endnormalize_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(xu, yv, t1)其中(u, v)是这一点的位移。对右端在(x, y)处做一阶泰勒展开并丢掉高阶项就得到光流约束方程I_x · u I_y · v I_t 0I_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; endconf这个掩膜是 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的(xv, yu)处采样放到(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 endsign(rho)这一步就是 L1 数据项的全部内容残差为正就往负梯度方向推lambda*tau为负就往正方向推残差恰好为零则不动。这就是著名的软阈值收缩它天然产生“稀疏残差”——大部分像素的残差被压到零只有遮挡和真实运动不连续处保留非零值这正好符合我们对帧间预测残差的预期。TV 项那边的px ./ max(1, |p|)是投影到单位球保证对偶变量满足 TV 的对偶约束|p| ≤ 1。这里的tau和sigma不是随便填的Chambolle-Pock 型算法的收敛条件是tau * sigma * ||K||² 1K是梯度算子其谱范数上界是√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。前两个衡量预测质量后一个衡量运动估计本身准不准三者不一致的时候恰恰是最有信息量的地方。方法族自由度大位移输出密度主要失效场景全搜索块匹配2可扩展强受搜索窗限制块级平坦区、周期纹理快速搜索块匹配2中易陷局部极小块级纹理复杂且位移大相位相关2全局强单一位移多运动物体、非刚性Lucas-Kanade2局部弱需金字塔稀疏/半稠密弱纹理、大位移Horn-Schunck2逐点弱需金字塔稠密运动边界模糊TV-L1 变分2逐点中需金字塔稠密参数敏感、阶梯效应仿射配准6强全局非平面、多物体单应性 RANSAC8强全局点少、近共线、多平面选型的判断链条其实很短先问场景是不是近似平面且相机有旋转是就上单应性不是就往下走再问需不需要逐像素的稠密流需要就看纹理和位移量纹理够、位移小用 LK位移大就套金字塔边界重要就上 TV-L1只关心全局平移且要求快相位相关是性价比最高的选择如果只是给编码器做块级预测块匹配加中值平滑和矢量预测工程上最稳。验证环节有一个习惯值得保持任何新方法进来先在一对人工合成的图上跑。左边一张右边把它按已知的平移、旋转、单应性搬过去这样真值完全可控估计出来的量能直接和参数对照。我见过太多次在真实序列上调半天参最后发现是 warp 方向写反了而这个问题在合成数据上五秒钟就能暴露。这份资源里的十一份 MATLAB 实现最大的用处也在这里它们把同一个问题的十一种视角摆在了一起你可以用同一个合成用例去喂它们看谁在哪种退化条件下先崩溃这比读十篇综述更能建立直觉。从那以后我每次碰运动估计的代码都强制先用合成图跑一遍单位测试确认位移方向和量级都对得上再往真实数据上放。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/27 7:32:21 · 阅读 1.2k 帧间预测十一种方法解析：运动模型、正则项与凸优化避坑
简介这份凸优化大作业代码包面向学习运动估计与计算机视觉课程的学生围绕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·(2S1)²·blk²)S16、blk16时单帧就要上亿次运算而这正是各种快速搜索算法出现的动机三步搜索、菱形搜索、六边形搜索本质都是在代价面上做贪心下降。这里有一个必须说清楚的取舍全搜索在整数格点上做的是穷举得到的是全局最小但那个“最小”只是离散格点上的最小快速搜索连这个保证都放弃了。代价面在纹理丰富区域通常是单峰的贪心还能收敛到附近一旦进入大面积平坦区或周期纹理栅栏、条纹、键盘代价面会出现多个几乎等高的极小值贪心搜索会随机落入其中一个这就是后面避坑章节要专门处理的问题。另一个常被忽略的点是亚像素。整数格点最优只说明真实位移落在某个半像素范围内直接用整数矢量做补偿残差里会残留系统性偏差。标准的修补方式是在最优整数点周围对代价面做一维或二维抛物线拟合% 在最优整数位移 (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; endnormalize_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(xu, yv, t1)其中(u, v)是这一点的位移。对右端在(x, y)处做一阶泰勒展开并丢掉高阶项就得到光流约束方程I_x · u I_y · v I_t 0I_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; endconf这个掩膜是 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的(xv, yu)处采样放到(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 endsign(rho)这一步就是 L1 数据项的全部内容残差为正就往负梯度方向推lambda*tau为负就往正方向推残差恰好为零则不动。这就是著名的软阈值收缩它天然产生“稀疏残差”——大部分像素的残差被压到零只有遮挡和真实运动不连续处保留非零值这正好符合我们对帧间预测残差的预期。TV 项那边的px ./ max(1, |p|)是投影到单位球保证对偶变量满足 TV 的对偶约束|p| ≤ 1。这里的tau和sigma不是随便填的Chambolle-Pock 型算法的收敛条件是tau * sigma * ||K||² 1K是梯度算子其谱范数上界是√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。前两个衡量预测质量后一个衡量运动估计本身准不准三者不一致的时候恰恰是最有信息量的地方。方法族自由度大位移输出密度主要失效场景全搜索块匹配2可扩展强受搜索窗限制块级平坦区、周期纹理快速搜索块匹配2中易陷局部极小块级纹理复杂且位移大相位相关2全局强单一位移多运动物体、非刚性Lucas-Kanade2局部弱需金字塔稀疏/半稠密弱纹理、大位移Horn-Schunck2逐点弱需金字塔稠密运动边界模糊TV-L1 变分2逐点中需金字塔稠密参数敏感、阶梯效应仿射配准6强全局非平面、多物体单应性 RANSAC8强全局点少、近共线、多平面选型的判断链条其实很短先问场景是不是近似平面且相机有旋转是就上单应性不是就往下走再问需不需要逐像素的稠密流需要就看纹理和位移量纹理够、位移小用 LK位移大就套金字塔边界重要就上 TV-L1只关心全局平移且要求快相位相关是性价比最高的选择如果只是给编码器做块级预测块匹配加中值平滑和矢量预测工程上最稳。验证环节有一个习惯值得保持任何新方法进来先在一对人工合成的图上跑。左边一张右边把它按已知的平移、旋转、单应性搬过去这样真值完全可控估计出来的量能直接和参数对照。我见过太多次在真实序列上调半天参最后发现是 warp 方向写反了而这个问题在合成数据上五秒钟就能暴露。这份资源里的十一份 MATLAB 实现最大的用处也在这里它们把同一个问题的十一种视角摆在了一起你可以用同一个合成用例去喂它们看谁在哪种退化条件下先崩溃这比读十篇综述更能建立直觉。从那以后我每次碰运动估计的代码都强制先用合成图跑一遍单位测试确认位移方向和量级都对得上再往真实数据上放。希望帮到你。本文还有配套的精品资源点击获取