From c4b0b0a46ccf16d57e7ed67b4a0c98b552f155e1 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Mon, 28 Feb 2022 19:45:04 +0100
Subject: [PATCH] rewrite adaptive denoise experiment
---
scripts/run_experiments.jl | 125 ++++++++++++++++++++++++-------------
1 file changed, 83 insertions(+), 42 deletions(-)
diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index c573932..6e83010 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -794,62 +794,103 @@ function norm_residual(st::L1L2TVState)
return sqrt(upart2 + ppart2)
end
-function denoise(img; params...)
- m = 1
- img = from_img(img) # coord flip
- #mesh = init_grid(img; type=:vertex)
- mesh = init_grid(img, 5, 5)
+# TODO: deduplicate, cf. optflow()
+function denoise(ctx)
+ # expect ctx.params.g_arr
+
+ project_image! = project_l2_lagrange!
+ eps_newton = 1e-5 # cauchy criterion for inner newton loop
+ n_refine = 6
+
+ # convert to cartesian coordinates
+ g_arr = from_img(ctx.params.g_arr)
+
+ mesh = init_grid(g_arr, floor.(Int, size(g_arr) ./ 2^(n_refine / 2))...)
+ mesh_area = area(mesh)
T(tdata, u) = u
+ T(::typeof(adjoint), tdata, v) = v
S(u, nablau) = u
- st = L1L2TVState{m}(mesh; T, tdata = nothing, S, params...)
+ st = L1L2TVState{1}(mesh;
+ T, tdata = nothing, S,
+ ctx.params.alpha1, ctx.params.alpha2,
+ ctx.params.lambda, ctx.params.beta,
+ ctx.params.gamma1, ctx.params.gamma2)
- project_l2_lagrange!(st.g, img)
- #interpolate!(st.g, x -> evaluate_bilinear(img, x))
- #m = (size(img) .- 1) ./ 2 .+ 1
- #interpolate!(st.g, x -> norm(x .- m) < norm(m .- 1) / 3)
+ function interpolate_image_data!()
+ println("interpolate image data ...")
+ project_image!(st.g, g_arr)
+ end
- save_denoise(st, i) =
- output(st, "output/$(st.name)_$(lpad(i, 5, '0')).vtu",
+ save_step(i) =
+ output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
st.g, st.u, st.p1, st.p2, st.est)
- pvd = paraview_collection("output/$(st.name).pvd")
- pvd[0] = save_denoise(st, 0)
-
- k = 0
- println("primal energy: $(primal_energy(st))")
- while true
- while true
- k += 1
- step!(st)
- estimate_pd!(st)
- pvd[k] = save_denoise(st, k)
- println()
-
- norm_step_ = norm_step(st)
- norm_residual_ = norm_residual(st)
-
- println("ndofs: $(ndofs(st.u.space)), est: $(norm_l2(st.est)))")
- println("primal energy: $(primal_energy(st))")
- println("norm_step: $(norm_step_)")
- println("norm_residual: $(norm_residual_)")
-
- norm_step_ <= 1e-1 && break
- end
- marked_cells = mark(st; theta = 0.5)
- println("refining ...")
- st, _ = refine(st, marked_cells)
- test_mesh(st.mesh)
+ pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
+
+ interpolate_image_data!()
+ pvd[0] = save_step(0)
+
+ i = 0
+ k_newton = 0
+ k_refine = 0
+ while true
+ # interior newton
+ k_newton += 1
+ step!(st)
+ norm_step_ = norm_step(st) / sqrt(mesh_area)
+ println("norm_step = $norm_step_")
- project_l2_lagrange!(st.g, img)
+ # interior newton stop criterion
+ norm_step_ > eps_newton && k_newton < 10 && continue
+ k_newton = 0
- k >= 100 && break
+ # plot
+ i += 1
+ display(plot(grayclamp.(to_img(sample(st.u)))))
+ pvd[i] = save_step(i)
+
+ # refinement stop criterion
+ k_refine += 1
+ k_refine > n_refine && break
+ 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 = L1L2TVState(st; mesh, st.tdata,
+ fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
+ interpolate_image_data!()
+ end
end
- vtk_save(pvd)
+ #CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
+
+ u_sampled = sample(st.u)
+ saveimg(joinpath(ctx.outdir, "g.png"), to_img(g_arr))
+ saveimg(joinpath(ctx.outdir, "output.png"), grayclamp.(to_img(u_sampled)))
+ savedata(joinpath(ctx.outdir, "data.tex");
+ eps_newton, n_refine,
+ st.alpha1, st.alpha2, st.lambda, st.beta, st.gamma1, st.gamma2,
+ width=size(u_sampled, 1), height=size(u_sampled, 2))
return st
end
+function experiment_denoise(ctx)
+ g_arr = loadimg(joinpath(ctx.indir, "input.png"))
+ mesh = init_grid(g_arr;)
+
+ df = DataFrame()
+ denoise(Util.Context(ctx; name = "test", df,
+ g_arr, mesh,
+ alpha1 = 0., alpha2 = 30., lambda = 1., beta = 1e-5,
+ gamma1 = 1e-4, gamma2 = 1e-4,
+ eps_newton = 1e-5, adaptive = true,
+ ))
+end
+
+
function denoise_pd(ctx)
params = ctx.params
m = 1
--
GitLab