根据光流计算发散度和旋度并将其绘制出来

Ben*_*agy 4 python opencv numpy matplotlib

我使用flow = cv2.calcOpticalFlowFarneback()来计算视频中的光流,它为我提供了一个形状为 (height, width, 2) 的 numpy 数组,其中包含每个像素的FxFy值 ( flow[: ,:,0] = Fx流量[:,:,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)])
Run Code Online (Sandbox Code Playgroud)

接下来,我要计算旋度。我知道sympy.physicals.vector中有一个curl函数,但我真的不明白它是如何工作的或者它如何应用于我的流程。所以我想我也可以使用 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
Run Code Online (Sandbox Code Playgroud)

这是正确的方法还是我错过了什么?

现在,如果我有curl,我想用matplotlib绘制两个图。我的观点是,我想用不同的颜色图显示中的向量。一个图将使用向量的幅度值作为颜色图,标准化为 (0-magnitude_max)。另一个图将使用旋度值作为颜色图,其中如果该位置的旋度为负,则箭头为蓝色,如果旋度为正,则箭头为红色。

这是我正在尝试使用的:

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
Run Code Online (Sandbox Code Playgroud)

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

我注意到,如果我没有显示绘图后面的原始框架,我必须反转 y 轴并在plt.quiver()中使用-fy,如果我想显示后面的框架,我必须反转 y 轴也可以使用fy,但之后我必须翻转整个图像。这有什么意义呢?我听不懂。

至于卷发,对我来说有点凌乱。几乎没有显示任何颜色,随机的红色和蓝色斑点,没有一堆流体清晰旋转的红色/蓝色箭头。就像这样: 凌乱卷发的 image1凌乱卷发的 image2

计算这种流量的旋度是不是一个不好的方法?我缺少什么?

Nic*_*ell 5

如果我理解您正确设置轴的方式,那么您在和flow = np.swapaxes(flow, 0, 1)的顶部都缺少 a 。divergence_npgradcurl_npgrad

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

例如,我尝试了这个功能:

F_x = x
F_y = y
Run Code Online (Sandbox Code Playgroud)

产生这个情节:

箭头从中心指向的绘图

对于这个函数,散度在任何地方都应该很高。应用你的函数,它告诉我散度是 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])
Run Code Online (Sandbox Code Playgroud)

我还尝试了一个测试curl的函数:

F_x = cos(x + y)
F_y = sin(x - y)
Run Code Online (Sandbox Code Playgroud)

这产生了这个情节:

箭头绕中心旋转的绘图

生成此数组的代码:

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])
Run Code Online (Sandbox Code Playgroud)

感谢密歇根大学数学系提供的这个例子

对于此功能,卷度应在中心漩涡周围最高。将你的卷曲函数应用于此,我得到的结果是卷曲在角落最高。

这是我尝试过的有效代码:

F_x = x
F_y = y
Run Code Online (Sandbox Code Playgroud)

为什么我认为这是正确的更改的一些解释:

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

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

  1. 这是第一个函数的散度计算结果。

    x 和 y 函数的图,到处都是 0.888

    结果:到处都是积极的、恒定的价值。

  2. 这是第二个函数的旋度计算结果。

    函数在 x=y=pi/4 处达到峰值并从那里开始下降

    结果:以函数的漩涡部分为中心的正值。(这里的轴发生了变化 - 0,0 位于左上角。)