为什么输出只包含 2 个值而不包含整个图像的位移？

Question

我已经被困在这里一段时间了。我不明白我在使用Lucas Kanade 方法计算沿 x 轴和 y 轴的位移向量时做错了什么。

我按照上面的维基百科链接实现了它。这是我所做的：

    import cv2
    import numpy as np


    img_a = cv2.imread("./images/1.png",0)
    img_b = cv2.imread("./images/2.png",0)


    # Calculate gradient along x and y axis
    ix = cv2.Sobel(img_a, cv2.CV_64F, 1, 0, ksize = 3, scale = 1.0/3.0)
    iy = cv2.Sobel(img_a, cv2.CV_64F, 0, 1, ksize = 3, scale = 1.0/3.0)

    # Calculate temporal difference between the 2 images
    it = img_b - img_a


    ix = ix.flatten()
    iy = iy.flatten()
    it = -it.flatten()

    A = np.vstack((ix, iy)).T


    atai = np.linalg.inv(np.dot(A.T,A))
    atb = np.dot(A.T, it)

    v = np.dot(np.dot(np.linalg.inv(np.dot(A.T,A)),A.T),it)

    print(v)

这段代码运行没有错误，但它打印了一个包含 2 个值的数组！我曾预计v 矩阵的大小与图像的大小相同。为什么会这样？我做错了什么？

PS：我知道 OpenCV 有直接可用的方法，但我想自己编写这个简单的算法（也可以在上面共享的 Wikipedia 链接中给出）。

Answer 1

要正确计算 Lucas–Kanade 光流估计，您需要使用来自其邻域的信息，而不是整个图像，为每个像素求解两个方程组。

这是配方（符号指的是the Wikipedia page上使用的）：

1. 使用任何方法计算第一张图像（OP 中的ix、iy）的图像梯度（A）（Sobel 可以，我更喜欢高斯导数；注意它是在 Sobel 中应用正确的缩放比例很重要：1/8）。

    ix = cv2.Sobel(img_a, cv2.CV_64F, 1, 0, ksize = 3, scale = 1.0/8.0)
    iy = cv2.Sobel(img_a, cv2.CV_64F, 0, 1, ksize = 3, scale = 1.0/8.0)
   

2. 计算结构张量（A^TWA）：Axx = ix * ix、Axy = ix * iy、Ayy = iy * iy。这三个图像中的每一个都必须使用高斯滤波器进行平滑（这是窗口化）。例如，

   Axx = cv2.GaussianBlur(ix * ix, (0,0), 5)
   Axy = cv2.GaussianBlur(ix * iy, (0,0), 5)
   Ayy = cv2.GaussianBlur(iy * iy, (0,0), 5)
   

   这三幅图像共同构成了结构张量，它是每个像素处的 2x2 对称矩阵。对于(i,j) 处的像素，矩阵为：

   |  Axx(i,j)  Axy(i,j)  |
   |  Axy(i,j)  Ayy(i,j)  |
   

3. 通过减去两个图像（OP 中的it）来计算时间梯度（b）。

   it = img_b - img_a
   

4. 计算A^TWb：Abx = ix * it、Aby = iy * it，并使用与上述相同的高斯滤波器平滑这两个图像。

   Abx = cv2.GaussianBlur(ix * it, (0,0), 5)
   Aby = cv2.GaussianBlur(iy * it, (0,0), 5)
   

5. 计算A^TWA（对称正定矩阵）的逆并乘以A^TWb时间>。请注意，这个逆矩阵是每个像素的 2x2 矩阵，而不是整个图像。您可以将其写为图像Axx、Axy、Ayy、Abx 和Aby 上的一组简单算术运算。

   矩阵A^TWA的逆由下式给出：

   |  Ayy -Axy  |  
   | -Axy  Axx  | / ( Axx*Ayy - Axy*Axy )
   

   所以你可以把解决方案写成

   norm = Axx*Ayy - Axy*Axy
   vx = ( Ayy * Abx - Axy * Aby ) / norm
   vy = ( Axx * Aby - Axy * Abx ) / norm
   

   如果图像是自然的，它至少会有一点点噪点，norm 不会有零。但是对于人造图像norm 可能有零，这意味着你不能除以它。只需向其添加一个小值即可避免除以零错误：norm += 1e-6。

选择高斯过滤器的大小作为精度和允许的运动速度之间的折衷：较大的过滤器将产生不太精确的结果，但适用于图像之间较大的偏移。

通常，vx 和 vy 仅在矩阵 A^TWA 的两个特征值足够大（如果至少一个很小）的情况下进行评估，结果不准确或可能错误）。

---

使用DIPlib（披露：我是作者）这一切都非常简单，因为它支持在每个像素处具有矩阵的图像。你可以这样做：

import diplib as dip

img_a = dip.ImageRead("./images/1.png")
img_b = dip.ImageRead("./images/2.png")

A = dip.Gradient(img_a, [1.0])
b = img_b - img_a
ATA = dip.Gauss(A * dip.Transpose(A), [5.0])
ATb = dip.Gauss(A * b, [5.0])
v = dip.Inverse(ATA) * ATb