add inpaint_denoise examples

scripts/run.jl

@@ -27,3 +27,5 @@ ctx(experiment_denoise, "fem/denoise")
 ctx(experiment_inpaint, "fem/inpaint")
 # adaptive optical flow
 ctx(experiment_optflow_middlebury_all, "fem/optflow/middlebury")
+
+# additional experiments for paper

scripts/run_experiments.jl

@@ -28,6 +28,38 @@ isdefined(Main, :Revise) && Revise.track(joinpath(@__DIR__, "util.jl"))
 
 using .Util
 
+module Operators
+
+struct TIdentity end
+(::TIdentity)(tdata, u) = u
+(::TIdentity)(::typeof(adjoint), tdata, v) = v
+Base.show(io, ::MIME"text/latex", ::TIdentity) =
+    print(io, raw"\ensuremath{I}")
+
+struct TInpaint end
+(::TInpaint)(tdata, u) = abs(tdata[begin] - 1.) < 1e-8 ? u : zero(u)
+(T::TInpaint)(::typeof(adjoint), tdata, v) = T(tdata, v)
+Base.show(io, ::MIME"text/latex", ::TInpaint) =
+    print(io, raw"\ensuremath{T_{\mathrm{inpaint}}}")
+
+struct TOptflow end
+(::TOptflow)(tdata, u) = tdata * u
+(::TOptflow)(::typeof(adjoint), tdata, v) = tdata' * v
+Base.show(io, ::MIME"text/latex", ::TOptflow) =
+    print(io, raw"\ensuremath{T_{\mathrm{optflow}}}")
+
+struct SIdentity end
+(::SIdentity)(u, nablau) = u
+Base.show(io, ::MIME"text/latex", ::SIdentity) =
+    print(io, raw"\ensuremath{I}")
+
+struct SNabla end
+(::SNabla)(u, nablau) = nablau'
+Base.show(io, ::MIME"text/latex", ::SNabla) =
+    print(io, raw"\ensuremath{\nabla}")
+
+end
+
 grayclamp(value) = Gray(clamp(value, 0., 1.))
 loadimg(x) = Float64.(FileIO.load(x))
 saveimg(io, x::Array{<:Gray}) = FileIO.save(io, grayclamp.(x))
