diff --git a/scripts/run.jl b/scripts/run.jl index ba63920a42adc6062675f298d6056d41c0331bcc..64f564b2261b44b4171d0b29c0d4dc8a3c894c53 100644 --- a/scripts/run.jl +++ b/scripts/run.jl @@ -7,6 +7,7 @@ import DataFrames: DataFrame import FileIO using OpticalFlowUtils using WriteVTK: paraview_collection +using Plots using SemiSmoothNewton using SemiSmoothNewton: HMesh, ncells, refine @@ -691,12 +692,13 @@ function optflow(ctx) # coord flip imgf0 = from_img(ctx.params.imgf0) imgf1 = from_img(ctx.params.imgf1) + maxflow = 4.6157 # Dimetrodon size(imgf0) == size(imgf1) || throw(ArgumentError("non-matching dimensions")) m = 2 #mesh = init_grid(imgf0; type=:vertex) - mesh = init_grid(imgf0, 1, 1) + mesh = init_grid(imgf0, (size(imgf0) .รท 16)...) #mesh = init_grid(imgf0) # optflow specific stuff @@ -706,25 +708,33 @@ function optflow(ctx) fw = FeFunction(Vg, name="fw") T(tdata, u) = tdata * u # tdata = nablafw - #S(u, nablau) = nablau - S(u, nablau) = u + S(u, nablau) = nablau + #S(u, nablau) = u st = L1L2TVContext("run", mesh, m; T, tdata = nabla(fw), S, ctx.params.alpha1, ctx.params.alpha2, ctx.params.lambda, ctx.params.beta, ctx.params.gamma1, ctx.params.gamma2) - function reproject!() - project_img!(f0, imgf0) - project_img!(f1, imgf1) + function warp!() + imgfw = warp_backwards(imgf1, sample(st.u)) + project_img!(fw, imgfw) - # TODO: should be warping from imgfw - fw.data .= f1.data + # replace new tdata + st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fw), S, + st.alpha1, st.alpha2, st.beta, st.lambda, st.gamma1, st.gamma2, + st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2) g_optflow(x; u, f0, fw, nablafw) = nablafw * u - (fw - f0) - interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = nabla(fw)) + interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = st.tdata) + end + + function reproject!() + project_img2!(f0, imgf0) + project_img2!(f1, imgf1) end reproject!() + warp!() save_step(i) = output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"), @@ -734,14 +744,26 @@ function optflow(ctx) pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd pvd[i] = save_step(i) while true - for k in 1:10 + for j in 1:4 + norm_g_old = norm_l2(st.g) + for k in 1:5 + i += 1 + step!(st) + estimate!(st) + pvd[i] = save_step(i) + println() + end + warp!() + norm_g = norm_l2(st.g) i += 1 - step!(st) - estimate!(st) pvd[i] = save_step(i) - println() + display(plot(colorflow(to_img(sample(st.u)); maxflow))) end + i >= 50 && break + #continue + marked_cells = mark(st; theta = 0.5) + #marked_cells = Set(axes(mesh.cells, 2)) println("refining ...") @@ -752,12 +774,16 @@ function optflow(ctx) st.alpha1, st.alpha2, st.beta, st.lambda, st.gamma1, st.gamma2, fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2) f0, f1, fw = (fs.f0, fs.f1, fs.fw) + i += 1 + pvd[i] = save_step(i) println("reprojecting ...") reproject!() + i += 1 + pvd[i] = save_step(i) end end - display(plot(colorflow(to_img(sample(st.u))))) + display(plot(colorflow(to_img(sample(st.u)); maxflow))) #CSV.write(joinpath(ctx.outdir, "energies.csv"), df) #saveimg(joinpath(ctx.outdir, "output_glob.png"), fetch_u(states.glob)) diff --git a/src/image.jl b/src/image.jl index a8bebd6a4cbd4865e999f546faa469742a3a2157..081d5e5eebc3c9ec10050e1818f5d45b2d0d0f37 100644 --- a/src/image.jl +++ b/src/image.jl @@ -1,4 +1,4 @@ -export interpolate_bilinear, halve +export interpolate_bilinear, halve, warp_backwards eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...]