从光流计算散度和旋度并绘制它

Question

我正在使用 flow = cv2.calcOpticalFlowFarneback() 来计算视频中的光流，它为我提供了一个形状为 (height, width) 的 numpy 数组, 2) 包含每个像素的 Fx 和 Fy 值 (flow[:,:,0] = Fx 和 flow[:,:,1] = Fy).

为了计算分歧，我使用 np.gradient 像这样：

def divergence_npgrad(flow):
    Fx, Fy = flow[:, :, 0], flow[:, :, 1]
    F = [Fx, Fy]
    d = len(F)
    return np.ufunc.reduce(np.add, [np.gradient(F[i], axis=i) for i in range(d)])

接下来，我要计算旋度。我知道 sympy.physics.vector 中有一个 curl 函数，但我真的不明白它是如何工作的，或者它如何应用于我的 flow。所以我想我也可以为此使用 np.gradient 。在 2D 中，我需要为每个像素计算 dFy/dx - dFx/dy，所以我会像这样：

def curl_npgrad(flow):
    Fx, Fy = flow[:, :, 0], flow[:, :, 1]
    dFx_dy = np.gradient(Fx, axis=1)
    dFy_dx = np.gradient(Fy, axis=0)
    curl = np.ufunc.reduce(np.subtract, [dFy_dx, dFx_dy])
    return curl

这样做是正确的方法还是我遗漏了什么？

现在，如果我有 curl，我想用 matplotlib 绘制两个图。我的意思是我想用不同的颜色图显示来自 flow 的向量。一个图将使用矢量的幅度值作为颜色图，归一化为 (0-magnitude_max)。另一个图将使用 curl 值作为颜色图，其中如果 curl 为负则箭头为蓝色，如果 则箭头为红色curl 在那个位置是正的。

这是我正在尝试使用的：

def flow_plot(flow, frame):
    frame = cv2.cvtColor(frame, cv2.COLOR_BGR2RGB)
    h, w = flow.shape[:2]
    dpi = 72
    xinch = w / dpi
    yinch = h / dpi

    step = 24

    y, x = np.mgrid[step / ((h % step) // 2):h:step, step / ((w % step) // 2):w:step].reshape(2, -1).astype(int)
    fx, fy = flow[y, x].T
    mag = np.sqrt(np.power(fx, 2) + np.power(fy, 2))
    fx = fx / mag
    fy = fy / mag

    curl = curl_npgrad(flow)
    curl_map = curl[y, x]

    quiver_params = dict(cmap='Oranges', # for magnitude
                         #cmap='seismic', # for curl
                         norm=colors.Normalize(vmin=0.0, vmax=1.0), # for magnitude
                         #norm = colors.CenteredNorm(), # for curl
                         units='xy',
                         scale=0.03,
                         headwidth=3,
                         headlength=5,
                         minshaft=1.5,
                         pivot='middle')

    fig = plt.figure(figsize=(xinch, yinch), dpi=dpi)
    plt.imshow(frame)
    plt.quiver(x, y, fx, fy, mag, **quiver_params)
    plt.gca().invert_yaxis()
    plt.gca().set_aspect('equal', 'datalim')
    plt.axis('off')
    fig.tight_layout(pad=0)
    fig.canvas.draw()
    img = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')
    img = img.reshape(fig.canvas.get_width_height()[::-1] + (3,))
    img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)
    img = cv2.flip(img, 0)
    plt.close(fig)
    return img

我正在将绘图转换为 cv2 图像，以便我可以将其用于 opencv 视频编写器。

我注意到如果我不显示情节背后的原始框架，我必须反转 y 轴并在 [= 中使用 -fy 77=]()，如果我想显示后面的框架，我也必须反转y轴，可以使用fy，但我必须翻转之后的整个图像。它有什么意义？听不懂。

至于卷曲，对我来说有点乱。几乎没有显示任何颜色，随机的红色和蓝色斑点，而不是一堆 red/blue 箭头，其中流体明显旋转。就像这些： image1 of the messy curl, image2 of the messy curl

计算这种流的旋度是不是一种不好的方法？我错过了什么？

Answer 1

如果我理解你正确设置坐标轴的方式，你会在 divergence_npgrad 和 curl_npgrad 的顶部缺少一个 flow = np.swapaxes(flow, 0, 1)。

我尝试将你的旋度和散度函数应用于我已经知道正确旋度和散度的简单函数。

例如，我试过这个函数：

F_x = x
F_y = y

产生这个情节：

对于这个函数，处处散度都应该很高。应用你的函数，它告诉我散度为 0。

生成此数组的代码：

import numpy as np
import matplotlib.pyplot as plt

x,y = np.meshgrid(np.linspace(-2,2,10),np.linspace(-2,2,10))

u = x
v = y
field2 = np.stack((u, v), axis=-1)

plt.quiver(x, y, field2[:, :, 0], field2[:, :, 1])

我也尝试了一个测试curl的函数：

F_x = cos(x + y)
F_y = sin(x - y)

产生这个情节的是：

生成此数组的代码：

import numpy as np
import matplotlib.pyplot as plt

x,y = np.meshgrid(np.linspace(0,2,10),np.linspace(0,2,10))

u = np.cos(x + y)
v = np.sin(x - y)
field = np.stack((u, v), axis=-1)

plt.quiver(x, y, field[:, :, 0], field[:, :, 1])

_{感谢U of Mich. Math department for this example!}

对于此函数，旋涡中心周围的旋度应该最高。将你的 curl 函数应用于此，我得到的结果是 curl 在角落最高。

这是我试过的有效代码：

def divergence_npgrad(flow):
    flow = np.swapaxes(flow, 0, 1)
    Fx, Fy = flow[:, :, 0], flow[:, :, 1]
    dFx_dx = np.gradient(Fx, axis=0)
    dFy_dy = np.gradient(Fy, axis=1)
    return dFx_dx + dFy_dy

def curl_npgrad(flow):
    flow = np.swapaxes(flow, 0, 1)
    Fx, Fy = flow[:, :, 0], flow[:, :, 1]
    dFx_dy = np.gradient(Fx, axis=1)
    dFy_dx = np.gradient(Fy, axis=0)
    curl = dFy_dx - dFx_dy
    return curl

我认为这是正确更改的一些解释：

np.gradient(..., axis=0) 提供跨行的渐变，因为那是轴 0。在您的输入图像中，您具有形状 (height, width, 2)，因此轴 0 实际上表示高度或 y。使用 np.swapaxes(flow, 0, 1) 交换该数组的 x 轴和 y 轴的顺序。
不需要使用 np.ufunc.reduce - 您可以改用广播，它会正常工作。它也更容易阅读。

以下是使用此代码的结果。

这里是第一个函数的散度计算结果

结果：到处都是正值，恒定值。
这里是第二个函数的curl计算结果

结果：以函数的旋涡部分为中心的正值。（这里有一个轴的变化 - 0,0 在左上角。）

从光流计算散度和旋度并绘制它

Calculating divergence and curl from optical flow and plotting it

python

opencv

numpy

matplotlib