From cf4aa0e632071627f798af139bbe2b843d0ed1a2 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Sun, 30 Jan 2022 16:38:02 +0100
Subject: [PATCH] refactor optflow examples
---
scripts/run_experiments.jl | 98 +++++++++++++++++++++-----------------
1 file changed, 54 insertions(+), 44 deletions(-)
diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 477c081..f2e80e9 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -1051,6 +1051,7 @@ function optflow(ctx)
size(ctx.params.imgf0) == size(ctx.params.imgf1) ||
throw(ArgumentError("non-matching image sizes"))
+ project_image! = project_l2_lagrange!
tol_newton = 1e-5 # cauchy criterion for inner newton loop
m = 2 # range dim of vector field
@@ -1073,14 +1074,15 @@ function optflow(ctx)
ctx.params.gamma1, ctx.params.gamma2)
function warp!()
+ println("warp and reproject ...")
# warp image into imgfw / fw
imgfw = warp_backwards(imgf1, sample(st.u))
- project_l2_pixel!(fw, imgfw)
+ project_image!(fw, imgfw)
# recompute optflow operator T based on u0 and fw
- # TODO: investigate what julia does for singular matrices here
function tdata_optflow(x; u0_deriv, nablafw)
- res = nablafw / (I + u0_deriv')
+ #res = nablafw / (I + u0_deriv')
+ res = nablafw
all(isfinite, res) || throw(DivideError("singular optflow matrix"))
return res
end
@@ -1088,15 +1090,17 @@ function optflow(ctx)
u0_deriv = nabla(st.u), nablafw = nabla(fw))
# recompute optflow data g
+ # note that it is important here that the space of g is general enough
+ # to avoid any information loss that could screw up image warping
g_optflow(x; u0, f0, fw, tdata) =
- #-(fw - f0)
st.T(tdata, u0) - (fw - f0)
interpolate!(st.g, g_optflow; u0 = st.u, f0, fw, st.tdata)
end
function interpolate_image_data!()
- project_l2_pixel!(f0, imgf0)
- project_l2_pixel!(f1, imgf1)
+ println("interpolate image data ...")
+ project_image!(f0, imgf0)
+ project_image!(f1, imgf1)
end
@@ -1104,55 +1108,59 @@ function optflow(ctx)
output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw)
- i = 0
pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
- pvd[i] = save_step(i)
+ interpolate_image_data!()
+ warp!() # in the first step just to fill st.g
+ pvd[0] = save_step(0)
- #norm_g_old = norm_l2(st.g)
-
- for j in 1:5
- println("interpolate image data ...")
- interpolate_image_data!()
- println("warp ...")
- warp!() # in the first step just to fill st.g
- i += 1
+ i = 0
+ k = 0
+ while true
+ step!(st)
+ k += 1
+ #estimate_pd!(st)
- # interior newton loop
- norm_step_ = Inf
- k = 0
- while norm_step_ > tol_newton && k < 20
- k += 1
- println()
- step!(st)
- #estimate_pd!(st)
+ norm_step_ = norm_step(st) / sqrt(mesh_area)
+ println("norm_step = $norm_step_")
- 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) / sqrt(mesh_area)
+ #println("norm_residual = $norm_residual_")
- # residual is basically unusable: takes long to compute and
- # oscillates
- #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
+ # interior newton stop criterion
+ norm_step_ > tol_newton && k < 10 && continue
+ k = 0
display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
- estimate_res!(st)
+ i += 1
pvd[i] = save_step(i)
- println("refine ...")
- marked_cells = mark(st; theta = 0.5)
- #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,
- st.tdata, f0, f1, fw)
- st = L1L2TVState(mesh, st.d, st.m, st.T, fs.tdata, st.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)
+ # warping stop criterion
+ warp!()
+ if i > 10
+ break
+ end
+
+ # refinement stop criterion
+ if false
+ println("refine ...")
+ estimate_res!(st)
+ marked_cells = mark(st; theta = 0.5)
+ #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,
+ st.tdata, f0, f1, fw)
+ st = L1L2TVState(mesh, st.d, st.m, st.T, fs.tdata, st.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)
+ interpolate_image_data!()
+ end
end
end
#display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
@@ -1167,9 +1175,11 @@ end
function experiment_optflow_middlebury(ctx)
imgf0 = loadimg(joinpath(ctx.indir, "frame10.png"))
imgf1 = loadimg(joinpath(ctx.indir, "frame11.png"))
- maxflow = 4.6157 # Dimetrodon
+ gtflow = FileIO.load(joinpath(ctx.indir, "flow10.flo"))
+ maxflow = OpticalFlowUtils.maxflow(gtflow)
ctx = Util.Context(ctx; imgf0, imgf1, maxflow)
+ mkpath(ctx.outdir)
return optflow(ctx)
end
--
GitLab