I am currently trying to simulate an optical flow using the following equation:

Below is a basic example where I have a 7x7 image where the central pixel is illuminated. The velocity I am applying is a uniform x-velocity of 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

What I expect is that the illuminated pixel in x will shift two pixels to the right in x_predicted. What I am seeing is the following:

where the original illuminated pixel's value is moved to the neighboring pixel twice rather than being shifted two pixels to the right. I.e. the neighboring pixel goes from being 0 to 2 and the original pixel goes from a value of 1 to -1. I'm not sure if I'm messing up the equation or if I'm thinking of velocity in the wrong way here. Any ideas?

CodePudding user response：

Without looking into it too deeply, I think there are a couple of potential issues here:

Violation of the Courant Condition

The code you originally posted (I've edited it now) simulates a single timestep. I would not expect a cell 2 units away from your source cell to be activated in a single timestep. Doing so would voilate the Courant condition. From wikipedia:

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.

The Courant condition requires that uΔt/Δx <= 1 (for an explicit time-marching solver such as the one you've implemented). Plugging in u=2, Δt=1, Δx=1 gives 2, which is greater than 1, so you have a mathematical problem. The general way of fixing this problem is to make Δt smaller. You probably want something like:

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

Missing gradients?

I'm a little concerned about what's going on here:

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];

You're able to pull two gradients out of both Gvx and Gvy, but you're only using one from each of them. Does that mean you're throwing information away?

https://scicomp.stackexchange.com/ is likely to provide better help with this.