Skip to content
Snippets Groups Projects
Commit 97374898 authored by Stephan Hilb's avatar Stephan Hilb
Browse files

add first experiment

parent c8780f1c
No related branches found
No related tags found
No related merge requests found
......@@ -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"
......
......@@ -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"
......
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
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"))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment