Lucas-Kanade光流法详解

简介：个人学习分享，如有错误，欢迎批评指正。

光流（Optical Flow）描述的是图像序列中各像素点随时间的运动情况，是计算机视觉中的基本问题之一。光流问题涉及尝试找出一幅图像中的许多点在第二幅图像中移动的位置–通常是以视频序列完成的，因此可以假定第一幅图像中的大部分点框架都可以在第二幅图像中找到。光流可以用于场景中物体的运动估计，设置用于相机相对于整个场景的自运动估计。光流算法的理想输出是两帧图像中每个像素的速度的估计关联，或者等效的，一幅图像中每个像素的位移矢量，指示该像素在另一幅图像中的相对位置，如果图像中的每个像素都采用这种方法通常称为“稠密光流”。还有一种算法称为“稀疏光流”，仅仅只跟踪图像中某些点的子集，该算法通常是快速且可靠的，因为其将注意力只放在容易跟踪的特定点上，稀疏跟踪的计算成本远远低于稠密跟踪，这也导致了后者的研究基本被限制在了学术圈。

一、基本假设

Lucas-Kanade方法基于以下三个主要假设，这些假设简化了光流估计问题，使其在实际应用中可行：

• 亮度恒定假设（Brightness Constancy Assumption）：

  • 定义：同一物体点在连续两帧图像中的亮度值保持不变。

  • 数学表达：对于任意时间$t$，点$(x,y)$在$t$和$t+\Delta t$时刻的亮度满足
    $$I(x,y,t)=I(x+\Delta x,y+\Delta y,t+\Delta t)$$

  • 意义：将物体的运动与亮度变化联系起来，为光流估计提供约束条件。

• 小位移假设（Small Motion Assumption）：

  • 定义：物体的运动在时间间隔$\Delta t$内是微小的，位移$(\Delta x,\Delta y)$较小。

  • 数学表达：$\Delta x$和$\Delta y$相对于图像尺度而言足够小。

  • 意义：允许对图像亮度进行泰勒展开，并忽略高阶项，简化光流方程的线性化过程。

• 局部光流一致性假设（Local Flow Constancy Assumption）：

  • 定义：在一个小的邻域窗口内，所有像素点的光流向量相同。
  • 数学表达：对于窗口内的所有点$(x_i,y_i)$，光流向量$(u,v)$相同。
  • 意义：将每个像素的光流问题转化为局部一致的运动估计，简化求解过程。

这些假设为光流估计提供了必要的约束，使得原本不可解的方程组变为可解的形式。然而，这些假设在实际应用中可能不完全成立，导致算法在某些场景下性能下降。

二、数学原理

1.亮度恒定方程的推导

基于亮度恒定假设，Lucas-Kanade方法首先建立亮度恒定方程，连接图像的空间变化与时间变化。

1.1 亮度恒定假设

设$I(x,y,t)$表示在时间$t$的图像，其中$(x,y)$是像素坐标。对于一个物体点$(x,y)$，其在时间$t$和$t+\Delta t$的亮度满足：
$$I(x,y,t)=I(x+\Delta x,y+\Delta y,t+\Delta t)$$

其中，$\Delta x$和$\Delta y$分别是物体点在$x$和$y$方向上的位移。

1.2 泰勒展开

由于假设$\Delta t$较小，位移$\Delta x$和$\Delta y$也较小，可以对$I(x+\Delta x,y+\Delta y,t+\Delta t)$在$(x,y,t)$处进行泰勒展开，并忽略高阶项（即$O(\Delta x^2,\Delta y^2,\Delta t^2)$）：
$$I(x+\Delta x,y+\Delta y,t+\Delta t)\approx I(x,y,t)+I_x\Delta x+I_y\Delta y+I_t\Delta t$$

其中，$I_x=\frac{\partial I}{\partial x}$，$I_y=\frac{\partial I}{\partial y}$，$I_t=\frac{\partial I}{\partial t}$分别表示图像在$x$、$y$方向的空间梯度和时间梯度。

1.3 代入亮度恒定方程

根据亮度恒定假设，将泰勒展开式代入，得到：

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

两边减去$I(x,y,t)$，得到：
$$I_x\Delta x+I_y\Delta y+I_t\Delta t=0$$

1.4 引入光流向量

定义光流向量$(u,v)$为位移$(\Delta x,\Delta y)$与时间间隔$\Delta t$的比值：

$$u=\frac{\Delta x}{\Delta t},v=\frac{\Delta y}{\Delta t}$$

将其代入上式，得到：
$$I_{x}u+I_{y}v+I_t=0$$

这就是光流约束方程，它将图像的空间梯度、时间梯度与光流向量联系起来。其中$(u,v)$代表两个方向(x方向和y方向)的移动速度,把$(u,v)$计算出来，咱们的光流也就算出来了。

拿到当前帧，假设我们要计算点p的光流。其中$I_x$,$I_y$都可以通过当前帧计算出来，而 $I_t$ 可以通过两帧的差分计算出来。所以，对于公式 $I_{x}u+I_{y}v+I_t=0$ 而言，未知数只有 $u$ 和 $v$。

