implement dörfler marking

src/run.jl

@@ -233,9 +233,26 @@ function refine(ctx::L1L2TVContext, marked_cells)
     return new_ctx
 end
 
+function mark(ctx::L1L2TVContext; theta=0.5)
+    n = ncells(ctx.mesh)
+    esttotal = sum(ctx.est.data)
+
+    cellerrors = collect(pairs(ctx.est.data))
+    cellerrors_sorted = sort(cellerrors; lt = (x, y) -> x.second > y.second)
+
+    marked_cells = Int[]
+    estacc = 0.
+    for (cell, error) in cellerrors_sorted
+	estacc >= theta * esttotal &&  break
+	push!(marked_cells, cell)
+	estacc += error
+    end
+    return marked_cells
+end
+
 
 function save(ctx::L1L2TVContext, filename, fs...)
-    print("save ... ")
+    print("save \"$filename\" ... ")
     vtk = vtk_mesh(filename, ctx.mesh)
     vtk_append!(vtk, fs...)
     vtk_save(vtk)
@@ -318,32 +335,31 @@ function denoise(img; name, params...)
     m = size(img) ./ 2
     interpolate!(ctx.g, x -> norm(x .- m) < norm(m) / 3)
 
-    save_denoise(i) =
+    save_denoise(ctx, i) =
 	save(ctx, "$(ctx.name)_$(lpad(i, 5, '0')).vtu",
 	    ctx.g, ctx.u, ctx.p1, ctx.p2, ctx.est)
 
     pvd = paraview_collection("$(ctx.name).pvd")
-    pvd[0] = save_denoise(0)
-    for i in 1:10
-	step!(ctx)
-	estimate!(ctx)
-	pvd[i] = save_denoise(i)
-	println()
-    end
-    estmean = mean(ctx.est.data)
-    marked_cells = findall(>(2*estmean), ctx.est.data)
-    test_mesh(ctx.mesh)
-    println(marked_cells)
-    ctx = refine(ctx, marked_cells)
-    test_mesh(ctx.mesh)
-    interpolate!(ctx.g, x -> norm(x .- m) < norm(m) / 3)
-    for i in 6:20
-        step!(ctx)
-        estimate!(ctx)
-        pvd[i] = save_denoise(i)
-        println()
+    pvd[0] = save_denoise(ctx, 0)
+
+    k = 0
+    while true
+	for i in 1:10
+	    k += 1
+	    step!(ctx)
+	    estimate!(ctx)
+	    pvd[k] = save_denoise(ctx, k)
+	    println()
+	end
+	marked_cells = mark(ctx; theta = 0.7)
+	#println(marked_cells)
+	println("refining ...")
+	ctx = refine(ctx, marked_cells)
+	test_mesh(ctx.mesh)
+	interpolate!(ctx.g, x -> norm(x .- m) < norm(m) / 3)
+	k >= 50 && break
     end
-
+    vtk_save(pvd)
     return ctx
 end