diff --git a/Manifest.toml b/Manifest.toml index cd18d24b6f12066267c61746e37422a441b3d2cf..af12a4e49fd603a6c067a4a9fb14c7211d48ee73 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -33,6 +33,12 @@ git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" version = "0.11.0" +[[Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.12.8" + [[CommonSubexpressions]] deps = ["MacroTools", "Test"] git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" diff --git a/Project.toml b/Project.toml index d7d50c63b65079138f569595295b6f4bf8dfe5bb..7aef2295931f67d2eae3bcc0ce8da1a25bed06c9 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" diff --git a/src/image.jl b/src/image.jl index ab14b3745eb8bc4de25dac19857e39f4c0c8380f..a8bebd6a4cbd4865e999f546faa469742a3a2157 100644 --- a/src/image.jl +++ b/src/image.jl @@ -1,4 +1,4 @@ -export interpolate_bilinear +export interpolate_bilinear, halve eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...] @@ -28,3 +28,18 @@ function warp_backwards(img, u) end return res end + +function halve(img) + all(iseven.(size(img))) || throw(ArgumentError("uneven size")) + + res = similar(img, size(img) .÷ 2) + fill!(res, zero(eltype(res))) + + for i in CartesianIndices(img) + j = CartesianIndex((Tuple(i) .+ 1) .÷ 2) + res[j] += img[i] + end + res ./= 2 ^ ndims(img) + + return res +end diff --git a/src/run.jl b/src/run.jl index ad33d2b4b12c34650aec09e5338aa61d44d1291b..2478e3e345d91d5bc2dbc57abf35b522862dbdfb 100644 --- a/src/run.jl +++ b/src/run.jl @@ -1,6 +1,7 @@ export myrun, denoise, denoise_pd, inpaint, optflow, solve_primal!, estimate!, loadimg, saveimg using LinearAlgebra: norm +using Colors: Gray # avoid world-age-issues by preloading ColorTypes import ColorTypes @@ -43,7 +44,7 @@ function L1L2TVContext(name, mesh, m; T, tdata, S, Vest = FeSpace(mesh, DP0(), (1,)) Vg = FeSpace(mesh, P1(), (1,)) Vu = FeSpace(mesh, P1(), (m,)) - Vp1 = FeSpace(mesh, P1(), (1,)) + Vp1 = FeSpace(mesh, DP0(), (1,)) Vp2 = FeSpace(mesh, DP1(), (m, d)) est = FeFunction(Vest, name="est") @@ -145,13 +146,15 @@ function step!(ctx::L1L2TVContext) ctx.u, ctx.nablau, ctx.p2, ctx.du, ctx.nabladu) # newton update - ctx.u.data .+= ctx.du.data - ctx.p1.data .+= ctx.dp1.data - ctx.p2.data .+= ctx.dp2.data + theta = 1. + ctx.u.data .+= theta * ctx.du.data + ctx.p1.data .+= theta * ctx.dp1.data + ctx.p2.data .+= theta * ctx.dp2.data # reproject p1 function p1_project!(p1, alpha1) p1.space.element == DP0() || p1.space.element == P1() || + p1.space.element == DP1() || throw(ArgumentError("element unsupported")) p1.data .= clamp.(p1.data, -alpha1, alpha1) end @@ -286,6 +289,7 @@ function step_pd2!(ctx::L1L2TVContext; sigma, tau, theta = 1.) # reproject p1 function p1_project!(p1, alpha1) p1.space.element == DP0() || p1.space.element == P1() || + p1.space.element == DP1() || throw(ArgumentError("element unsupported")) p1.data .= clamp.(p1.data, -alpha1, alpha1) end @@ -365,11 +369,11 @@ function estimate!(ctx::L1L2TVContext) # FIXME: sign? function estf(x_; g, u, p1, p2, nablau, w, nablaw, tdata) alpha1part = iszero(ctx.alpha1) ? 0. : ctx.alpha1 * ( - huber(norm(ctx.T(tdata, u) - g), ctx.gamma1) + + huber(norm(ctx.T(tdata, u) - g), ctx.gamma1) - dot(ctx.T(tdata, u) - g, p1 / ctx.alpha1) + ctx.gamma1 / 2 * norm(p1 / ctx.alpha1)^2) lambdapart = iszero(ctx.lambda) ? 0. : ctx.lambda * ( - huber(norm(nablau), ctx.gamma2) + + huber(norm(nablau), ctx.gamma2) - dot(nablau, p2 / ctx.lambda) + ctx.gamma2 / 2 * norm(p2 / ctx.lambda)^2) bpart = 1 / 2 * ( @@ -448,7 +452,7 @@ end toimg(arr) = Gray.(clamp.(arr, 0., 1.)) loadimg(x) = reverse(transpose(Float64.(FileIO.load(x))); dims = 2) -saveimg(io, x) = FileIO.save(io, toimg(x)) +saveimg(io, x) = FileIO.save(io, transpose(reverse(toimg(x); dims = 2))) function primal_energy(ctx::L1L2TVContext) function integrand(x; g, u, nablau, tdata) @@ -469,7 +473,18 @@ function norm_residual(ctx::L1L2TVContext) w = FeFunction(ctx.u.space) solve_primal!(w, ctx) w.data .-= ctx.u.data - return norm_l2(w) + upart2 = norm_l2(w)^2 + + function integrand(x; g, u, nablau, p1, p2, tdata) + p1part = p1 * max(ctx.gamma1, norm(ctx.T(tdata, u) - g)) - + ctx.alpha1 * (ctx.T(tdata, u) - g) + p2part = p2 * max(ctx.gamma2, norm(nablau)) - + ctx.lambda * nablau + return norm(p1part)^2 + norm(p2part)^2 + end + ppart2 = integrate(ctx.mesh, integrand; ctx.g, ctx.u, ctx.nablau, ctx.p1, ctx.p2, ctx.tdata) + + return sqrt(upart2 + ppart2) end function denoise(img; name, params...) @@ -527,7 +542,8 @@ function denoise(img; name, params...) return ctx end -function denoise_pd(img; name, params...) + +function denoise_pd(img; df=nothing, name, algorithm, params...) m = 1 mesh = init_grid(img; type=:vertex) #mesh = init_grid(img, 5, 5) @@ -539,19 +555,30 @@ function denoise_pd(img; name, params...) # semi-implicit primal dual parameters gamma = ctx.alpha2 + ctx.beta # T = I, S = I - sigma = 1e-2 - tau = 1e-2 + gamma /= 100 # kind of arbitrary? + + tau = 1e-1 + L = 100 + sigma = inv(tau * L^2) theta = 1. - project_img!(ctx.g, img) - #interpolate!(ctx.g, x -> interpolate_bilinear(img, x)) - #m = (size(img) .- 1) ./ 2 .+ 1 - #interpolate!(ctx.g, x -> norm(x .- m) < norm(m .- 1) / 3) + #project_img!(ctx.g, img) + interpolate!(ctx.g, x -> interpolate_bilinear(img, x)) ctx.u.data .= ctx.g.data save_denoise(ctx, i) = output(ctx, "output/$(ctx.name)_$(lpad(i, 5, '0')).vtu", - ctx.g, ctx.u, ctx.p1, ctx.p2, ctx.est, ctx.du) + ctx.g, ctx.u, ctx.p1, ctx.p2) + + log!(x::Nothing; kwargs...) = x + function log!(df::DataFrame; k, norm_step, norm_residual) + push!(df, (; + k, + primal_energy = primal_energy(ctx), + norm_step, + norm_residual)) + println(NamedTuple(last(df))) + end pvd = paraview_collection("output/$(ctx.name).pvd") pvd[0] = save_denoise(ctx, 0) @@ -559,44 +586,70 @@ function denoise_pd(img; name, params...) k = 0 println("primal energy: $(primal_energy(ctx))") - algorithm = :pd2 while true k += 1 - #println("") - if algorithm == :pd2 #|| primal_energy(ctx) < 2409 - println("step_pd: theta = $theta, tau = $tau, sigma = $sigma") - #theta = 1 / sqrt(1 + 2 * gamma * tau) - #tau *= theta - #sigma /= theta - step_pd2!(ctx; sigma, tau, theta) - elseif algorithm == :pd1 + if algorithm == :pd1 # no step size control - step_pd!(ctx; sigma, tau, theta) + step_pd2!(ctx; sigma, tau, theta) + elseif algorithm == :pd2 + theta = 1 / sqrt(1 + 2 * gamma * tau) + tau *= theta + sigma /= theta + step_pd2!(ctx; sigma, tau, theta) elseif algorithm == :newton step!(ctx) end - #estimate!(ctx) #pvd[k] = save_denoise(ctx, k) - println() domain_factor = 1 / sqrt(area(mesh)) norm_step_ = norm_step(ctx) * domain_factor - residual = -1. - residual = norm_residual(ctx) * domain_factor + norm_residual_ = norm_residual(ctx) * domain_factor - #println("ndofs: $(ndofs(ctx.u.space)), est: $(norm_l2(ctx.est)))") - println("primal energy: $(primal_energy(ctx))") - println("norm_step: $(norm_step_)") - println("norm_residual: $(residual)") + log!(df; k, norm_step = norm_step_, norm_residual = norm_residual_) - #return ctx - #residual <= 1e-3 && break - #k >= 5 && break + norm_residual_ < 1e-6 && norm_step_ < 1e-6 && break + norm_step_ < 1e-13 && break end + pvd[1] = save_denoise(ctx, 1) vtk_save(pvd) return ctx end +""" + logfilter(dt; a=20) + +Reduce dataset such that the resulting dataset would have approximately +equidistant datapoints on a log-scale. `a` controls density. +""" +logfilter(dt; a=20) = filter(:k => k -> round(a*log(k+1)) - round(a*log(k)) > 0, dt) + +function experiment1(img) + path = "data/pd-comparison" + + img = loadimg("data/denoising/input.png") + + algparams = (alpha1=0., alpha2=30., lambda=1., beta=0., gamma1=1e-3, gamma2=1e-3) + + df1 = DataFrame() + df2 = DataFrame() + df3 = DataFrame() + + denoise_pd(img; name="test", algorithm=:pd1, df = df1, algparams...); + denoise_pd(img; name="test", algorithm=:pd2, df = df2, algparams...); + denoise_pd(img; name="test", algorithm=:newton, df = df3, algparams...); + + energy_min = min(minimum(df1.primal_energy), minimum(df2.primal_energy), + minimum(df3.primal_energy)) + + df1.primal_energy .-= energy_min + df2.primal_energy .-= energy_min + df3.primal_energy .-= energy_min + + CSV.write(joinpath(path, "semiimplicit.csv"), logfilter(df1)) + CSV.write(joinpath(path, "semiimplicit-accelerated.csv"), logfilter(df2)) + CSV.write(joinpath(path, "newton.csv"), logfilter(df3)) +end + function inpaint(img, imgmask; name, params...) size(img) == size(imgmask) || throw(ArgumentError("non-matching dimensions"))