diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl index 5bbf6eb01a815b30416e65409ff0723a5964b7e0..8a6075f592e8d6ba0ec2b473793a684258123949 100644 --- a/scripts/run_experiments.jl +++ b/scripts/run_experiments.jl @@ -1,4 +1,4 @@ -using LinearAlgebra: norm, dot +using LinearAlgebra: norm, dot, I using Colors: Gray # avoid world-age-issues by preloading ColorTypes @@ -77,7 +77,8 @@ function L1L2TVContext(name, mesh, m; T, tdata, S, d = ndims_domain(mesh) Vest = FeSpace(mesh, DP0(), (1,)) - Vg = FeSpace(mesh, P1(), (1,)) + # DP1 only for optical flow + Vg = FeSpace(mesh, DP1(), (1,)) Vu = FeSpace(mesh, P1(), (m,)) Vp1 = FeSpace(mesh, P1(), (1,)) Vp2 = FeSpace(mesh, DP0(), (m, d)) @@ -376,6 +377,8 @@ end huber(x, gamma) = abs(x) < gamma ? x^2 / (2 * gamma) : abs(x) - gamma / 2 +# this computes the primal-dual error indicator which is not really useful +# if not computed on a finer mesh than `u` was solved on function estimate!(ctx::L1L2TVContext) function estf(x_; g, u, p1, p2, nablau, w, nablaw, tdata) alpha1part = @@ -870,11 +873,13 @@ function optflow(ctx) f0 = FeFunction(Vg, name="f0") f1 = FeFunction(Vg, name="f1") fw = FeFunction(Vg, name="fw") + Vdg = FeSpace(mesh, DP0(), (1, 2)) + tdata = FeFunction(Vdg, name="tdata") T(tdata, u) = tdata * u # here tdata will be nabla(fw) S(u, nablau) = nablau - st = L1L2TVContext("run", mesh, m; T, tdata = nabla(fw), S, + st = L1L2TVContext("run", mesh, m; T, tdata, S, ctx.params.alpha1, ctx.params.alpha2, ctx.params.lambda, ctx.params.beta, ctx.params.gamma1, ctx.params.gamma2) @@ -882,13 +887,25 @@ function optflow(ctx) u_acc.data .= 0 function warp!() + # warp image into imgfw / fw imgfw = warp_backwards(imgf1, sample(u_acc)) project_l2_pixel!(fw, imgfw) - g_optflow(x; u, f0, fw, nablafw) = - -(fw - f0) - #nablafw * u - (fw - f0) - interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = st.tdata) + # 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') + all(isfinite, res) || throw(DivideError("singular optflow matrix")) + return res + end + interpolate!(tdata, tdata_optflow; + u0_deriv = nabla(st.u), nablafw = nabla(fw)) + + # recompute optflow data g + g_optflow(x; u0, f0, fw, tdata) = + #-(fw - f0) + T(tdata, u0) - (fw - f0) + interpolate!(st.g, g_optflow; u0 = st.u, f0, fw, st.tdata) end function interpolate_image_data!() @@ -910,15 +927,17 @@ function optflow(ctx) #norm_g_old = norm_l2(st.g) - for j in 1:3 + for j in 1:5 + i += 1 + # interior newton loop norm_step_ = Inf - while norm_step_ > tol_newton + k = 0 + while norm_step_ > tol_newton && k < 20 + k += 1 println() - i += 1 step!(st) #estimate!(st) - pvd[i] = save_step(i) norm_step_ = norm_step(st) / sqrt(mesh_area) println("norm_step = $norm_step_") @@ -931,13 +950,17 @@ function optflow(ctx) #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow))) end - u_acc.data .+= st.u.data + #estimate!(st) + #pvd[i] = save_step(i) + + 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 = 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, @@ -947,13 +970,19 @@ function optflow(ctx) # 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 #println("interpolate image data ...") #interpolate_image_data!() + + println("warp ...") + ##warp!() # yay + estimate!(st) + pvd[i] = save_step(i) + + #u_acc.data .= st.u.data warp!() # yay end