LK算法就是用来求解公式 $I_{x}u+I_{y}v+I_t=0$ 的。LK有一个window的概念，即我先划定一块区域比如(5x5)的像素区域，我们可以认为这块区域每个点的移动速度$u$ 和 $v$是一致的。
首先，$I_x$,$I_y$是怎么得到的呢？对于光流法，咱们有个理想的假定就是：运动物体只会做平移。所以，亮度梯度咱们只需考虑当前帧的梯度即可。对于 $I_t$ 我们在两帧做一个差分就可以得到。

咱们看当前帧，也就是右边那个。$I_x(3,3)=0$ 是因为在x轴的数值左右都是3，没有梯度变化。$I_y(3,3)=1$ 是因为在y轴数值变化幅度为1。上图显示的情况，咱们只能算出y轴速度v，没有办法算出水平（x轴）速度u。这是因为该移动目标本身在x轴方向上就没有亮度变化。这也是一种典型的问题，叫孔径问题(Aperture Problem)。
孔径问题是讲，如果我们通过一个小孔来看全局，很多情况下的移动信息我们是看不出来的，这个用一张图就可以很好理解：

假设墙上破了一个洞，咱们通过这个洞来看一个图形的移动情况。假设，我们看到的是上图这种情况，绿色部分就是我们的视野，我们无法判断这个图形是否是沿着切线方向移动或者是静止。
那么基于小区域的LK光流法也可能遇到Aperture Problem，所以我们在追光流的时候，选点通常会选目标的角点(corner)。角点的情况如下图：

如果角点在视野内的话，咱们就可以判断这个图形的运动方向。
接着公式 $I_{x}u+I_{y}v+I_t=0$ ，咱们如果通过window方式求解u、v，那还是很好办的。假设咱们取的是5x5的window，那么window内的每个点，我们都认为有一样的移动方向，咱们可以构建出25个等式。求解二元一次方程，通常两个等式就可以求解。但那是理想情况，实际情况是没有一组( u , v ) 能同时满足这25个等式，咱们要做的是最小化这个差异。

写成矩阵形式：

只需找到一组( u , v )，即上图中的x，满足

公式(5)
这里就可以用最小二乘法来进行优化了，分别求偏导得出导数为0的点的值就是最优值。可推导出，

这便是Lucas-Kanade光流算法的公式了。

2.光流约束方程的建立

对于上述问题进行重新表述：在单个像素点，光流约束方程 $I_{x}u+I_{y}v+I_t=0$ 包含两个未知数$u$和$v$，无法直接求解。Lucas-Kanade方法通过在一个小窗口内（如$3\times 3$像素）对多个像素点的光流约束方程进行组合，形成一个过定的线性方程组，从而估计出统一的$u$和$v$。

2.1 局部窗口内的光流方程

假设在一个$N\times N$的窗口内（共$n=N^2$个像素点），每个像素点$(x_i,y_i)$都有一个光流约束方程：

$$I_{x_i}u+I_{y_i}v+I_{t_i}=0,i=1,2,\dots ,n$$

将这些方程表示为矩阵形式：

$$\left[\begin{array}{cc}{I}_{x_1}& I_{y_1}\\ I_{x_2}& I_{y_2}\\ ⋮& ⋮\\ I_{x_n}& I_{y_n}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}-I_{t_1}\\ -I_{t_2}\\ ⋮\\ -I_{t_n}\end{array}\right]$$

记作：
$$A\cdot \mathbf{d}=\mathbf{b}$$

其中：
$$A=\left[\begin{array}{cc}{I}_{x_1}& I_{y_1}\\ I_{x_2}& I_{y_2}\\ ⋮& ⋮\\ I_{x_n}& I_{y_n}\end{array}\right],\mathbf{d}=\left[\begin{array}{c}u\\ v\end{array}\right],\mathbf{b}=\left[\begin{array}{c}-I_{t_1}\\ -I_{t_2}\\ ⋮\\ -I_{t_n}\end{array}\right]$$

2.2 方程组的特性

• 过定性：由于$n>2$，方程组为过定的，存在多个方程与两个未知数。
• 最小二乘解：通过最小化误差平方和，求解最优的$u$和$v$。

3.最小二乘法求解

为了求解光流向量$(u,v)$，Lucas-Kanade方法采用最小二乘法，具体步骤如下：

3.1 最小二乘法基本原理

目标是找到$\mathbf{d}=\left[\begin{array}{c}uv\end{array}\right]$，使得误差平方和最小：

$$E(\mathbf{d})=\sum _{i=1}^n(I_{x_i}u+I_{y_i}v+I_{t_i}{)}^2=\Vert A\mathbf{d}-\mathbf{b}{\Vert }^2$$

3.2 正规方程（Normal Equations）

通过对$E(\mathbf{d})$关于$\mathbf{d}$求导并设为零，可以得到正规方程：

$$A^{T}A\mathbf{d}=A^T\mathbf{b}$$

3.3 解的存在性与唯一性

• 解存在的条件：矩阵$A^{T}A$是对称且半正定的。
• 解唯一性的条件：$A^{T}A$是正定的，即窗口内的梯度信息足够丰富（线性无关），矩阵$A$的秩为2。

如果$A^{T}A$是非奇异的（可逆），则正规方程有唯一解：

$$\mathbf{d}=(A^{T}A{)}^{-1}{A}^T\mathbf{b}$$

