Select Git revision
run_experiments.jl
run_experiments.jl 58.39 KiB
using LinearAlgebra: I, det, dot, norm, normalize
using SparseArrays: sparse, ishermitian
using Statistics: mean, stdm
using Random: MersenneTwister
using Colors: Gray
# avoid world-age-issues by preloading ColorTypes
import ColorTypes
import CSV
import DataFrames: DataFrame
import FileIO
using ForwardDiff: jacobian
using ImageQualityIndexes: assess_psnr, assess_ssim
using OpticalFlowUtils
using WriteVTK: paraview_collection
using Plots
using StaticArrays: MVector, SArray, SMatrix, SVector, SA
using SemiSmoothNewton
using SemiSmoothNewton: HMesh, ncells, refine, area, mesh_size, ndofs, diam,
elmap, assemble_lumped, norm_gradient, integrate_lumped
using SemiSmoothNewton: project!, project_l2_lagrange!, project_qi_lagrange!,
project_l2_pixel!, compute_local_error
using SemiSmoothNewton: vtk_mesh, vtk_append!, vtk_save
include("util.jl")
isdefined(Main, :Revise) && Revise.track(joinpath(@__DIR__, "util.jl"))
using .Util
module Operators
struct TIdentity end
(::TIdentity)(tdata, u) = u
(::TIdentity)(::typeof(adjoint), tdata, v) = v
Base.show(io, ::MIME"text/latex", ::TIdentity) =
print(io, raw"\ensuremath{I}")
struct TInpaint end
(::TInpaint)(tdata, u) = abs(tdata[begin] - 1.) < 1e-8 ? u : zero(u)
(T::TInpaint)(::typeof(adjoint), tdata, v) = T(tdata, v)
Base.show(io, ::MIME"text/latex", ::TInpaint) =
print(io, raw"\ensuremath{T_{\mathrm{inpaint}}}")
struct TOptflow end
(::TOptflow)(tdata, u) = tdata * u
(::TOptflow)(::typeof(adjoint), tdata, v) = tdata' * v
Base.show(io, ::MIME"text/latex", ::TOptflow) =
print(io, raw"\ensuremath{T_{\mathrm{optflow}}}")
struct SIdentity end
(::SIdentity)(u, nablau) = u
Base.show(io, ::MIME"text/latex", ::SIdentity) =
print(io, raw"\ensuremath{I}")
struct SNabla end
(::SNabla)(u, nablau) = nablau'
Base.show(io, ::MIME"text/latex", ::SNabla) =
print(io, raw"\ensuremath{\nabla}")
end
grayclamp(value) = Gray(clamp(value, 0., 1.))
loadimg(x) = Float64.(FileIO.load(x))
saveimg(io, x::Array{<:Gray}) = FileIO.save(io, grayclamp.(x))
saveimg(io, x) = FileIO.save(io, x)
function saveimgdiff(io, f0, f1)
n = 2^8 # colors
cmap = colormap("RdBu", n)
k(v0, v1) = cmap[clamp(ceil(Int, n * ((v1 - v0) / 2 + 0.5)), 1, n)]
saveimg(io, k.(f0, f1))
end
# convert image to/from image coordinate system
from_img(arr) = permutedims(reverse(arr; dims = 1))
to_img(arr) = permutedims(reverse(arr; dims = 2))
function to_img(arr::AbstractArray{<:Any,3})
# for flow fields handle flow direction in first dimension too
out = permutedims(reverse(arr; dims = (1,3)), (1, 3, 2))
out[1, :, :] .= .-out[1, :, :]
return out
end
# Optical Flow measures
function map_flow(f, xs...)
ax = axes(first(xs))[2:end]
return map(CartesianIndices(ax)) do i
vs = (view(x, :, i) for x in xs)
any(v -> any(ismissing, v), vs) && return missing
return f(vs...)
end
end
function angular_error(x, y)
z = (1. + dot(x, y)) / sqrt((1. + dot(x, x)) * (1. + dot(y, y)))
# gracefully handle z > 1 in case of rounding errors
return acos(clamp(z, 0., 1.))
end
function endpoint_error(x, y)
return norm(Iterators.map((x, y) -> x - y, x, y))
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)
struct L1L2TVState{M, m, Ttype, Stype}
mesh::M
T::Ttype
tdata
S::Stype
alpha1::Float64
alpha2::Float64
beta::Float64
lambda::Float64
gamma1::Float64
gamma2::Float64
est::FeFunction
g::FeFunction
u::FeFunction
p1::FeFunction
p2::FeFunction
du::FeFunction
dp1::FeFunction
dp2::FeFunction
end
# for swapping out mesh and discrete functions
# TODO: any of state, mesh or functions should probably be mutable instead
function L1L2TVState(st::L1L2TVState; mesh, tdata = nothing, est, g, u, p1, p2, du, dp1, dp2)
isnothing(tdata) || mesh === tdata.space.mesh ||
throw(ArgumentError("functions have different meshes"))
mesh == est.space.mesh == g.space.mesh ==
u.space.mesh == p1.space.mesh == p2.space.mesh ==
du.space.mesh == dp1.space.mesh == dp2.space.mesh ||
throw(ArgumentError("functions have different meshes"))
return L1L2TVState{typeof(mesh),ndims_u_codomain(st),typeof(st.T),typeof(st.S)}(
mesh, st.T, tdata, st.S,
st.alpha1, st.alpha2, st.beta, st.lambda, st.gamma1, st.gamma2,
est, g, u, p1, p2, du, dp1, dp2)
end
# usual constructor
function L1L2TVState{m}(mesh; T, tdata, S,
alpha1, alpha2, beta, lambda, gamma1, gamma2) where m
alpha2 > 0 || beta > 0 ||
throw(ArgumentError("operator B is singular with these parameters"))
d = ndims_domain(mesh)
Vest = FeSpace(mesh, DP0(), (1,))
Vg = FeSpace(mesh, P1(), (1,))
Vu = FeSpace(mesh, P1(), (m,))
Vp1 = FeSpace(mesh, P1(), (1,))
Vp2 = FeSpace(mesh, DP0(), (d, m))
est = FeFunction(Vest, name = "est")
g = FeFunction(Vg, name = "g")
u = FeFunction(Vu, name = "u")
p1 = FeFunction(Vp1, name = "p1")
p2 = FeFunction(Vp2, name = "p2")
du = FeFunction(Vu; name = "du")
dp1 = FeFunction(Vp1; name = "dp1")
dp2 = FeFunction(Vp2; name = "dp2")
est.data .= 0
g.data .= 0
u.data .= 0
p1.data .= 0
p2.data .= 0
du.data .= 0
dp1.data .= 0
dp2.data .= 0
return L1L2TVState{typeof(mesh),m,typeof(T),typeof(S)}(
mesh, T, tdata, S,
alpha1, alpha2, beta, lambda, gamma1, gamma2,
est, g, u, p1, p2, du, dp1, dp2)
end
function OptFlowState(mesh;
alpha1, alpha2, beta, lambda, gamma1, gamma2, Schoice)
alpha2 > 0 || beta > 0 ||
throw(ArgumentError("operator B is singular with these parameters"))
d = ndims_domain(mesh)
m = 2
Vest = FeSpace(mesh, DP0(), (1,))
# DP1 only for optical flow
Vg = FeSpace(mesh, DP1(), (1,))
Vu = FeSpace(mesh, P1(), (m,))
Vp1 = FeSpace(mesh, P1(), (1,))
Vp2 = FeSpace(mesh, DP0(), (m, d))
Vdg = FeSpace(mesh, DP0(), (1, d))
# tdata will be something like nabla(fw)
T(tdata, u) = tdata * u
T(::typeof(adjoint), tdata, v) = tdata' * v
S(u, nablau) = Schoice == :nabla ? nablau' : u
est = FeFunction(Vest, name = "est")
g = FeFunction(Vg, name = "g")
u = FeFunction(Vu, name = "u")
p1 = FeFunction(Vp1, name = "p1")
p2 = FeFunction(Vp2, name = "p2")
du = FeFunction(Vu; name = "du")
dp1 = FeFunction(Vp1; name = "dp1")
dp2 = FeFunction(Vp2; name = "dp2")
tdata = FeFunction(Vdg, name = "tdata")
est.data .= 0
g.data .= 0
u.data .= 0
p1.data .= 0
p2.data .= 0
du.data .= 0
dp1.data .= 0
dp2.data .= 0
return L1L2TVState{typeof(mesh),m,typeof(T),typeof(S)}(
mesh, T, tdata, S,
alpha1, alpha2, beta, lambda, gamma1, gamma2,
est, g, u, p1, p2, du, dp1, dp2)
end
ndims_u_codomain(st::L1L2TVState{<:Any,m}) where m = m
"calculate l1 norm on reference triangle from langrange dofs"
function p1_ref_l1_expr(ldofs)
u0, u1, u2 = ldofs
return (
(u0 - u1) * abs(u2) * u2 ^ 2 -
((u0 - u2) * abs(u1) * u1 ^ 2 -
(u1 - u2) * abs(u0) * u0 ^ 2)) /
(6 * ((u0 ^ 2 + u1 * u2) * (u1 - u2) - (u1 ^ 2 - u2 ^ 2) * u0))
end
p1_project!(p1, alpha1) = p1_project!(p1, alpha1, p1.space.element)
p1_project!(p1, alpha1, _) = throw(ArgumentError("element unsupported"))
# FIXME: probably only correct for DP0?
# TODO: finish correct discrete projection?!
function p1_project!(p1, alpha1, ::Union{DP0, P1, DP1})
p1.data .= clamp.(p1.data, -alpha1, alpha1)
end
#function p1_project!(p1, alpha1, ::Union{P1, DP1})
# mesh = p1.space.mesh
# for cell in cells(mesh)
# bind!(p1, cell)
# F = elmap(mesh, cell)
#
# ForwardDiff.gradient(p1_ref_l1_expr, p1.ldata)
#
# end
#end
function p2_project!(p2, lambda)
p2.space.element == DP0() || p2.space.element == DP1() ||
p2.space.element == P1() ||
throw(ArgumentError("element unsupported"))
p2d = reshape(p2.data, prod(p2.space.size), :) # no copy
for i in axes(p2d, 2)
p2in = norm(p2d[:, i])
if p2in > lambda
p2d[:, i] .*= lambda ./ p2in
end
end
end
function step!(st::L1L2TVState)
T = st.T
S = st.S
alpha1 = st.alpha1
alpha2 = st.alpha2
beta = st.beta
lambda = st.lambda
gamma1 = st.gamma1
gamma2 = st.gamma2
function du_a(x_, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
cond1 = norm(T(tdata, u) - g) > gamma1 ?
dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
zero(p1)
a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
dot(cond1, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
cond2 = norm(nablau) > gamma2 ?
dot(nablau, nabladu) / norm(nablau)^2 * p2 :
zero(p2)
a2 = lambda / m2 * dot(nabladu, nablaphi) -
dot(cond2, nablaphi')
aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
beta * dot(S(du, nabladu), S(phi, nablaphi))
return a1 + a2 + aB
end
function du_au(x_, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
beta * dot(S(du, nabladu), S(phi, nablaphi))
return aB
end
function du_ap(x_, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
cond1 = norm(T(tdata, u) - g) > gamma1 ?
dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
zero(p1)
a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
dot(cond1, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
cond2 = norm(nablau) > gamma2 ?
dot(nablau, nabladu) / norm(nablau)^2 * p2 :
zero(p2)
a2 = lambda / m2 * dot(nabladu, nablaphi) -
dot(cond2, nablaphi')
return a1 + a2
end
function du_l(x_, phi, nablaphi; g, u, nablau, p1, p2, tdata)
aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
beta * dot(S(u, nablau), S(phi, nablaphi))
m1 = max(gamma1, norm(T(tdata, u) - g))
p1part = alpha1 / m1 * dot(T(tdata, u) - g, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
p2part = lambda / m2 * dot(nablau, nablaphi)
gpart = alpha2 * dot(g, T(tdata, phi))
return -aB - p1part - p2part + gpart
end
function du_lu(x_, phi, nablaphi; g, u, nablau, p1, p2, tdata)
aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
beta * dot(S(u, nablau), S(phi, nablaphi))
gpart = alpha2 * dot(g, T(tdata, phi))
return -aB + gpart
end
function du_lp(x_, phi, nablaphi; g, u, nablau, p1, p2, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
p1part = alpha1 / m1 * dot(T(tdata, u) - g, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
p2part = lambda / m2 * dot(nablau, nablaphi)
return - p1part - p2part
end
# solve du
print("assemble ... ")
Au, bu = assemble(st.du.space, du_au, du_lu;
st.g, st.u, nablau = nabla(st.u), st.p1, st.p2, st.tdata)
Ap, bp = assemble_lumped(st.du.space, du_ap, du_lp;
st.g, st.u, nablau = nabla(st.u), st.p1, st.p2, st.tdata)
print("solve ... ")
A = Au + Ap
b = bu + bp
#println(norm(A - A'))
#println(ishermitian(A))
#A, b = assemble(st.du.space, du_a, du_l;
# st.g, st.u, nablau = nabla(st.u), st.p1, st.p2, st.tdata)
st.du.data .= A \ b
# solve dp1
function dp1_update(x_; g, u, p1, du, tdata)
m1 = max(gamma1, norm(T(tdata, u) - g))
cond = norm(T(tdata, u) - g) > gamma1 ?
dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
zero(p1)
return -p1 + alpha1 / m1 * (T(tdata, u) + T(tdata, du) - g) - cond
end
interpolate!(st.dp1, dp1_update;
st.g, st.u, st.p1, st.du, st.tdata)
# solve dp2
function dp2_update(x_; u, nablau, p2, du, nabladu)
m2 = max(gamma2, norm(nablau))
cond = norm(nablau) > gamma2 ?
dot(nablau, nabladu) / norm(nablau)^2 * p2 :
zero(p2)
return -p2 + lambda / m2 * (nablau + nabladu)' - cond
end
interpolate!(st.dp2, dp2_update;
st.u, nablau = nabla(st.u), st.p2, st.du, nabladu = nabla(st.du))
# newton update
theta = 1.
st.u.data .+= theta * st.du.data
st.p1.data .+= theta * st.dp1.data
st.p2.data .+= theta * st.dp2.data
# reproject p1, p2
# FIXME: the p1 projection destroys the primal reconstruction
p1_project!(st.p1, st.alpha1)
p2_project!(st.p2, st.lambda)
end
"
2010, Chambolle and Pock: primal-dual semi-implicit algorithm
2017, Alkämper and Langer: fem dualisation
"
function step_pd!(st::L1L2TVState; sigma, tau, theta = 1.)
# note: ignores gamma1, gamma2, beta and uses T = I, lambda = 1, m = 1!
# changed alpha2 -> alpha2 / 2
# chambolle-pock require: sigma * tau * L^2 <= 1, L = |grad|
if ndims_u_codomain(st) != 1 || st.lambda != 1. || st.beta != 0.
error("unsupported parameters")
end
beta = tau * st.alpha1 / (1 + tau * st.alpha2)
# u is P1
# p2 is essentially DP0 (technically may be DP1)
# 1. y update
u_new = FeFunction(st.u.space) # x_bar only used here
u_new.data .= st.u.data .+ theta .* st.du.data
p2_next = FeFunction(st.p2.space)
function p2_update(x_; p2, nablau)
return p2 + sigma * nablau
end
interpolate!(p2_next, p2_update; st.p2, nablau = nabla(u_new))
p2_project!(p2_next, st.lambda)
st.dp2.data .= p2_next.data .- st.p2.data
st.p2.data .+= st.dp2.data
# 2.
u_a(x, z, nablaz, phi, nablaphi; g, u, p2) =
dot(z, phi)
u_l(x, phi, nablaphi; u, g, p2) =
(dot(u + tau * st.alpha2 * g, phi) + tau * dot(p2, -nablaphi)) /
(1 + tau * st.alpha2)
# z = 1 / (1 + tau * alpha2) *
# (u + tau * alpha2 * g + tau * div(p))
z = FeFunction(st.u.space)
A, b = assemble(z.space, u_a, u_l; st.g, st.u, st.p2)
z.data .= A \ b
function u_update!(u, z, g, beta)
u.space.element::P1
g.space.element::P1
for i in eachindex(u.data)
if z.data[i] - beta >= g.data[i]
u.data[i] = z.data[i] - beta
elseif z.data[i] + beta <= g.data[i]
u.data[i] = z.data[i] + beta
else
u.data[i] = g.data[i]
end
end
end
u_update!(u_new, z, st.g, beta)
# 3.
# note: step-size control not implemented, since \nabla G^* is not 1/gamma continuous
st.du.data .= u_new.data .- st.u.data
st.u.data .= u_new.data
#st.du.data .= z.data
any(isnan, st.u.data) && error("encountered nan data")
any(isnan, st.p2.data) && error("encountered nan data")
return st
end
"
2010, Chambolle and Pock: accelerated primal-dual semi-implicit algorithm
"
function step_pd2!(st::L1L2TVState; sigma, tau, theta = 1.)
# chambolle-pock require: sigma * tau * L^2 <= 1, L = |grad|
# u is P1
# p2 is essentially DP0 (technically may be DP1)
# 1. y update
u_new = FeFunction(st.u.space) # x_bar only used here
u_new.data .= st.u.data .+ theta .* st.du.data
p1_next = FeFunction(st.p1.space)
p2_next = FeFunction(st.p2.space)
function p1_update(x_; p1, g, u, tdata)
st.alpha1 == 0 && return zero(p1)
return (p1 + sigma * (st.T(tdata, u) - g)) /
(1 + st.gamma1 * sigma / st.alpha1)
end
interpolate!(p1_next, p1_update; st.p1, st.g, st.u, st.tdata)
function p2_update(x_; p2, nablau)
st.lambda == 0 && return zero(p2)
return (p2 + sigma * nablau') /
(1 + st.gamma2 * sigma / st.lambda)
end
interpolate!(p2_next, p2_update; st.p2, nablau = nabla(u_new))
# reproject p1, p2
p1_project!(p1_next, st.alpha1)
p2_project!(p2_next, st.lambda)
st.dp1.data .= p1_next.data .- st.p1.data
st.dp2.data .= p2_next.data .- st.p2.data
st.p1.data .+= st.dp1.data
st.p2.data .+= st.dp2.data
# 2. x update
u_a(x, w, nablaw, phi, nablaphi; g, u, p1, p2, tdata) =
dot(w, phi) +
tau * st.alpha2 * dot(st.T(tdata, w), st.T(tdata, phi)) +
tau * st.beta * dot(st.S(w, nablaw), st.S(phi, nablaphi))
u_au(x, w, nablaw, phi, nablaphi; g, u, p1, p2, tdata) =
dot(w, phi) +
tau * st.alpha2 * dot(st.T(tdata, w), st.T(tdata, phi)) +
tau * st.beta * dot(st.S(w, nablaw), st.S(phi, nablaphi))
u_ap(x, w, nablaw, phi, nablaphi; g, u, p1, p2, tdata) =
0.
u_l(x, phi, nablaphi; g, u, p1, p2, tdata) =
dot(u, phi) - tau * (
dot(p1, st.T(tdata, phi)) +
dot(p2, nablaphi) -
st.alpha2 * dot(g, st.T(tdata, phi)))
u_lu(x, phi, nablaphi; g, u, p1, p2, tdata) =
dot(u, phi) - tau * (-
st.alpha2 * dot(g, st.T(tdata, phi)))
u_lp(x, phi, nablaphi; g, u, p1, p2, tdata) =
- tau * (
dot(p1, st.T(tdata, phi)) +
dot(p2, nablaphi'))
Au, bu = assemble(u_new.space, u_au, u_lu; st.g, st.u, st.p1, st.p2, st.tdata)
Ap, bp = assemble_lumped(u_new.space, u_ap, u_lp; st.g, st.u, st.p1, st.p2, st.tdata)
A = Au + Ap
b = bu + bp
u_new.data .= A \ b
st.du.data .= u_new.data .- st.u.data
st.u.data .= u_new.data
#st.du.data .= z.data
any(isnan, st.u.data) && error("encountered nan data")
any(isnan, st.p1.data) && error("encountered nan data")
any(isnan, st.p2.data) && error("encountered nan data")
return st
end
"
2004, Chambolle: dual semi-implicit algorithm
"
function step_d!(st::L1L2TVState; tau)
# u is P1
# p2 is essentially DP0 (technically may be DP1)
# TODO: this might not be implementable without higher order elements
# need grad(div(p))
return st
end
function solve_primal!(u::FeFunction, st::L1L2TVState)
u_au(x, u, nablau, phi, nablaphi; g, p1, p2, tdata) =
st.alpha2 * dot(st.T(tdata, u), st.T(tdata, phi)) +
st.beta * dot(st.S(u, nablau), st.S(phi, nablaphi))
u_ap(x, u, nablau, phi, nablaphi; g, p1, p2, tdata) =
0.
u_lu(x, phi, nablaphi; g, p1, p2, tdata) =
st.alpha2 * dot(g, st.T(tdata, phi))
u_lp(x, phi, nablaphi; g, p1, p2, tdata) =
-dot(p1, st.T(tdata, phi)) - dot(p2, nablaphi')
# u = -B^{-1} * (T^* p_1 + nabla^* p_2 - alpha2 * T^* g)
# solve:
# <B u, phi> = -<p_1, T phi> - <p_2, nablaphi> + alpha_2 * <g, T phi>
Au, bu = assemble(u.space, u_au, u_lu; st.g, st.p1, st.p2, st.tdata)
Ap, bp = assemble_lumped(u.space, u_ap, u_lp; st.g, st.p1, st.p2, st.tdata)
A = Au + Ap
b = bu + bp
u.data .= A \ b
end
huber(x, gamma) = abs(x) < gamma ? x^2 / (2 * gamma) : abs(x) - gamma / 2
# TODO: finish!
function refine_and_estimate_pd(st::L1L2TVState)
# globally refine
marked_cells = Set(axes(st.mesh.cells, 2))
mesh_new, fs_new = refine(st.mesh, marked_cells;
st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
st_new = L1L2TVState(st; mesh = mesh_new,
fs_new.est, fs_new.g, fs_new.u, fs_new.p1, fs_new.p2,
fs_new.du, fs_new.dp1, fs_new.dp2)
# compute primal dual error indicators
estimate_pd!(st_new)
# transfer data to old state
#
end
# this computes the primal-dual error indicator which is not really useful
# if not computed on a finer mesh than `u` was solved on and has problems with
# alpha1 > 0
function estimate_pd!(st::L1L2TVState)
function estf(x_; g, u, p1, p2, nablau, w, nablaw, tdata)
alpha1part =
st.alpha1 * huber(norm(st.T(tdata, u) - g), st.gamma1) -
dot(st.T(tdata, u) - g, p1) +
(iszero(st.alpha1) ? 0. :
st.gamma1 / (2 * st.alpha1) * norm(p1)^2)
lambdapart =
st.lambda * huber(norm(nablau), st.gamma2) -
dot(nablau', p2) +
(iszero(st.lambda) ? 0. :
st.gamma2 / (2 * st.lambda) * norm(p2)^2)
# avoid non-negative rounding errors
alpha1part = max(0, alpha1part)
lambdapart = max(0, lambdapart)
bpart = 1 / 2 * (
st.alpha2 * dot(st.T(tdata, w - u), st.T(tdata, w - u)) +
st.beta * dot(st.S(w, nablaw) - st.S(u, nablau),
st.S(w, nablaw) - st.S(u, nablau)))
res = alpha1part + lambdapart + bpart
@assert isfinite(res)
return res
end
w = FeFunction(st.u.space)
solve_primal!(w, st)
#w.data .= .-w.data
# TODO: find better name: is actually a cell-wise integration
project!(st.est, estf; st.g, st.u, st.p1, st.p2,
nablau = nabla(st.u), w, nablaw = nabla(w), st.tdata)
st.est.data .= sqrt.(st.est.data)
end
# this computes the residual error indicator which has missing theoretical
# support for alpha1 > 0
function estimate_res!(st::L1L2TVState)
# FIXME: dirty hack
isSnabla = st.S(0, 1) == 1'
norm2(v) = dot(v, v)
cellf(hcell; g, u, nablau, p1, p2, tdata) =
st.alpha2 * st.T(adjoint, tdata, st.T(tdata, u) - g) +
st.T(adjoint, tdata, p1) +
(isSnabla ? zero(u) : st.beta * u)
facetf(hfacet, n; g, u, nablau, p1, p2, tdata) =
# norm2 is calculated later after both facet contributions are
# consolidated
n' * ((isSnabla ? st.beta * st.S(zero(u), nablau) : zero(p2)) + p2)
# manual method
mesh = st.mesh
fs = (;st.g, st.u, nablau = nabla(st.u), st.p1, st.p2, st.tdata)
space = st.est.space
facetv = Dict{Tuple{Int, Int}, MVector{ndims_u_codomain(st), Float64}}()
st.est.data .= 0.
for cell in cells(mesh)
A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
view(mesh.vertices, :, view(mesh.cells, :, cell)))
for f in fs
bind!(f, cell)
end
# cell contribution
hcell = diam(mesh, cell)
cellmap = elmap(mesh, cell)
delmap = jacobian(cellmap, SA[0., 0.])
intel = abs(det(delmap))
qw, qx = SemiSmoothNewton.quadrature()
val = 0.
for k in axes(qx, 2)
xhat = qx[:, k]
x = cellmap(xhat)
opvalues = map(f -> evaluate(f, xhat), fs)
val += qw[k] * intel *
(isSnabla ? hcell^2 : hcell) * norm2(cellf(hcell; opvalues...))
end
gdofs = space.dofmap[:, 1, cell]
st.est.data[gdofs] .= val
# facet contributions
centroid = SArray{Tuple{ndims_domain(mesh)}}(mean(A, dims = 2))
lcentroid = SA[1/3, 1/3]
opvalues = map(f -> evaluate(f, lcentroid), fs)
cross(v) = SA[-v[2], v[1]]
for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3)
fi = mesh.cells[i, cell]
fj = mesh.cells[j, cell]
v = A[:, j] - A[:, i]
hfacet = norm(v)
normal = -normalize(cross(v))
# we don't need to recompute `opvalues` since only piecewise
# constant data is used in `facetf`
facetres = facetf(hfacet, normal; opvalues...)
# average facet contributions
v = SVector{ndims_u_codomain(st)}(facetres / 2)
facetv[(fi, fj)] = get(facetv, (fi, fj), zero(v)) + v
facetv[(fj, fi)] = get(facetv, (fj, fi), zero(v)) + v
end
end
# add facet contributions to cells
for cell in cells(mesh)
for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3)
A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
view(mesh.vertices, :, view(mesh.cells, :, cell)))
v = A[:, j] - A[:, i]
hfacet = norm(v)
intel = hfacet
fi = mesh.cells[i, cell]
fj = mesh.cells[j, cell]
gdofs = space.dofmap[:, 1, cell]
st.est.data[gdofs] .+= (isSnabla ? hfacet : inv(hfacet)) *
norm2(facetv[(fi, fj)]) * intel
end
end
st.est.data .= sqrt.(st.est.data)
end
estimate_error(st::L1L2TVState) =
sqrt(sum(x -> x^2, st.est.data) / area(st.mesh))
# minimal Dörfler marking
function mark(st::L1L2TVState; theta=0.5)
n = ncells(st.mesh)
esttotal = sum(st.est.data)
cellerrors = collect(pairs(st.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 output(st::L1L2TVState, filename, fs...)
println("save \"$filename\" ... ")
vtk = vtk_mesh(filename, st.mesh)
vtk_append!(vtk, fs...)
vtk_save(vtk)
return vtk
end
function primal_energy(st::L1L2TVState)
function integrand(x; g, u, nablau, tdata)
return st.alpha1 * huber(norm(st.T(tdata, u) - g), st.gamma1) +
st.alpha2 / 2 * norm(st.T(tdata, u) - g)^2 +
st.beta / 2 * norm(st.S(u, nablau))^2 +
st.lambda * huber(norm(nablau), st.gamma2)
end
return integrate(st.mesh, integrand; st.g, st.u,
nablau = nabla(st.u), st.tdata)
end
function dual_energy(st::L1L2TVState)
# primal reconstruction
w = FeFunction(st.u.space)
solve_primal!(w, st)
# 1 / 2 * \|w\|_B^2 - alpha2 / 2 * \|g\| + <g, p_1> +
# gamma1 / 2 / alpha1 * \|p_1\|^2 + gamma2 / 2 / lambda * \|p_2\|^2
function integrand(x; g, tdata, w, nablaw, p1, p2)
return 1 / 2 * (
st.alpha2 * dot(st.T(tdata, w), st.T(tdata, w)) +
st.beta * dot(st.S(w, nablaw), st.S(w, nablaw))) +
- st.alpha2 / 2 * dot(g, g) +
dot(g, p1) +
(iszero(st.alpha1) ? 0 : st.gamma1 / 2 / st.alpha1 * dot(p1, p1)) +
(iszero(st.lambda) ? 0 : st.gamma2 / 2 / st.lambda * dot(p2, p2))
end
return integrate(st.mesh, integrand; st.g, st.tdata,
w, nablaw = nabla(w), st.p1, st.p2)
end
norm_l2(f) = norm_l2((x; f) -> f, f.space.mesh; f)
function norm_l2(f::Function, mesh; params...)
f_l2 = function(x; params...)
fx = f(x; params...)
return dot(fx, fx)
end
res = integrate(mesh, f_l2; params...)
return sqrt(res)
end
norm_step(st::L1L2TVState) =
sqrt((norm_l2(st.du)^2 + norm_l2(st.dp1)^2 + norm_l2(st.dp2)^2) / area(st.mesh))
function norm_residual(st::L1L2TVState)
w = FeFunction(st.u.space)
solve_primal!(w, st)
w.data .-= st.u.data
upart2 = norm_l2(w)^2
function integrand(x; g, u, nablau, p1, p2, tdata)
p1part = p1 * max(st.gamma1, norm(st.T(tdata, u) - g)) -
st.alpha1 * (st.T(tdata, u) - g)
p2part = p2 * max(st.gamma2, norm(nablau)) -
st.lambda * nablau'
return norm(p1part)^2 + norm(p2part)^2
end
ppart2 = integrate_lumped(st.mesh, integrand; st.g, st.u,
nablau = nabla(st.u), st.p1, st.p2, st.tdata)
return sqrt(upart2 + ppart2)
end
# 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 = 0
# 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)
st = L1L2TVState{1}(mesh;
T = Operators.TIdentity(), tdata = nothing, S = Operators.SIdentity(),
ctx.params.alpha1, ctx.params.alpha2,
ctx.params.lambda, ctx.params.beta,
ctx.params.gamma1, ctx.params.gamma2)
function interpolate_image_data!()
println("interpolate image data ...")
project_image!(st.g, g_arr)
end
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(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_")
# interior newton stop criterion
norm_step_ > eps_newton && k_newton < 10 && continue
k_newton >= 10 && @warn "Newton reached maximum number of iterations"
k_newton = 0
# plot
i += 1
display(plot(grayclamp.(to_img(sample(st.u)))))
estimate_res!(st)
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
#CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
u_sampled = sample(st.u)
saveimg(joinpath(ctx.outdir, "g.png"), grayclamp.(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,
ctx.params.noise_sigma, ctx.params.noise_p,
width=size(u_sampled, 1), height=size(u_sampled, 2))
return st
end
function add_noise(img; noise_sigma, noise_p)
rng = (seed = 42; MersenneTwister(seed))
#img_noisy = clamp.(img .+ noise_sigma * randn(rng, Float64, size(img)), 0., 1.)
img_noisy = img .+ noise_sigma * randn(rng, Float64, size(img))
for i in eachindex(img)
if rand(rng) < noise_p
img_noisy[i] = rand(rng, Bool) ? 1. : 0.
end
end
return img_noisy
end
function experiment_denoise(ctx)
g_arr = loadimg(joinpath(ctx.indir, "input.png"))
noise_sigma = 0.1
noise_p = 0.02
g_arr = add_noise(g_arr; noise_sigma, noise_p)
saveimg(joinpath(ctx.outdir, "input_noisy.png"), grayclamp.(g_arr))
mesh = init_grid(g_arr;)
df = DataFrame()
denoise(Util.Context(ctx; name = "test", df,
g_arr, mesh,
alpha1 = 0.2, alpha2 = 8., lambda = 1., beta = 0.,
gamma1 = 1e-4, gamma2 = 1e-4,
eps_newton = 1e-5, adaptive = true,
noise_sigma, noise_p
))
end
function denoise_pd(ctx)
t_0 = NaN # initialized later
params = ctx.params
m = 1
mesh = ctx.params.mesh
T(tdata, u) = u
S(u, nablau) = u
st = L1L2TVState{m}(mesh;
T, tdata = nothing, S,
params.alpha1, params.alpha2, params.lambda, params.beta,
params.gamma1, params.gamma2)
# semi-implicit primal dual parameters
mu = st.alpha2 + st.beta # T = I, S = I
#mu /= 100 # kind of arbitrary?
L = norm_gradient(mesh)
#tau = 1e-1
tau = 1 / L
#L = 100
sigma = inv(tau * L^2)
theta = NaN # initialized later
#project_l2_lagrange!(st.g, ctx.params.g_arr)
interpolate!(st.g, x -> evaluate_bilinear(ctx.params.g_arr, x))
if haskey(ctx.params, :init_zero) && ctx.params.init_zero
st.u.data .= 0
else
st.u.data .= st.g.data
end
save_step(st, i) =
output(st, joinpath(ctx.outdir, "output_$(ctx.params.algorithm)_$(lpad(i, 5, '0')).vtu"),
st.g, st.u, st.p1, st.p2)
log!(x::Nothing; kwargs...) = x
function log!(df::DataFrame; kwargs...)
row = NamedTuple(kwargs)
push!(df, row)
#println(df)
println(row)
end
pvd = paraview_collection(joinpath(ctx.outdir, "output_$(ctx.params.algorithm).pvd")) do pvd
k = 0
println("primal energy: $(primal_energy(st))")
while true
k += 1
if ctx.params.algorithm == :pd1
theta = 1 / sqrt(1 + 2 * mu * tau) # constant
# no step size control
step_pd2!(st; sigma, tau, theta)
elseif ctx.params.algorithm == :pd2
theta = 1 / sqrt(1 + 2 * mu * tau)
tau *= theta
sigma /= theta
step_pd2!(st; sigma, tau, theta)
elseif ctx.params.algorithm == :newton
step!(st)
end
#pvd[k] = save_step(st, k)
domain_factor = 1 / sqrt(area(mesh))
norm_step_ = norm_step(st) * domain_factor
#estimate_pd!(st)
t_0 = k == 1 ? time() - 1. : t_0 # avoid 0-start
log!(ctx.params.df; k,
t_elapsed = time() - t_0,
norm_step = norm_step_,
norm_residual = norm_residual(st),
#est = estimate_error(st),
primal_energy = primal_energy(st))
#norm_residual_ < params.eps_newton && norm_step_ < params.eps_newton && break
haskey(params, :eps_newton) && norm_step_ < params.eps_newton && break
haskey(params, :max_iters) && k >= params.max_iters && break
end
pvd[1] = save_step(st, 0)
end
return st
end
function experiment_convergence_rate(ctx)
img = loadimg(joinpath(ctx.indir, "input.png"))
#img = [0. 0.; 1. 0.]
g_arr = from_img(img)
mesh = init_grid(g_arr;)
algparams = (
alpha1=0., alpha2=30., lambda=1., beta=0.,
gamma1=1e-2, gamma2=1e-3,
eps_newton = 1e-10,
max_iters = 10000,
)
algctx = Util.Context(ctx; g_arr, mesh, algparams...)
df1 = DataFrame()
df2 = DataFrame()
df3 = DataFrame()
st1 = denoise_pd(Util.Context(algctx; algorithm = :pd1, df = df1));
st2 = denoise_pd(Util.Context(algctx; algorithm = :pd2, df = df2));
st3 = denoise_pd(Util.Context(algctx; algorithm = :newton, df = df3));
df1[!, :steptime_mean] .= mean(diff(df1.t_elapsed))
df2[!, :steptime_mean] .= mean(diff(df2.t_elapsed))
df3[!, :steptime_mean] .= mean(diff(df3.t_elapsed))
energy_min = min(minimum(df1.primal_energy), minimum(df2.primal_energy),
minimum(df3.primal_energy))
df1[!, :energy_min] .= energy_min
df2[!, :energy_min] .= energy_min
df3[!, :energy_min] .= energy_min
CSV.write(joinpath(ctx.outdir, "semi-implicit.csv"),
logfilter(df1))
CSV.write(joinpath(ctx.outdir, "semi-implicit-accelerated.csv"),
logfilter(df2))
CSV.write(joinpath(ctx.outdir, "newton.csv"),
logfilter(df3))
saveimg(joinpath(ctx.outdir, "input.png"),
img)
saveimg(joinpath(ctx.outdir, "semi-implicit.png"),
to_img(sample(st1.u)))
saveimg(joinpath(ctx.outdir, "semi-implicit-accelerated.png"),
to_img(sample(st2.u)))
saveimg(joinpath(ctx.outdir, "newton.png"),
to_img(sample(st3.u)))
savedata(joinpath(ctx.outdir, "data.tex");
energy_min, algparams...)
end
function experiment_convergence_difference(ctx)
g_arr = from_img([0. 0.; 1. 0.])
mesh = init_grid(g_arr;)
algparams = (
init_zero = true,
alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
gamma1 = 1e-7, gamma2 = 1e-7,
eps_newton = 1e-10,
max_iters = 10000,
)
algctx = Util.Context(ctx; g_arr, mesh, algparams...)
df1 = DataFrame()
df2 = DataFrame()
st1 = denoise_pd(Util.Context(algctx; algorithm = :pd2, df = df1));
st2 = denoise_pd(Util.Context(algctx; algorithm = :newton, df = df2));
p1 = plot(df1.k, df1.norm_step; axis = :log)
p2 = plot(df2.k, df2.norm_step; axis = :log)
p3 = plot(df1.k, df1.norm_residual; axis = :log)
p4 = plot(df2.k, df2.norm_residual; axis = :log)
display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
display(st1.u.data)
display(norm_residual(st1))
display(primal_energy(st1))
display(st2.u.data)
display(norm_residual(st2))
display(primal_energy(st2))
return
end
function experiment_convergence_oscillations(ctx)
g_arr = from_img([0. 0.; 1. 0.])
mesh = init_grid(g_arr;)
algparams = (
#init_zero = true,
#alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
#alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
alpha1 = 0.5, alpha2 = 0., lambda = 1., beta = 1., # oscillating
# DP0 for p1 has less oscillations
gamma1 = 1e-7, gamma2 = 1e-7,
eps_newton = 1e-10,
max_iters = 10000,
)
algctx = Util.Context(ctx; g_arr, mesh, algparams...)
df1 = DataFrame()
df2 = DataFrame()
st1 = denoise_pd(Util.Context(algctx; algorithm = :pd2, df = df1));
st2 = denoise_pd(Util.Context(algctx; algorithm = :newton, df = df2));
p1 = plot(df1.k, df1.norm_step; axis = :log)
p2 = plot(df2.k, df2.norm_step; axis = :log)
p3 = plot(df1.k, df1.primal_energy; xaxis = :log)
p4 = plot(df2.k, df2.primal_energy; xaxis = :log)
display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
display(st1.u.data)
display(primal_energy(st1))
display(st2.u.data)
display(primal_energy(st2))
return
end
function denoise_approximation(ctx)
domain_factor = 1 / sqrt(area(ctx.params.mesh))
m = 1
T(tdata, u) = u
T(::typeof(adjoint), tdata, v) = v
#S(u, nablau) = u
S(u, nablau) = nablau'
st = L1L2TVState{m}(ctx.params.mesh;
T, tdata = nothing, S,
ctx.params.alpha1, ctx.params.alpha2,
ctx.params.lambda, ctx.params.beta,
ctx.params.gamma1, ctx.params.gamma2)
# primal reconstruction
w = FeFunction(st.u.space, name = "w")
w.data .= 0
#project_image! = project_l2_lagrange!
project_image!(f, img) = interpolate!(f, x -> evaluate_bilinear(img, x))
function interpolate_data!()
project_image!(st.g, ctx.params.g_arr)
interpolate!(st.g, x -> float(norm(x .- 1., Inf) <= 3/4))
end
interpolate_data!()
st.u.data .= st.g.data
#st.u.data .= rand(size(st.u.data))
save_vtk(st, i) =
output(st,
joinpath(ctx.outdir, "$(ctx.params.name)_$(lpad(i, 5, '0')).vtu"),
st.g, st.u, st.p1, st.p2, st.est, w)
log!(x::Nothing; kwargs...) = x
function log!(df::DataFrame; kwargs...)
row = NamedTuple(kwargs)
push!(df, row)
println(df)
#println(row)
end
pvd_path = joinpath(ctx.outdir, "$(ctx.params.name).pvd")
pvd = paraview_collection(pvd_path) do pvd
pvd[0] = save_vtk(st, 0)
println("primal energy: $(primal_energy(st))")
# initial refinement
for _ in 1:6
marked_cells = Set(axes(st.mesh.cells, 2))
mesh, fs = refine(st.mesh, marked_cells;
st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2, w)
st = L1L2TVState(st; mesh,
fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
w = fs.w
end
#marked_cells = Set(axes(st.mesh.cells, 2))
#mesh2, fs = refine(st.mesh, marked_cells;
# st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
#st2 = L1L2TVState(st; mesh = mesh2,
# fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
k = 0
k_refine = 0
while true && ndofs(st.u.space) < 2000 && k < 300
k += 1
step!(st)
solve_primal!(w, st)
norm_step_ = norm_step(st) * domain_factor
norm_residual_ = norm_residual(st) * domain_factor
estimate_res!(st)
pvd[k] = save_vtk(st, k)
log!(ctx.params.df; k,
ndofs = ndofs(st.u.space),
hmax = mesh_size(st.mesh),
norm_step = norm_step_, norm_residual = norm_residual_,
primal_energy = primal_energy(st),
dual_energy = dual_energy(st),
est = estimate_error(st),
)
if haskey(ctx.params, :eps_newton) && norm_step_ < ctx.params.eps_newton
k += 1
k_refine += 1
norm_residual_ = norm_residual(st) * domain_factor
estimate_res!(st)
pvd[k] = save_vtk(st, k)
log!(ctx.params.df; k,
ndofs = ndofs(st.u.space),
hmax = mesh_size(st.mesh),
norm_step = norm_step_, norm_residual = norm_residual_,
primal_energy = primal_energy(st),
dual_energy = dual_energy(st),
est = estimate_error(st),
)
marked_cells = ctx.params.adaptive ?
mark(st; theta = 0.5) : # estimator + dörfler
Set(axes(st.mesh.cells, 2))
mesh, fs = refine(st.mesh, marked_cells;
st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2, w)
st = L1L2TVState(st; mesh,
fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
w = fs.w
interpolate_data!()
norm_residual_ = norm_residual(st) * domain_factor
estimate_res!(st)
pvd[k] = save_vtk(st, k)
log!(ctx.params.df; k,
ndofs = ndofs(st.u.space),
hmax = mesh_size(st.mesh),
norm_step = norm_step_, norm_residual = norm_residual_,
primal_energy = primal_energy(st),
dual_energy = dual_energy(st),
est = estimate_error(st),
)
end
end
end
return st
end
function experiment_approximation(ctx)
phi = (sqrt(5) + 1) / 2
g_arr = from_img([0. 0.; 1. 0.])
#g_arr = rand(30, 30)
mesh = init_grid(g_arr;)
df = DataFrame()
denoise_approximation(Util.Context(ctx; name = "test", df,
g_arr, mesh,
alpha1 = 0., alpha2 = 100., lambda = 1., beta = 1e-5,
#alpha1 = 1. - inv(phi), alpha2 = 0., lambda = 0., beta = 1.,
#alpha1 = 2., alpha2 = 0., lambda = 1., beta = 1e-5,
gamma1 = 1e-4, gamma2 = 1e-4,
eps_newton = 1e-5, adaptive = true,
))
end
# TODO: deduplicate, cf. optflow()
# NOTE: used also by inpaint_denoise experiment
function inpaint(ctx)
# expect ctx.params.g_arr
duration = @elapsed begin
#project_image! = project_l2_lagrange!
project_image! = project_qi_lagrange!
# convert to cartesian coordinates
g_arr = from_img(ctx.params.g_arr)
mask_arr = from_img(ctx.params.mask_arr)
dim_g = size(g_arr)
dim_refine = floor.(Int, size(g_arr) ./ 2^(ctx.params.n_refine / 2))
N_refine = floor.(Int, log2( (dim_g[1]-1) * (dim_g[2]-1) / ((dim_refine[1]-1)*(dim_refine[2]-1)) ))
println("Number of refinements: ", N_refine)
mesh = init_grid(g_arr, floor.(Int, size(g_arr) ./ 2^(ctx.params.n_refine / 2))...)
mesh_area = area(mesh)
# inpaint specific stuff
Vg = FeSpace(mesh, P1(), (1,))
_mask = FeFunction(Vg, name = "mask")
st = L1L2TVState{1}(mesh;
T = Operators.TInpaint(), tdata = _mask, ctx.params.S,
ctx.params.alpha1, ctx.params.alpha2,
ctx.params.lambda, ctx.params.beta,
ctx.params.gamma1, ctx.params.gamma2)
function interpolate_image_data!()
println("interpolate image data ...")
if ctx.params.n_refine == 0
# if we use a grid at image resolution, we wan't to avoid
# information loss at all cost.
# project_l2_pixel is not suited for inpainting as it is a global
# operator, which may leak data from the inpainting domain
interpolate!(st.g, x -> evaluate_bilinear(g_arr, x))
interpolate!(st.tdata, x -> evaluate_bilinear(mask_arr, x))
else
project_image!(st.g, g_arr)
# TODO: maybe a more conservative choice would be good
project_image!(st.tdata, mask_arr)
end
end
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(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_")
# interior newton stop criterion
norm_step_ > ctx.params.eps_newton && k_newton < 30 && continue
k_newton >= 30 && @warn "Newton reached maximum number of iterations"
k_newton = 0
# plot
i += 1
display(plot(grayclamp.(to_img(sample(st.u)))))
estimate_res!(st)
pvd[i] = save_step(i)
#break
# refinement stop criterion
k_refine += 1
# k_refine > ctx.params.n_refine && break
k_refine > N_refine && break
println("refine ...")
#estimate_res!(st)
marked_cells = Set(mark(st; theta = 0.5))
## EXPERIMENTAL: mark according to projection error
#g_err = compute_local_error(st.g, g_arr)
## st.est was written out above and used for marking and is no
## longer needed
#st.est.data .= g_err.data
#union!(marked_cells, Set(mark(st; theta = 0.5)))
# manually mark all cell within inpainting domain, since the
# estimator is not reliable there
for cell in cells(mesh)
bind!(st.tdata, cell)
maskv = evaluate(st.tdata, SA[1/3, 1/3])
if abs(maskv[begin] - 1.) > 1e-8
push!(marked_cells, cell)
end
end
#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.tdata)
st = L1L2TVState(st; mesh, fs.tdata,
fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
interpolate_image_data!()
end
end
end # @elapsed
u_sampled = sample(st.u)
saveimg(joinpath(ctx.outdir, "g.png"), grayclamp.(to_img(g_arr)))
saveimg(joinpath(ctx.outdir, "output.png"), grayclamp.(to_img(u_sampled)))
psnr = assess_psnr(u_sampled, g_arr)
ssim = assess_ssim(u_sampled, g_arr)
#save_csv(joinpath(ctx.outdir, "mesh.csv"), st.u)
savedata(joinpath(ctx.outdir, "data.tex");
ctx.params.eps_newton, ctx.params.n_refine,
st.alpha1, st.alpha2, st.lambda, st.beta, st.gamma1, st.gamma2,
width=size(u_sampled, 1), height=size(u_sampled, 2), duration, psnr, ssim, N_refine)
return st
end
function experiment_inpaint_adaptive(ctx)
g_arr = loadimg(joinpath(ctx.indir, "input.png"))
mask_arr = loadimg(joinpath(ctx.indir, "mask.png"))
mesh = init_grid(g_arr;)
for coarsen in 0:15
params = (
name = "test",
n_refine = coarsen,
g_arr, mask_arr, mesh,
#alpha1 = 0.2, alpha2 = 8., lambda = 1.,
alpha1 = 0., alpha2 = 50., lambda = 1.,
beta = 1e-5, S = Operators.SIdentity(),
gamma1 = 1e-4, gamma2 = 1e-4,
eps_newton = 1e-4)
savedata(joinpath(ctx.outdir, "data.tex");
params.eps_newton, params.n_refine,
params.alpha1, params.alpha2, params.lambda,
params.gamma1, params.gamma2,
params.beta, params.S,
width=size(g_arr, 1), height=size(g_arr, 2))
ctx(inpaint, "s-identity_alpha1_$(round(Int, params.alpha1))_alpha2_$(round(Int, params.alpha2))_coarse_$(params.n_refine)"; params...)
end
end
function experiment_inpaint_denoise(ctx)
g_arr = loadimg(joinpath(ctx.indir, "input.png"))
mask_arr = loadimg(joinpath(ctx.indir, "mask.png"))
mesh = init_grid(g_arr;)
params = (
name = "test",
n_refine = 0,
g_arr, mask_arr, mesh,
noise_sigma = 0.1, noise_p = 0.02,
alpha1 = 0.2, alpha2 = 8., lambda = 1., #= beta see below =#
gamma1 = 1e-4, gamma2 = 1e-4,
eps_newton = 1e-4)
params.g_arr .= add_noise(g_arr; params.noise_sigma, params.noise_p)
saveimg(joinpath(ctx.outdir, "input_noisy.png"), grayclamp.(g_arr))
savedata(joinpath(ctx.outdir, "data.tex");
params.eps_newton, params.n_refine,
params.noise_sigma, params.noise_p,
params.alpha1, params.alpha2, params.lambda,
params.gamma1, params.gamma2,
width=size(g_arr, 1), height=size(g_arr, 2))
beta = 3e-1
ctx(inpaint, "s-identity_beta-large"; params...,
S = Operators.SIdentity(), beta = beta * 10)
ctx(inpaint, "s-identity_beta-medium"; params...,
S = Operators.SIdentity(), beta = beta)
ctx(inpaint, "s-identity_beta-small"; params...,
S = Operators.SIdentity(), beta = beta / 10)
beta = 5e+1
ctx(inpaint, "s-nabla_beta-large"; params...,
S = Operators.SNabla(), beta = beta * 10)
ctx(inpaint, "s-nabla_beta-medium"; params...,
S = Operators.SNabla(), beta = beta)
ctx(inpaint, "s-nabla_beta-small"; params...,
S = Operators.SNabla(), beta = beta / 10)
end
function optflow(ctx)
size(ctx.params.imgf0) == size(ctx.params.imgf1) ||
throw(ArgumentError("non-matching image sizes"))
#project_image! = project_l2_lagrange!
project_image!(f, img) =
ctx.params.n_refine == 0 ?
interpolate!(f, x -> evaluate_bilinear(img, x)) : # fine mesh
project_l2_lagrange!(f, img) # coarse adaptive mesh
eps_newton = 1e-3 # cauchy criterion for inner newton loop
eps_warp = 0.05
# convert to cartesian coordinates
imgf0 = from_img(ctx.params.imgf0)
imgf1 = from_img(ctx.params.imgf1)
mesh = init_grid(imgf0, floor.(Int, size(imgf0) ./ 2^(ctx.params.n_refine / 2))...)
mesh_area = area(mesh)
# optflow specific stuff
Vg = FeSpace(mesh, P1(), (1,))
f0 = FeFunction(Vg, name = "f0")
f1 = FeFunction(Vg, name = "f1")
fw = FeFunction(Vg, name = "fw")
# initialized later
f0.data .= NaN
f1.data .= NaN
fw.data .= NaN
st = OptFlowState(mesh;
ctx.params.alpha1, ctx.params.alpha2,
ctx.params.lambda, ctx.params.beta,
ctx.params.gamma1, ctx.params.gamma2, ctx.params.Schoice)
function warp!()
println("warp and reproject ...")
# warp image into imgfw / fw
imgfw = warp_backwards(imgf1, sample(st.u))
project_image!(fw, imgfw)
# recompute optflow operator T based on u0 and fw
function tdata_optflow(x; u0_deriv, nablafw)
#res = nablafw / (I + u0_deriv')
res = nablafw
all(isfinite, res) || error("singular optflow matrix `$res`")
return res
end
interpolate!(st.tdata, tdata_optflow;
u0_deriv = nabla(st.u), nablafw = nabla(fw))
# recompute optflow data g
# note that it is important here that the space of g is general enough
# to avoid any information loss that could screw up image warping
g_optflow(x; u0, f0, fw, tdata) =
st.T(tdata, u0) - (fw - f0)
interpolate!(st.g, g_optflow; u0 = st.u, f0, fw, st.tdata)
end
function interpolate_image_data!()
println("interpolate image data ...")
project_image!(f0, imgf0)
project_image!(f1, imgf1)
end
calc_fdist() = norm_l2((x; f0, fw) -> f0 - fw, mesh; f0, fw)
save_step(i) =
output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw)
pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
interpolate_image_data!()
warp!() # in the first step just to fill st.g
pvd[0] = save_step(0)
i = 0
k_newton = 0
k_warp = 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_")
# interior newton stop criterion
norm_step_ > eps_newton && k_newton < ctx.params.newton_max_iters && continue
k_newton >= ctx.params.newton_max_iters && @warn "Newton reached maximum number of iterations"
ctx.params.warp || break
k_newton = 0
# plot
i += 1
display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
pvd[i] = save_step(i)
# warp
k_warp += 1
fdist0 = calc_fdist()
warp!()
fdist1 = calc_fdist()
rel_datachange = (fdist1 - fdist0) / fdist0
println("rel data change: $(rel_datachange)")
# warping stop criterion
rel_datachange < -eps_warp && continue
ctx.params.refine || break
# refinement stop criterion
k_refine += 1
k_refine > ctx.params.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.tdata, f0, f1, fw)
st = L1L2TVState(st; mesh, fs.tdata,
fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
f0, f1, fw = (fs.f0, fs.f1, fs.fw)
interpolate_image_data!()
end
end
#CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
u_sampled = sample(st.u)
saveimg(joinpath(ctx.outdir, "f0.png"), to_img(imgf0))
saveimg(joinpath(ctx.outdir, "f1.png"), to_img(imgf1))
saveimgdiff(joinpath(ctx.outdir, "g.png"), to_img(imgf0), to_img(imgf1))
u_flow = to_img(u_sampled)
saveimg(joinpath(ctx.outdir, "output.png"), colorflow(u_flow; ctx.params.maxflow))
endpoint_errors = collect(skipmissing(map_flow(endpoint_error, u_flow, ctx.params.gtflow)))
angular_errors = collect(skipmissing(map_flow(angular_error, u_flow, ctx.params.gtflow)))
endpoint_error_mean = mean(endpoint_errors)
endpoint_error_stddev = stdm(endpoint_errors, endpoint_error_mean)
angular_error_mean = mean(angular_errors)
angular_error_stddev = stdm(angular_errors, angular_error_mean)
imgfw = warp_backwards(imgf1, u_sampled)
saveimg(joinpath(ctx.outdir, "fw.png"), to_img(imgfw))
savedata(joinpath(ctx.outdir, "data.tex");
eps_newton, eps_warp, ctx.params.n_refine,
st.alpha1, st.alpha2, st.lambda, st.beta, st.gamma1, st.gamma2, ctx.params.Schoice,
width=size(u_sampled, 2), height=size(u_sampled, 3),
endpoint_error_mean, endpoint_error_stddev,
angular_error_mean, angular_error_stddev,
)
return st
end
function experiment_optflow_middlebury(ctx)
imgf0 = loadimg(joinpath(ctx.indir, "frame10.png"))
imgf1 = loadimg(joinpath(ctx.indir, "frame11.png"))
gtflow = FileIO.load(joinpath(ctx.indir, "flow10.flo"))
maxflow = OpticalFlowUtils._maxflow(gtflow)
ctx = Util.Context(ctx; imgf0, imgf1, maxflow, gtflow)
saveimg(joinpath(ctx.outdir, "ground_truth.png"), colorflow(gtflow; maxflow))
return optflow(ctx)
end
function experiment_optflow_middlebury_all_benchmarks(ctx)
for example in ["Dimetrodon", "Grove2", "Grove3", "Hydrangea",
"RubberWhale", "Urban2", "Urban3", "Venus"]
ctx(experiment_optflow_middlebury, example;
alpha1 = 10., alpha2 = 0., lambda = 1., beta = 1e-5,
gamma1 = 1e-4, gamma2 = 1e-4, Schoice = :nabla, newton_max_iters = 30)
end
end
# FIXME: legacy
function experiment_optflow_middlebury_all(ctx)
experiment_optflow_middlebury_all_benchmarks(
Util.Context(ctx; warp = true, refine = true, n_refine = 6))
end
function experiment_optflow_middlebury_warping_comparison(ctx)
ctx(experiment_optflow_middlebury_all_benchmarks, "vanilla";
warp = false, refine = false, n_refine = 0)
ctx(experiment_optflow_middlebury_all_benchmarks, "warping";
warp = true, refine = false, n_refine = 0)
end
function experiment_optflow_middlebury_warping_comparison_adaptive(ctx)
ctx(experiment_optflow_middlebury_all_benchmarks, "vanilla";
warp = false, refine = false, n_refine = 0)
ctx(experiment_optflow_middlebury_all_benchmarks, "warping";
warp = true, refine = false, n_refine = 0)
ctx(experiment_optflow_middlebury_all_benchmarks, "adaptive-warping";
warp = true, refine = true, n_refine = 6)
end
function experiment_optflow_schoice(ctx)
function run(beta, Schoice)
example = "Dimetrodon"
duration = @elapsed begin
st = ctx(experiment_optflow_middlebury, example;
alpha1 = 10., alpha2 = 0., lambda = 1., beta,
gamma1 = 1e-4, gamma2 = 1e-4, Schoice,
warp = true, refine = false, n_refine = 0, newton_max_iters = 15)
end
u_sampled = sample(st.u)
u_flow = to_img(u_sampled)
# not yet in ctx
gtflow = FileIO.load(joinpath(ctx.indir, example, "flow10.flo"))
endpoint_errors = collect(skipmissing(map_flow(endpoint_error, u_flow, gtflow)))
angular_errors = collect(skipmissing(map_flow(angular_error, u_flow, gtflow)))
endpoint_error_mean = mean(endpoint_errors)
endpoint_error_stddev = stdm(endpoint_errors, endpoint_error_mean)
angular_error_mean = mean(angular_errors)
angular_error_stddev = stdm(angular_errors, angular_error_mean)
return (;Schoice, beta, duration, endpoint_error_mean, endpoint_error_stddev, angular_error_mean, angular_error_stddev)
end
df = DataFrame()
for Schoice in (:id, :nabla)
for beta in 10 .^ range(-6,3,10)
res = run(beta, Schoice)
push!(df, res)
end
end
CSV.write(joinpath(ctx.outdir, "schoice.csv"), df)
end
function test_image(n = 2^6; supersample_factor = 16)
q = supersample_factor
imgs = zeros(n, n)
f((a, b)) = (sin(0.5 * 6 / (a^2 + b^2 + 1e-1)) + 1) / 2
for I in CartesianIndices((n, n))
for J in CartesianIndices((q, q))
imgs[I] += f((Tuple(I * q - J) .+ 0.5) ./ (n * q))
end
imgs[I] /= q^2
end
return imgs
end
function experiment_image_mesh_interpolation(ctx)
imgf = from_img(loadimg(joinpath(ctx.indir, "input.png")))
df_psnr = DataFrame()
df_ssim = DataFrame()
for mesh_size in (32, 16, 13)
mesh = init_grid(imgf, mesh_size)
space = FeSpace(mesh, P1(), (1,))
u = FeFunction(space)
# all methods use bilinear interpolation for image evaluations
methods = [
"nodal" => interpolate!,
"l2_lagrange" => project_l2_lagrange!,
#"clement" => projec_clement!,
"qi_lagrange" => project_qi_lagrange!,
#"qi_lagrange_avg" => project_qi_lagrange!,
"l2_pixel" => project_l2_pixel!,
]
qualities = map(methods) do (method, f!)
u.data .= false
f!(u, imgf)
save_csv(joinpath(ctx.outdir, "$(mesh_size)_$(method).csv"), u)
imgu = sample(u)
clamp01 = x -> clamp(x, 0., 1.)
saveimg(joinpath(ctx.outdir, "$(mesh_size)_$(method).png"), map(clamp01, to_img(imgu)))
return method => (
psnr = assess_psnr(imgu, imgf),
ssim = assess_ssim(imgu, imgf))
end
psnr = map(x -> Symbol(first(x)) => last(x).psnr, qualities)
ssim = map(x -> Symbol(first(x)) => last(x).ssim, qualities)
push!(df_psnr, (;mesh_size, psnr...))
push!(df_ssim, (;mesh_size, ssim...))
end
CSV.write(joinpath(ctx.outdir, "psnr.csv"), df_psnr)
CSV.write(joinpath(ctx.outdir, "ssim.csv"), df_ssim)
#savedata(joinpath(ctx.outdir, "data.tex"); energy_min, algparams...)
end