From 8beee843bdf599dc103c13194c6cfc3cbf00bd27 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Sat, 15 Jan 2022 17:39:08 +0100
Subject: [PATCH] implement accumulated warping
---
scripts/run_experiments.jl | 63 +++++++++++++++++++++-----------------
1 file changed, 35 insertions(+), 28 deletions(-)
diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 93ccb9d..5bbf6eb 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -851,18 +851,19 @@ function inpaint(img, imgmask; name, params...)
end
function optflow(ctx)
- tol_newton = 1e-6 # cauchy criterion
-
size(ctx.params.imgf0) == size(ctx.params.imgf1) ||
throw(ArgumentError("non-matching image sizes"))
+ tol_newton = 1e-5 # cauchy criterion for inner newton loop
+
m = 2 # range dim of vector field
+
+ # convert to cartesian coordinates
imgf0 = from_img(ctx.params.imgf0)
imgf1 = from_img(ctx.params.imgf1)
- #mesh = init_grid(imgf0)
- mesh = init_grid(imgf0, (size(imgf0) .÷ 2^3)...)
- domain_factor = 1 / sqrt(area(mesh))
+ mesh = init_grid(imgf0, (size(imgf0) .÷ 2^2)...)
+ mesh_area = area(mesh)
# optflow specific stuff
Vg = FeSpace(mesh, P1(), (1,))
@@ -872,19 +873,21 @@ function optflow(ctx)
T(tdata, u) = tdata * u # here tdata will be nabla(fw)
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)
+ u_acc = FeFunction(st.u.space, name="u_acc")
+ u_acc.data .= 0
+
function warp!()
- imgfw = warp_backwards(imgf1, sample(st.u))
- project_l2_lagrange!(fw, imgfw)
+ imgfw = warp_backwards(imgf1, sample(u_acc))
+ project_l2_pixel!(fw, imgfw)
g_optflow(x; u, f0, fw, nablafw) =
- #nablafw * u
- nablafw * u - (fw - f0)
+ -(fw - f0)
+ #nablafw * u - (fw - f0)
interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = st.tdata)
end
@@ -898,7 +901,7 @@ function optflow(ctx)
save_step(i) =
output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
- st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw)
+ st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw, u_acc)
i = 0
pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
@@ -917,34 +920,38 @@ function optflow(ctx)
#estimate!(st)
pvd[i] = save_step(i)
- norm_step_ = norm_step(st) * domain_factor
+ norm_step_ = norm_step(st) / sqrt(mesh_area)
println("norm_step = $norm_step_")
# residual is basically unusable: takes long to compute and
# oscillates
- #norm_residual_ = norm_residual(st) * domain_factor
+ #norm_residual_ = norm_residual(st) / sqrt(mesh_area)
#println("norm_residual = $norm_residual_")
- display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
+ #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
end
+
+ u_acc.data .+= st.u.data
+ display(plot(colorflow(to_img(sample(u_acc)); ctx.params.maxflow)))
+
# poor man's hash
#println(integrate(st.mesh, (x; tdata) -> tdata; st.tdata))
- println("refine ...")
- marked_cells = Set(axes(mesh.cells, 2))
- mesh, fs = refine(mesh, marked_cells;
- st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2,
- f0, f1, fw)
- st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fs.fw), S,
- 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)
+ #println("refine ...")
+ #marked_cells = Set(axes(mesh.cells, 2))
+ #mesh, fs = refine(mesh, marked_cells;
+ # st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2,
+ # f0, f1, fw, u_acc)
+ #st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fs.fw), S,
+ # 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, u_acc = (fs.f0, fs.f1, fs.fw, fs.u_acc)
- #i += 1
- pvd[i] = save_step(i)
+ ##i += 1
+ #pvd[i] = save_step(i)
- println("interpolate image data ...")
- interpolate_image_data!()
+ #println("interpolate image data ...")
+ #interpolate_image_data!()
println("warp ...")
warp!() # yay
@@ -978,7 +985,7 @@ function optflow(ctx)
# i += 1
# pvd[i] = save_step(i)
end
- display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
+ #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
#CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
#saveimg(joinpath(ctx.outdir, "output_glob.png"), fetch_u(states.glob))
--
GitLab