3.4 解析解与数值求解方法

3.4.1 直接求逆

对于小窗口（如$3\times 3$），可以直接计算$(A^{T}A{)}^{-1}{A}^T\mathbf{b}$，得到光流向量。

3.4.2 矩阵分解方法

• QR分解：

  • 分解$A$为$QR$，其中$Q$是正交矩阵，$R$是上三角矩阵。
  • 解方程$R\mathbf{d}=Q^T\mathbf{b}$。

• 奇异值分解（SVD）：

  • 将$A$分解为$U\Sigma V^T$，其中$U$和$V$是正交矩阵，$\Sigma$是对角矩阵。

  • 最小二乘解为$\mathbf{d}=V{\Sigma }^{-1}{U}^T\mathbf{b}$。

  • 优点：提高数值稳定性，尤其在$A^{T}A$接近奇异时。

• Cholesky分解：

  • 适用于正定矩阵$A^{T}A$，通过分解为下三角矩阵和其转置的乘积，快速求解。

3.4.3 数值稳定性

在实际计算中，$A^{T}A$可能接近奇异，导致数值不稳定。为提高稳定性，可以采取以下措施：

• 增加窗口内像素点数量：扩大窗口尺寸，提供更多约束条件。
• 正则化：引入Tikhonov正则化（岭回归），在正规方程中加入正则项：
  $$\mathbf{d}=(A^{T}A+\lambda I{)}^{-1}{A}^T\mathbf{b}$$

其中，$\lambda$是正则化参数，$I$是单位矩阵。正则化可以改善矩阵$A^{T}A$的条件数，增强解的稳定性。

鲁棒估计方法：如RANSAC，通过迭代剔除异常值，提高求解的鲁棒性。

3.5 最终光流向量

通过最小二乘法求解，得到光流向量：

$$\mathbf{d}=\left[\begin{array}{c}u\\ v\end{array}\right]=(A^{T}A{)}^{-1}{A}^T\mathbf{b}$$

其中，$u$和$v$分别表示光流在$x$和$y$方向上的分量。

4.矩阵性质与数值稳定性

为了确保最小二乘解的有效性和稳定性，需要深入理解矩阵$A$及其转置乘积$A^{T}A$的性质。

4.1 矩阵$A$的结构

矩阵$A$的每一行对应窗口内一个像素点的空间梯度：

$$A=\left[\begin{array}{cc}{I}_{x_1}& I_{y_1}\\ I_{x_2}& I_{y_2}\\ ⋮& ⋮\\ I_{x_n}& I_{y_n}\end{array}\right]$$

4.2 矩阵$A^{T}A$的性质

• 对称性：$A^{T}A$是对称矩阵。
• 正定性：
  • 如果窗口内的梯度向量$\left[\begin{array}{c}{I}_{xi}{I}_{yi}\end{array}\right]$线性无关，则$A^{T}A$是正定的。
  • 否则，$A^{T}A$是半正定的。

4.3 条件数与数值稳定性

• 条件数：矩阵$A^{T}A$的条件数衡量了其数值稳定性。条件数越大，矩阵越接近奇异，数值求解越不稳定。
• 梯度分布：窗口内梯度向量的分布影响$A^{T}A$的条件数。例如：
  • 高纹理区域：梯度分布丰富，$A^{T}A$的条件数较小，数值稳定。
  • 低纹理区域：梯度分布稀疏或线性相关，$A^{T}A$的条件数较大，可能导致数值不稳定。

4.4 解决数值不稳定的方法

• 窗口大小调整：在低纹理区域使用更大的窗口，增加方程组的约束。
• 正则化：如前述，加入正则项改善条件数。
• 奇异值分解：使用SVD求解最小二乘问题，避免直接求逆带来的数值不稳定。
• 多尺度方法：通过构建图像金字塔，从粗到细估计光流，缓解局部条件数较差的问题。

三、算法步骤

Lucas-Kanade方法的具体实现包括以下主要步骤，每个步骤都包含若干子步骤和技术细节：

1. 图像预处理

图像预处理是Lucas-Kanade方法中的首要步骤，其主要目的是提高光流估计的准确性和鲁棒性。预处理通常包括图像平滑和梯度计算。

1.1. 图像平滑

目的：减小图像噪声对光流估计的影响，提升梯度计算的稳定性。

方法：

• 高斯滤波：最常用的平滑方法，通过卷积一个高斯核来平滑图像，减少高频噪声。
  $$G(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

其中，$\sigma$是高斯核的标准差，控制平滑的程度。

• 其他平滑方法：如中值滤波、双边滤波等，根据具体应用需求选择。

实现细节：

• 核大小选择：核大小通常选择为奇数（如3×3、5×5、7×7），较大的核提供更强的平滑效果，但可能模糊细节。
• 边缘处理：常用的方法包括镜像扩展、零填充等，以处理图像边缘的卷积操作。

1.2. 梯度计算

目的：计算图像在$x$、$y$方向的空间梯度 $I_x$ 和 $I_y$，以及时间梯度$I_t$，为光流约束方程提供必要的信息。

方法：

• 空间梯度：

  • Sobel算子：
    $$I_x=\left[\begin{array}{ccc}-1& 0& +1\\ -2& 0& +2\\ -1& 0& +1\end{array}\right]\ast I{I}_y=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ +1& +2& +1\end{array}\right]\ast I$$

  • Scharr算子：类似于Sobel算子，但具有更好的旋转对称性和梯度响应。

  • Prewitt算子：另一种常见的梯度算子，与Sobel类似。

