模拟一个简单的光流

Question

我目前正在尝试使用以下等式模拟光流：

下面是一个基本示例，其中我有一个 7x7 图像，其中中心像素被照亮。我应用的速度是均匀的 x 速度 2。

using Interpolations
using PrettyTables

# Generate grid
nx = 7  # Image will be 7x7 pixels
x = zeros(nx, nx)
yy = repeat(1:nx, 1, nx)  # grid of y-values
xx = yy'  # grid of x-values

# In this example x is the image I in the above equation
x[(nx-1)÷2 + 1, (nx-1)÷2 + 1] = 1.0  # set central pixel equal to 1

# Initialize velocity
velocity = 2;
vx = velocity .* ones(nx, nx); # vx=2
vy = 0.0      .* ones(nx, nx); # vy=0

for t in 1:1
    # create 2d grid interpolator of the image
    itp = interpolate((collect(1:nx), collect(1:nx)), x, Gridded(Linear()));

    # create 2d grid interpolator of vx and vy
    itpvx = interpolate((collect(1:nx), collect(1:nx)), vx, Gridded(Linear()));
    itpvy = interpolate((collect(1:nx), collect(1:nx)), vy, Gridded(Linear()));

    ∇I_x  = Array{Float64}(undef, nx, nx);   # Initialize array for ∇I_x
    ∇I_y  = Array{Float64}(undef, nx, nx);   # Initialize array for ∇I_y
    ∇vx_x = Array{Float64}(undef, nx, nx);  # Initialize array for ∇vx_x
    ∇vy_y = Array{Float64}(undef, nx, nx);  # Initialize array for ∇vy_y
    for i=1:nx
        for j=1:nx
            # gradient of image in x and y directions
            Gx = Interpolations.gradient(itp, i, j);
            ∇I_x[i, j] = Gx[2];
            ∇I_y[i, j] = Gx[1];
            Gvx = Interpolations.gradient(itpvx, i, j)  # gradient of vx in both directions
            Gvy = Interpolations.gradient(itpvy, i, j)  # gradient of vy in both directions
            ∇vx_x[i, j] = Gvx[2];
            ∇vy_y[i, j] = Gvy[1];
        end
    end

    v∇I = (vx .* ∇I_x) .+ (vy .* ∇I_y)  # v dot ∇I
    I∇v = x .* (∇vx_x .+ ∇vy_y) # I dot ∇v
    x = x .- (v∇I .+ I∇v)  # I(x, y, t+dt)
    pretty_table(x)
end

我期望 x 中的发光像素会在 x_predicted 中向右移动两个像素。我看到的是以下内容：

原始照明像素的值两次移动到相邻像素，而不是向右移动两个像素。 IE。相邻像素从 0 变为 2，原始像素从 1 变为 -1。我不确定我是否在搞乱方程式，或者我在这里以错误的方式考虑速度。有什么想法吗？

Answer 1

虽然没有深入研究，但我认为这里有几个潜在的问题：

违反库朗条件

您最初发布的代码（我现在已经对其进行了编辑）模拟了单个时间步长。我不希望在单个时间步内激活距离源单元格 2 个单位的单元格。这样做会破坏 Courant condition。来自维基百科：

The principle behind the condition is that, for example, if a wave is moving across a discrete spatial grid and we want to compute its amplitude at discrete time steps of equal duration, then this duration must be less than the time for the wave to travel to adjacent grid points.

Courant 条件要求 uΔt/Δx <= 1（对于显式时间推进求解器，例如您已实现的求解器）。代入u=2，Δt=1，Δx=1得到2，大于1，所以你有一道数学题。解决这个问题的一般方法是使 Δt 变小。你可能想要这样的东西：

x = x .- Δt*(v∇I .+ I∇v)  # I(x, y, t+dt)

缺少梯度？

我有点担心这里发生的事情：

Gvx = Interpolations.gradient(itpvx, i, j)  # gradient of vx in both directions
Gvy = Interpolations.gradient(itpvy, i, j)  # gradient of vy in both directions
∇vx_x[i, j] = Gvx[2];
∇vy_y[i, j] = Gvy[1];

您可以从 Gvx 和 Gvy 中提取两个渐变，但您只使用其中的一个。这是否意味着您正在丢弃信息？

https://scicomp.stackexchange.com/ 可能会为此提供更好的帮助。

模拟一个简单的光流

Simulating a simple optical flow

image-processing

julia