diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index c5739323fdd36bcbf8043c8dbad96267968c1acf..6e830108f4e74e2356990ef7a1e0c40f08915ffe 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