• 时间梯度：
  • 计算连续两帧图像$I_t(x,y)=I(x,y,t+\Delta t)-I(x,y,t)$。
  • 可以使用更精确的方法，如中心差分法，计算时间梯度。

实现细节：

• 边缘处理：与平滑步骤类似，需要处理图像边缘的梯度计算。

• 梯度精度：选择合适的梯度算子和参数，以平衡梯度的精度和计算复杂度。

2. 窗口选择

窗口选择是Lucas-Kanade方法中的关键步骤，决定了局部光流估计的区域范围和光流一致性的假设。

2.1. 窗口类型

• 固定窗口：常用的固定大小窗口（如3×3、5×5、7×7），适用于光流估计假设有效的区域。
• 自适应窗口：根据图像内容动态调整窗口大小，如在纹理丰富区域使用较小窗口，在纹理稀疏区域使用较大窗口。

2.2. 窗口形状

• 矩形窗口：最常用的窗口形状，覆盖一个矩形区域。
• 其他形状：如圆形窗口，可以更好地适应某些应用需求。

2.3. 窗口中心的选择

• 稠密光流：为图像中的每个像素点选择窗口，生成稠密的光流场。
• 稀疏光流：仅为图像中的特征点（如角点、边缘点）选择窗口，减少计算量。

2.4. 窗口大小的选择

影响因素：

• 运动范围：较大的窗口能够捕捉更大范围的运动，但可能降低局部精度。
• 纹理信息：较大的窗口在纹理稀疏区域更有效，较小窗口在纹理丰富区域更精准。

选择策略：

• 经验法则：根据应用场景和图像特性选择合适的窗口大小。
• 多尺度方法：结合金字塔Lucas-Kanade方法，通过多尺度处理适应不同窗口大小。

3. 构建局部方程组

在选定的窗口内，构建光流约束方程组，以求解统一的光流向量。

3.1. 收集窗口内的梯度和时间梯度

对于窗口内的每个像素点$(x_i,y_i)$，收集其空间梯度 $I_{xi}$、$I_{yi}$和时间梯度 $I_{ti}$。

3.2. 构建矩阵$ A $和向量$ b $

根据光流约束方程 $I_{x}u+I_{y}v=-I_t$，为窗口内的每个像素点构建一个方程，形成过定的线性方程组：

$$A=\left[\begin{array}{cc}{I}_{x1}& I_{y1}\\ I_{x2}& I_{y2}\\ ⋮& ⋮\\ I_{xn}& I_{yn}\end{array}\right],b=\left[\begin{array}{c}-I_{t1}\\ -I_{t2}\\ ⋮\\ -I_{tn}\end{array}\right]$$

其中，$ n $是窗口内像素点的数量。

3.3. 数据质量检查

在构建方程组之前，需检查数据质量：

• 梯度的非零性：确保窗口内存在足够的梯度信息，避免平坦区域导致矩阵$A$退化。
• 异常值检测：剔除异常的梯度值或时间梯度，提升方程组的稳定性。

4. 最小二乘解法

通过最小二乘法求解过定的线性方程组，获得光流向量$(u,v)$。

4.1. 最小二乘法基本原理

目标是找到 $(u,v)$，使得 $|A\cdot \left[\begin{array}{c}uv\end{array}\right]-b{|}^2$最小，即最小化误差平方和。

4.2. 解析解

当矩阵$A^{T}A$可逆时，解析解为：

4.3 数值求解方法

直接求逆：在小窗口情况下（如3×3），可以直接计算$(A^{T}A{)}^{-1}{A}^{T}b$。

矩阵分解：

• QR分解：分解$A$为$QR$，然后解$R\left[\begin{array}{c}uv\end{array}\right]=Q^{T}b$。
• 奇异值分解（SVD）：将$A$分解为$U\Sigma V^T$，避免直接求逆，提高数值稳定性。
• Cholesky分解：用于正定矩阵$A^{T}A$的高效分解。

梯度下降法：在大规模问题或需要迭代优化时使用，但在小窗口情况下通常不必要。

4.4 正则化（可选）

为处理矩阵$A^{T}A$接近奇异的情况，可以引入正则化项，如Tikhonov正则化：

$$\left[\begin{array}{c}u\\ v\end{array}\right]=(A^{T}A+\lambda I{)}^{-1}{A}^{T}b$$

其中，$\lambda$是正则化参数，$I$是单位矩阵。

4.5. 解的有效性检查

• 可逆性检查：确保$A^{T}A$的条件数合理，避免数值不稳定。
• 异常值检测：在解出$(u,v)$后，检查是否存在异常值，如过大的光流向量，可能需要进一步处理。

5. 光流向量估计

将求解得到的光流向量$(u,v)$赋值给窗口中心的像素点，完成局部光流估计。

5.1 光流场的构建

通过遍历整个图像，对每个需要计算光流的像素点执行上述步骤，构建整个图像的光流场。

5.2 数据结构选择

