diff --git a/src/run.jl b/src/run.jl index ef2083db5cad02cb6be724d5eb9af70a07ee6359..b4c75c9162d874270c5cfe56a1cc6708fea40a16 100644 --- a/src/run.jl +++ b/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