分层Lucas-Kanade光流算法实现异常问题求助

问题现象 #

• 应用于旋转球体数据集时，光流x方向表现正确，但y方向结果异常，怀疑可视化代码存在问题。
• 静止像素区域出现非零光流值。

相关代码 #

可视化函数 #

def show_flow(img, flow, filename=None):
    x = np.arange(0, img.shape[1], 1)
    y = np.arange(0, img.shape[0], 1)
    x, y = np.meshgrid(x, y)
    plt.figure(figsize=(10, 10))
    fig = plt.imshow(img, cmap="gray", interpolation="bicubic")
    plt.axis("off")
    fig.axes.get_xaxis().set_visible(False)
    fig.axes.get_yaxis().set_visible(False)

    num_points_per_axis = 32
    step = int(img.shape[0] / num_points_per_axis)
    plt.quiver(
        x[::step, ::step],
        y[::step, ::step],
        flow[::step, ::step, 0],
        flow[::step, ::step, 1],
        color="r",
        pivot="tail",
        headwidth=2,
        headlength=3,
    )
    if filename is not None:
        plt.savefig(filename, bbox_inches="tight", pad_inches=0)

完整算法实现 #

基础Lucas-Kanade实现

def lucas_kanade(img1, img2):
    img1 = np.copy(img1).astype(np.float32)
    img2 = np.copy(img2).astype(np.float32)

    window_size = min(max(3, min(img1.shape[:2]) / 6), 31)
    window_size = int(2 * (window_size // 2) + 1)
    print("window size: ", window_size)

    Ix = np.zeros(img1.shape, dtype=np.float32)
    Iy = np.zeros(img1.shape, dtype=np.float32)
    Ix[1:-1, 1:-1] = (img1[1:-1, 2:] - img1[1:-1, :-2]) / 2
    Iy[1:-1, 1:-1] = (img1[2:, 1:-1] - img1[:-2, 1:-1]) / 2

    It = np.zeros(img1.shape, dtype=np.float32)
    It = img1 - img2

    kernel = np.ones((window_size, window_size), dtype=np.float32)

    Ix2 = convolve2d(Ix**2, kernel, mode="same", boundary="fill", fillvalue=0)
    Iy2 = convolve2d(Iy**2, kernel, mode="same", boundary="fill", fillvalue=0)
    Ixy = convolve2d(Ix * Iy, kernel, mode="same", boundary="fill", fillvalue=0)
    Ixt = convolve2d(Ix * It, kernel, mode="same", boundary="fill", fillvalue=0)
    Iyt = convolve2d(Iy * It, kernel, mode="same", boundary="fill", fillvalue=0)

    det = Ix2 * Iy2 - Ixy**2
    u = np.where((det > 1e-6), (Iy2 * Ixt - Ixy * Iyt) / det, 0)
    v = np.where((det > 1e-6), (Ix2 * Iyt - Ixy * Ixt) / det, 0)

    optical_flow = np.stack((u, v), axis=2)
    return optical_flow.astype(np.float32)

生成高斯金字塔

def gen_gaussian_pyramid(im, max_level):
    gauss_pyr = [im]
    for i in range(max_level):
        gauss_pyr.append(cv2.pyrDown(gauss_pyr[-1]))
    return gauss_pyr

上采样光流

def expand(img, dst_size, interpolation=None):
    height, width = dst_size[:2]
    return cv2.GaussianBlur(
        cv2.resize(
            img, (width, height), interpolation=interpolation or cv2.INTER_LINEAR
        ),
        (5, 5),
        0,
    )

图像扭曲

def remap(a, flow):
    height, width = flow.shape[:2]

    y, x = np.meshgrid(np.arange(height), np.arange(width), indexing="ij")

    flow_map = np.column_stack(
        (x.flatten() + -flow[:, :, 0].flatten(), y.flatten() + -flow[:, :, 1].flatten())
    )

    flow_map = flow_map.reshape((height, width, 2))
    flow_map[:, :, 0] = np.clip(flow_map[:, :, 0], 0, width - 1)
    flow_map[:, :, 1] = np.clip(flow_map[:, :, 1], 0, height - 1)
    flow_map = flow_map.astype(np.float32)

    warped = cv2.remap(a, flow_map, None, cv2.INTER_LINEAR)
    return warped

分层Lucas-Kanade主函数

def hierarchical_lucas_kanade(im1, im2, max_level):
    gauss_pyr_1 = gen_gaussian_pyramid(im1, max_level)
    gauss_pyr_2 = gen_gaussian_pyramid(im2, max_level)

    g_L = [0 for _ in range(max_level + 1)]
    d_L = [0 for _ in range(max_level + 1)]
    assert len(g_L) == 5
    g_L[max_level] = np.zeros(gauss_pyr_1[-1].shape[:2] + (2,)).astype(np.float32)

    for level in range(max_level, -1, -1):
        warped = remap(gauss_pyr_1[level], g_L[level])
        d_L[level] = lucas_kanade(warped, gauss_pyr_2[level])

        g_L[level - 1] = 2.0 * expand(
            g_L[level] + d_L[level],
            gauss_pyr_2[level - 1].shape[:2] + (2,),
            interpolation=cv2.INTER_LINEAR,
        )
    return g_L[0] + d_L[0]

测试代码

sphere_seq = []
for fname in natsorted(Path("./input/sphere/").rglob("*.ppm")):
    sphere_seq.append(cv2.imread(str(fname), cv2.IMREAD_GRAYSCALE))

flows = []
for i in range(len(sphere_seq) - 1):
    flows.append(hierarchical_lucas_kanade(sphere_seq[i], sphere_seq[i + 1], max_level=4))
    show_flow(sphere_seq[i], flows[i], f"./output/sphere/flow-{i}.png")

问题排查与解决建议 #

针对y方向光流异常 #

1. 坐标系匹配问题：Matplotlib的quiver函数使用的坐标系中y轴向上，而OpenCV图像的y轴向下，两者方向相反。可视化时需将y方向光流分量取反，修改quiver调用中的v参数：

   plt.quiver(
       x[::step, ::step],
       y[::step, ::step],
       flow[::step, ::step, 0],
       -flow[::step, ::step, 1],  # 取反y分量
       color="r",
       pivot="tail",
       headwidth=2,
       headlength=3,
   )
   

2. 验证光流计算正确性：选取几个已知运动方向的点，打印其光流的y分量值，对比可视化箭头方向，确认是计算错误还是显示错误。
3. 检查扭曲函数的坐标映射：remap函数中，光流(u,v)表示像素从img1到img2的偏移，因此扭曲img1到img2的位置时，应该是x + u、y + v，当前代码中用了x + -flow[:, :, 0]，需确认光流定义和扭曲逻辑是否一致。

针对静止区域非零光流 #

1. 修正时间梯度It的符号：光流基本方程为Ix*u + Iy*v + It = 0，其中It = img2 - img1（图像随时间的变化）。当前代码中It = img1 - img2符号完全相反，这会导致静止区域的It不为零，进而计算出错误的非零光流。修正如下：

   It = img2 - img1
   

2. 提高行列式阈值：当前det > 1e-6的阈值过低，对于纹理稀疏的静止区域，梯度矩阵的行列式很小，这些区域的光流不可靠，应跳过。可将阈值提高至1e-4甚至1e-3：

   u = np.where((det > 1e-4), (Iy2 * Ixt - Ixy * Iyt) / det, 0)
   v = np.where((det > 1e-4), (Ix2 * Iyt - Ixy * Ixt) / det, 0)
   

3. 光流后处理：对计算得到的光流进行幅值过滤，将幅值小于极小值（如0.1）的光流设为0，去除静止区域的微小噪声：

   mag = np.sqrt(flow[...,0]**2 + flow[...,1]**2)
   flow[mag < 0.1] = 0
   

4. 优化窗口大小计算：当前窗口大小在金字塔底层（小分辨率）可能过大，导致过度平滑。可调整窗口大小的计算逻辑，比如根据金字塔层级动态调整，避免小图像上窗口过大。

内容的提问来源于stack exchange，提问作者Der Fänger im Roggen