• 稀疏光流：存储为特征点及其光流向量的列表，适用于稀疏光流估计。
• 稠密光流：存储为与图像同尺寸的二维向量场，适用于稠密光流估计。

5.3 插值与填充（可选）

在某些应用中，可能需要对光流场进行插值或填充，以处理缺失的光流向量或增强光流场的连续性。

5.4 后处理（可选）

• 光流平滑：对光流场进行平滑处理，减少噪声影响。
• 边缘保留：在光流估计中保留运动边界，避免模糊运动区域。
• 一致性检查：通过前向和后向光流的一致性，检测和剔除错误的光流向量。

6. 整体流程概述

以下是Lucas-Kanade方法的整体算法流程：

1.输入：连续两帧图像$I(t)$和$I(t+\Delta t)$。

2.预处理：

• 对$I(t)$和$I(t+\Delta t)$进行高斯平滑。
• 计算$I_x$、$I_y$、$I_t$。

3.光流估计：

对每个需要计算光流的像素点：

• 选择包含该点的窗口。
• 构建局部光流约束方程组。
• 通过最小二乘法求解$(u,v)$。
• 赋值给窗口中心点。

4.输出：光流场，包含每个像素点的光流向量$(u,v)$。

7.实例分析

1. 示例场景
假设在一个$3\times 3$的窗口内，有以下梯度信息：

2. 构建矩阵$A$和向量$b$
根据上表，矩阵$A$和向量$b$如下：

$$A=\left[\begin{array}{cc}1& 0\\ 2& 0\\ 1& 0\\ 0& 1\\ 0& 1\\ 0& 1\\ 1& 0\\ 2& 0\\ 1& 0\end{array}\right],b=\left[\begin{array}{c}-(-1)\\ -(-2)\\ -(-1)\\ -0\\ -0\\ -0\\ -(-1)\\ -(-2)\\ -(-1)\end{array}\right]=\left[\begin{array}{c}1\\ 2\\ 1\\ 0\\ 0\\ 0\\ 1\\ 2\\ 1\end{array}\right]$$

3. 计算$A^{T}A$和$A^{T}b$

$$A^{T}A=\left[\begin{array}{ccccccccc}1& 2& 1& 0& 0& 0& 1& 2& 1\\ 0& 0& 0& 1& 1& 1& 0& 0& 0\end{array}\right]\left[\begin{array}{cc}1& 0\\ 2& 0\\ 1& 0\\ 0& 1\\ 0& 1\\ 0& 1\\ 1& 0\\ 2& 0\\ 1& 0\end{array}\right]=\left[\begin{array}{cc}{1}^2+2^2+1^2+0^2+0^2+0^2+1^2+2^2+1^2& 1\cdot 0+2\cdot 0+1\cdot 0+0\cdot 1+0\cdot 1+0\cdot 1+1\cdot 0+2\cdot 0+1\cdot 0\\ 0\cdot 1+0\cdot 2+0\cdot 1+1\cdot 0+1\cdot 0+1\cdot 0+0\cdot 1+0\cdot 2+0\cdot 1& 0^2+0^2+0^2+1^2+1^2+1^2+0^2+0^2+0^2\end{array}\right]=\left[\begin{array}{cc}12& 0\\ 0& 3\end{array}\right]$$

$$A^{T}b=\left[\begin{array}{ccccccccc}1& 2& 1& 0& 0& 0& 1& 2& 1\\ 0& 0& 0& 1& 1& 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}1\\ 2\\ 1\\ 0\\ 0\\ 0\\ 1\\ 2\\ 1\end{array}\right]=\left[\begin{array}{c}1\cdot 1+2\cdot 2+1\cdot 1+0\cdot 0+0\cdot 0+0\cdot 0+1\cdot 1+2\cdot 2+1\cdot 1\\ 0\cdot 1+0\cdot 2+0\cdot 1+1\cdot 0+1\cdot 0+1\cdot 0+0\cdot 1+0\cdot 2+0\cdot 1\end{array}\right]=\left[\begin{array}{c}1+4+1+0+0+0+1+4+1\\ 0+0+0+0+0+0+0+0+0\end{array}\right]=\left[\begin{array}{c}12\\ 0\end{array}\right]$$

4. 求解正规方程
正规方程为：
$$A^{T}Ad=A^{T}b⟹\left[\begin{array}{cc}12& 0\\ 0& 3\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}12\\ 0\end{array}\right]$$

解得：

$$\begin{gathered} 12u=12⟹u=1 \\ 3v=0⟹v=0 \end{gathered}$$

因此，光流向量为$(u,v)=(1,0)$，表示图像在$x$方向上以单位速度移动，而在$y$方向上无运动。

5. 分析

• 矩阵$A^{T}A$的可逆性：在本例中，$A^{T}A$是对角矩阵，且对角元素均大于零，故可逆。
• 条件数：$A^{T}A$的条件数为$\frac{{\lambda }_{\text{max}}}{{\lambda }_{\text{min}}}=\frac{12}{3}=4$，较小，数值稳定。
• 光流方向：光流仅在$x$方向上有运动，符合光流向量的计算结果。

四、代码实例

简单版

import cv2
import numpy as np

