我目前正在嘗試使用以下等式模擬光流：
下面是一個基本示例，其中我有一個 7x7 影像，其中中心像素被照亮。我應用的速度是 2 的均勻 x 速度。
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。我所看到的是以下內容：
其中原始照明像素的值被移動到相鄰像素兩次，而不是向右移動兩個像素。即相鄰像素從 0 變為 2，原始像素從 1 的值變為 -1。我不確定我是否搞砸了等式，或者我是否在這里以錯誤的方式考慮速度。有任何想法嗎？
uj5u.com熱心網友回復：
在沒有深入研究的情況下，我認為這里有幾個潛在的問題：
違反柯朗條件
您最初發布的代碼（我現在已經對其進行了編輯）模擬了一個時間步長。我不希望在單個時間步長內激活距源單元格 2 個單位的單元格。這樣做會破壞Courant 條件。來自維基百科：
該條件背后的原理是，例如，如果一個波在一個離散的空間網格上移動，我們想在相等持續時間的離散時間步長上計算它的振幅，那么這個持續時間必須小于波傳播的時間到相鄰的網格點。
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/可能會為此提供更好的幫助。
轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/373242.html