@@ -156,7 +188,7 @@ function OptFlowState(mesh;
     # tdata will be something like nabla(fw)
     T(tdata, u) = tdata * u
     T(::typeof(adjoint), tdata, v) = tdata' * v
-    S(u, nablau) = nablau
+    S(u, nablau) = nablau'
 
     est = FeFunction(Vest, name = "est")
     g = FeFunction(Vg, name = "g")
@@ -817,12 +849,8 @@ function denoise(ctx)
     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{1}(mesh;
-        T, tdata = nothing, S,
+        T = Operators.TIdentity(), tdata = nothing, S = Operators.SIdentity(),
         ctx.params.alpha1, ctx.params.alpha2,
         ctx.params.lambda, ctx.params.beta,
         ctx.params.gamma1, ctx.params.gamma2)
@@ -1276,33 +1304,37 @@ end
 function inpaint(ctx)
     # expect ctx.params.g_arr
 
+    #project_image! = project_l2_lagrange!
     project_image! = project_qi_lagrange!
-    n_refine = 5
 
     # convert to cartesian coordinates
     g_arr = from_img(ctx.params.g_arr)
     mask_arr = from_img(ctx.params.mask_arr)
 
-    mesh = init_grid(g_arr, floor.(Int, size(g_arr) ./ 2^(n_refine / 2))...)
+    mesh = init_grid(g_arr, floor.(Int, size(g_arr) ./ 2^(ctx.params.n_refine / 2))...)
     mesh_area = area(mesh)
 
     # inpaint specific stuff
     Vg = FeSpace(mesh, P1(), (1,))
     _mask = FeFunction(Vg, name = "mask")
 
-    T(tdata, u) = abs(tdata[begin] - 1.) < 1e-8 ? u : zero(u)
-    T(::typeof(adjoint), tdata, v) = T(tdata, v)
-    S(u, nablau) = u
-
     st = L1L2TVState{1}(mesh;
-        T, tdata = _mask, S,
+        T = Operators.TInpaint(), tdata = _mask, ctx.params.S,
         ctx.params.alpha1, ctx.params.alpha2,
         ctx.params.lambda, ctx.params.beta,
         ctx.params.gamma1, ctx.params.gamma2)
 
     function interpolate_image_data!()
         println("interpolate image data ...")
+        if ctx.params.n_refine == 0
+            # if we use a grid at image resolution, we wan't to avoid
+            # information loss at all cost.
+            # projet_l2_pixel is not suited for inpainting as it is a global
+            # operator, which may leak data from the inpainting domain
+            interpolate!(st.g, x -> evaluate_bilinear(g_arr, x))
+        else
             project_image!(st.g, g_arr)
+        end
         project_image!(st.tdata, mask_arr)
     end
 
@@ -1313,6 +1345,7 @@ function inpaint(ctx)
     pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
 
         interpolate_image_data!()
+
         pvd[0] = save_step(0)
 
         i = 0
@@ -1338,7 +1371,7 @@ function inpaint(ctx)
 
             # refinement stop criterion
             k_refine += 1
-            k_refine > n_refine && break
+            k_refine > ctx.params.n_refine && break
             println("refine ...")
             #estimate_res!(st)
             marked_cells = Set(mark(st; theta = 0.5))
@@ -1362,11 +1395,11 @@ function inpaint(ctx)
     #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, "g.png"), grayclamp.(to_img(g_arr)))
     saveimg(joinpath(ctx.outdir, "output.png"), grayclamp.(to_img(u_sampled)))
     save_csv(joinpath(ctx.outdir, "mesh.csv"), st.u)
     savedata(joinpath(ctx.outdir, "data.tex");
-        ctx.params.eps_newton, n_refine,
+        ctx.params.eps_newton, ctx.params.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
@@ -1377,15 +1410,44 @@ function experiment_inpaint(ctx)
     mask_arr = loadimg(joinpath(ctx.indir, "mask.png"))
     mesh = init_grid(g_arr;)
 
-    df = DataFrame()
-    inpaint(Util.Context(ctx; name = "test", df,
+    inpaint(Util.Context(ctx; name = "test",
+        S = Operators.SIdentity,
+        n_refine = 5,
         g_arr, mask_arr, mesh,
         alpha1 = 0., alpha2 = 50., lambda = 1., beta = 1e-5,
         gamma1 = 1e-4, gamma2 = 1e-4,
-        eps_newton = 1e-4, adaptive = true,
+        eps_newton = 1e-4,
     ))
 end
 
+function experiment_inpaint_denoise(ctx)
+    g_arr = loadimg(joinpath(ctx.indir, "input.png"))
+    mask_arr = loadimg(joinpath(ctx.indir, "mask.png"))
+    mesh = init_grid(g_arr;)
+
+    params = (
+        name = "test",
+        n_refine = 0,
+        g_arr, mask_arr, mesh,
+        noise_sigma = 0.1, noise_p = 0.02,
+        alpha1 = 0.2, alpha2 = 8., lambda = 1., #= beta see below =#
+        gamma1 = 1e-4, gamma2 = 1e-4,
+        eps_newton = 1e-4)
+
+    params.g_arr .= add_noise(g_arr; params.noise_sigma, params.noise_p)
+    saveimg(joinpath(ctx.outdir, "input_noisy.png"), grayclamp.(g_arr))
+
+    ctx(inpaint, "s-identity_beta-large"; params...,
+        S = Operators.SIdentity(), beta = 3e-1)
+    ctx(inpaint, "s-identity_beta-small"; params...,
+        S = Operators.SIdentity(), beta = 3e-4)
+
+    ctx(inpaint, "s-nabla_beta-large"; params...,
+        S = Operators.SNabla(), beta = 5e1)
+    ctx(inpaint, "s-nabla_beta-small"; params...,
+        S = Operators.SNabla(), beta = 5e-2)
+end
+
 function optflow(ctx)
     size(ctx.params.imgf0) == size(ctx.params.imgf1) ||
 	throw(ArgumentError("non-matching image sizes"))