def lucas_kanade_optical_flow(I1, I2, window_size=5, max_iterations=10):
    # 1. Preprocessing: Convert to grayscale if necessary
    if len(I1.shape) == 3:
        I1 = cv2.cvtColor(I1, cv2.COLOR_BGR2GRAY)
    if len(I2.shape) == 3:
        I2 = cv2.cvtColor(I2, cv2.COLOR_BGR2GRAY)
    
    # 2. Smooth images
    I1_blur = cv2.GaussianBlur(I1, (5, 5), sigmaX=1.5)
    I2_blur = cv2.GaussianBlur(I2, (5, 5), sigmaX=1.5)
    
    # 3. Compute gradients
    Ix = cv2.Sobel(I1_blur, cv2.CV_64F, 1, 0, ksize=3)
    Iy = cv2.Sobel(I1_blur, cv2.CV_64F, 0, 1, ksize=3)
    It = I2_blur - I1_blur
    
    # 4. Initialize flow field
    flow_u = np.zeros(I1.shape)
    flow_v = np.zeros(I1.shape)
    
    # 5. Iterate over each pixel
    half_w = window_size // 2
    for y in range(half_w, I1.shape[0] - half_w):
        for x in range(half_w, I1.shape[1] - half_w):
            # Define window
            Ix_window = Ix[y - half_w : y + half_w + 1, x - half_w : x + half_w + 1].flatten()
            Iy_window = Iy[y - half_w : y + half_w + 1, x - half_w : x + half_w + 1].flatten()
            It_window = It[y - half_w : y + half_w + 1, x - half_w : x + half_w + 1].flatten()
            
            # Construct A and b
            A = np.vstack((Ix_window, Iy_window)).T  # Shape: (window_size^2, 2)
            b = -It_window  # Shape: (window_size^2,)
            
            # Check if A^T A is invertible
            ATA = A.T @ A
            if np.linalg.cond(ATA) < 1 / np.finfo(float).eps:
                # Solve for (u, v) using least squares
                nu, nv = np.linalg.lstsq(A, b, rcond=None)[0]
                flow_u[y, x] = nu
                flow_v[y, x] = nv
            else:
                flow_u[y, x] = 0
                flow_v[y, x] = 0
    
    # Combine flow fields
    flow = np.stack((flow_u, flow_v), axis=2)
    return flow

# 示例使用
if __name__ == "__main__":
    # 读取连续两帧图像
    I1 = cv2.imread('frame1.png')
    I2 = cv2.imread('frame2.png')
    
    # 计算光流
    flow = lucas_kanade_optical_flow(I1, I2, window_size=5)
    
    # 可视化光流
    hsv = np.zeros_like(I1)
    hsv = cv2.cvtColor(hsv, cv2.COLOR_BGR2HSV)
    mag, ang = cv2.cartToPolar(flow[...,0], flow[...,1])
    hsv[...,0] = ang * 180 / np.pi / 2
    hsv[...,1] = 255
    hsv[...,2] = cv2.normalize(mag, None, 0, 255, cv2.NORM_MINMAX)
    rgb_flow = cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)
    
    cv2.imshow('Optical Flow', rgb_flow)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

复杂版：图像对齐融合中的光流法

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


def read_dng(file_path, scale=1.0):
    """
    读取DNG图像并转换为RGB和灰度图像。

    Args:
        file_path (str): DNG图像文件路径。
        scale (float): 图像缩放比例，默认为1.0（不缩放）。

    Returns:
        tuple: (RGB图像, 灰度图像)
    """
    with rawpy.imread(file_path) as raw:
        # 使用raw.postprocess将RAW数据转换为RGB图像，设置输出位深为16位
        rgb = raw.postprocess(use_camera_wb=True, no_auto_bright=True, output_bps=16)

    if scale != 1.0:
        rgb = cv2.resize(rgb, None, fx=scale, fy=scale, interpolation=cv2.INTER_LINEAR)

    # 将RGB图像转换为灰度图
    gray = cv2.cvtColor(rgb, cv2.COLOR_RGB2GRAY)

    # 检查灰度图像的类型
    if gray.dtype == np.uint16:
        # 16-bit到8-bit的转换，归一化
        gray = cv2.convertScaleAbs(gray, alpha=(255.0 / 65535.0))
    elif gray.dtype != np.uint8:
        # 如果不是8位，进行归一化并转换为8位
        gray = cv2.normalize(gray, None, 0, 255, cv2.NORM_MINMAX)
        gray = gray.astype(np.uint8)

    # 高斯平滑，减少噪声
    gray = cv2.GaussianBlur(gray, (5, 5), 0)

    return rgb, gray


def detect_features(gray, maxCorners=500, qualityLevel=0.3, minDistance=7, blockSize=7):
    """
    使用Shi-Tomasi方法检测特征点。

    Args:
        gray (ndarray): 灰度图像。
        maxCorners (int): 最大特征点数。
        qualityLevel (float): 质量水平阈值。
        minDistance (float): 特征点之间的最小距离。
        blockSize (int): 阻塞大小。

    Returns:
        ndarray: 检测到的特征点坐标。
    """
    feature_params = dict(maxCorners=maxCorners,
                          qualityLevel=qualityLevel,
                          minDistance=minDistance,
                          blockSize=blockSize)
    p = cv2.goodFeaturesToTrack(gray, mask=None, **feature_params)
    return p


def track_features(old_gray, new_gray, p0, lk_params):
    """
    使用Lucas-Kanade方法跟踪特征点。

    Args:
        old_gray (ndarray): 上一帧的灰度图像。
        new_gray (ndarray): 当前帧的灰度图像。
        p0 (ndarray): 上一帧的特征点坐标。
        lk_params (dict): Lucas-Kanade方法参数。

    Returns:
        tuple: (新特征点, 上一帧的特征点, 状态)
    """
    p1, st, err = cv2.calcOpticalFlowPyrLK(old_gray, new_gray, p0, None, **lk_params)
    if p1 is None:
        return None, None, None
    good_new = p1[st == 1]
    good_old = p0[st == 1]
    return good_new, good_old, st


def estimate_global_transform(good_new, good_old):
    """
    使用RANSAC方法估计全局仿射变换矩阵。

    Args:
        good_new (ndarray): 当前帧的特征点。
        good_old (ndarray): 参考帧的特征点。

    Returns:
        tuple: (变换矩阵, inliers mask)
    """
    if len(good_new) < 3:
        return None, None
    M, mask = cv2.estimateAffinePartial2D(good_old, good_new, method=cv2.RANSAC, ransacReprojThreshold=5.0)
    return M, mask


def estimate_local_transforms(good_new, good_old, grid_size=(8, 8), block_size=50):
    """
    分块估计局部仿射变换矩阵。

    Args:
        good_new (ndarray): 当前帧的特征点。
        good_old (ndarray): 参考帧的特征点。
        grid_size (tuple): 网格大小（x, y）。
        block_size (int): 每个块的大小。

    Returns:
        dict: 每个块的变换矩阵。
    """
    transforms = {}
    for x in range(grid_size[0]):
        for y in range(grid_size[1]):
            # 定义块的区域
            mask = (good_new[:, 0] > x * block_size) & (good_new[:, 0] < (x + 1) * block_size) & \
                   (good_new[:, 1] > y * block_size) & (good_new[:, 1] < (y + 1) * block_size)
            block_new = good_new[mask]
            block_old = good_old[mask]
            if len(block_new) >= 3:
                M, _ = cv2.estimateAffinePartial2D(block_old, block_new, method=cv2.RANSAC, ransacReprojThreshold=5.0)
                if M is not None:
                    transforms[(x, y)] = M
    return transforms


def apply_global_transform(img, M):
    """
    应用全局仿射变换矩阵到图像。

    Args:
        img (ndarray): 目标图像。
        M (ndarray): 仿射变换矩阵。

    Returns:
        ndarray: 变换后的图像。
    """
    if M is None:
        return img
    rows, cols, _ = img.shape
    warped = cv2.warpAffine(img, M, (cols, rows), flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REFLECT)
    return warped


def apply_local_transforms(img, transforms, grid_size=(8, 8), block_size=50):
    """
    应用分块的局部仿射变换矩阵到图像。

    Args:
        img (ndarray): 目标图像。
        transforms (dict): 每个块的变换矩阵。
        grid_size (tuple): 网格大小（x, y）。
        block_size (int): 每个块的大小。

    Returns:
        ndarray: 变换后的图像。
    """
    warped = np.zeros_like(img)
    for x in range(grid_size[0]):
        for y in range(grid_size[1]):
            M = transforms.get((x, y))
            if M is not None:
                # 定义块的区域
                start_x = x * block_size
                start_y = y * block_size
                end_x = (x + 1) * block_size
                end_y = (y + 1) * block_size
                # 提取块
                block = img[start_y:end_y, start_x:end_x]
                if block.shape[0] == 0 or block.shape[1] == 0:
                    continue
                # 应用变换
                warped_block = cv2.warpAffine(block, M, (block.shape[1], block.shape[0]), flags=cv2.INTER_LINEAR,
                                              borderMode=cv2.BORDER_REFLECT)
                warped[start_y:end_y, start_x:end_x] = warped_block
    return warped


def laplacian_pyramid_fusion(img1, img2, levels=5):
    """
    使用拉普拉斯金字塔方法融合两张图像。

    Args:
        img1 (ndarray): 参考图像。
        img2 (ndarray): 对齐后的图像。
        levels (int): 金字塔层数。

    Returns:
        ndarray: 融合后的图像。
    """
    # 构建高斯金字塔
    G1 = img1.copy()
    G2 = img2.copy()
    gp1 = [G1]
    gp2 = [G2]
    for i in range(levels):
        G1 = cv2.pyrDown(G1)
        G2 = cv2.pyrDown(G2)
        gp1.append(G1)
        gp2.append(G2)

    # 构建拉普拉斯金字塔
    lp1 = [gp1[levels - 1]]
    lp2 = [gp2[levels - 1]]
    for i in range(levels - 1, 0, -1):
        GE1 = cv2.pyrUp(gp1[i])
        GE2 = cv2.pyrUp(gp2[i])
        L1 = cv2.subtract(gp1[i - 1], GE1)
        L2 = cv2.subtract(gp2[i - 1], GE2)
        lp1.append(L1)
        lp2.append(L2)

    # 融合拉普拉斯金字塔
    LS = []
    for l1, l2 in zip(lp1, lp2):
        rows, cols, dpt = l1.shape
        ls = np.hstack((l1, l2))
        LS.append(ls)

    # 重建图像
    fused_image = LS[0]
    for i in range(1, levels):
        fused_image = cv2.pyrUp(fused_image)
        fused_image = cv2.add(fused_image, LS[i])

    return fused_image


def main():
    # 多帧DNG图像路径列表（请替换为实际的DNG文件路径）
    dng_files = ['./input_file/payload_N000.dng', './input_file/payload_N001.dng', './input_file/payload_N002.dng', './input_file/payload_N003.dng']


    # 读取所有DNG图像
    images_rgb = []
    images_gray = []
    for file in dng_files:
        rgb, gray = read_dng(file, scale=1.0)  # 根据需要调整scale
        images_rgb.append(rgb)
        images_gray.append(gray)

    # 设定参考帧为第一帧
    ref_rgb = images_rgb[0]
    ref_gray = images_gray[0]

    # 特征点检测参数
    feature_params = dict(maxCorners=500,
                          qualityLevel=0.3,
                          minDistance=7,
                          blockSize=7)

    # Lucas-Kanade光流参数
    lk_params = dict(winSize=(15, 15),
                     maxLevel=4,  # 使用更高的金字塔层数处理更大位移
                     criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))

    # 检测参考帧的特征点
    p0 = detect_features(ref_gray, **feature_params)

    if p0 is None:
        print("参考帧未检测到任何特征点。")
        return

    # 初始化融合图像列表
    fused_images = [ref_rgb.copy()]

    # 依次处理后续帧
    for i in range(1, len(images_rgb)):
        print(f"Processing frame {i + 1}/{len(images_rgb)}...")
        current_rgb = images_rgb[i]
        current_gray = images_gray[i]

        # 跟踪特征点
        good_new, good_old, st = track_features(ref_gray, current_gray, p0, lk_params)
        if good_new is None or good_old is None:
            print(f"No features tracked in frame {i + 1}. Skipping...")
            continue

        # 全局变换矩阵估计
        M_global, mask_global = estimate_global_transform(good_new, good_old)
        if M_global is None:
            print(f"Global transform estimation failed for frame {i + 1}. Skipping...")
            continue

        # 应用全局变换对齐图像
        aligned_img = apply_global_transform(current_rgb, M_global)

        # 局部变换矩阵估计（分块）
        transforms_local = estimate_local_transforms(good_new, good_old, grid_size=(8, 8), block_size=50)

        # 应用局部变换对齐图像
        aligned_img_local = apply_local_transforms(aligned_img, transforms_local, grid_size=(8, 8), block_size=50)

        # 图像融合（拉普拉斯金字塔融合）
        fused_img = laplacian_pyramid_fusion(ref_rgb, aligned_img, levels=5)

        # 更新参考帧为当前融合结果
        ref_rgb = fused_img.copy()
        ref_gray = cv2.cvtColor(ref_rgb, cv2.COLOR_RGB2GRAY)
        ref_gray = cv2.GaussianBlur(ref_gray, (5, 5), 0)

        # 重新检测特征点
        # p0 = detect_features(ref_gray, **feature_params)
        if p0 is None:
            print(f"No features detected in fused frame {i + 1}.")

        # 添加到融合图像列表
        fused_images.append(fused_img)

    # 显示最终融合结果
    plt.figure(figsize=(10, 10))
    plt.imshow(cv2.cvtColor(fused_images[-1], cv2.COLOR_RGB2BGR))
    plt.title('Fused Image')
    plt.axis('off')
    plt.show()

    # 保存最终融合结果
    cv2.imwrite('fused_image.png', cv2.cvtColor(fused_images[-1], cv2.COLOR_RGB2BGR))
    print("Fused image saved as 'fused_image.png'.")


if __name__ == "__main__":
    main()

五、优缺点分析

1. 优点

• 计算效率高：Lucas-Kanade方法基于局部线性假设，计算简单，适合实时应用。
• 实现简单：算法流程明确，易于实现，适合嵌入式系统和实时视频处理。
• 鲁棒性强：通过窗口内多个像素点的约束，提高了对噪声和误差的鲁棒性。
• 适用性广：在小位移和亮度变化较小的情况下，具有较高的准确性。

2. 缺点

• 小位移限制：对于大位移或快速运动的场景，单次估计可能不准确，需多次迭代或金字塔方法。

• 纹理依赖性：依赖于图像的梯度信息，对于纹理稀疏或均匀区域，光流估计不准确。

• 遮挡处理困难：难以处理复杂的遮挡情况，可能导致错误的运动估计。

• 假设限制：局部恒定光流和亮度恒定假设在实际场景中可能不完全成立，影响估计精度。

Lucas-Kanade方法作为光流估计的经典算法，凭借其简洁高效的特点在计算机视觉领域占据重要地位。尽管存在一些局限性，但通过各种改进和结合现代技术，Lucas-Kanade方法及其衍生算法在实际应用中依然表现出色。随着深度学习和计算能力的提升，未来光流估计将朝着更加精确、高效和智能化的方向发展，Lucas-Kanade方法仍将在其中发挥基础性作用。

参考文献：
经典光流算法Lucas-Kanade（有图助理解）