From 02afebd29d4f810960f1ffd6d6494513742d3446 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Wed, 7 Jul 2021 18:02:54 +0200
Subject: [PATCH] implement inpainting
---
src/run.jl | 66 +++++++++++++++++++++++++++++++-----------------------
1 file changed, 38 insertions(+), 28 deletions(-)
diff --git a/src/run.jl b/src/run.jl
index 51153c7..906012a 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -6,7 +6,7 @@ using LinearAlgebra: norm
function myrun()
name = "test"
- mesh = init_grid(200)
+ mesh = init_grid(50)
d = 2
m = 1
@@ -16,6 +16,7 @@ function myrun()
Vp2 = FeSpace(mesh, DP1(), (m, d))
g = FeFunction(Vg, name="g")
+ mask = FeFunction(Vg, name="mask")
u = FeFunction(Vu, name="u")
p1 = FeFunction(Vp1, name="p1")
p2 = FeFunction(Vp2, name="p2")
@@ -41,18 +42,25 @@ function myrun()
gamma2 = 1e-3
interpolate!(g, x -> norm(x - SA[0.5, 0.5]) < 0.3)
+ interpolate!(mask, x -> abs(x[2] - 0.5) > 0.1)
+
+ #T(data, u) = u
+ #tdata = u # FIXME
+
+ @inline T(data, u) = iszero(data) ? zero(u) : u
+ tdata = mask
+
+ #T(data, u) = dot(data, u)
+ #tdata = nabla(fw)
- # TODO: fixme
- #T(u) = [u[1] + u[2]]
- T(u) = u
S(u) = u
- function dp1_update(x; g, u, p1, du)
- m1 = max(gamma1, norm(T(u) - g))
- cond = norm(T(u) - g) > gamma1 ?
- dot(T(u) - g, T(du)) / norm(T(u) - g)^2 * p1 :
+ 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 :
zeros(size(p1))
- return -p1 + alpha1 / m1 * (T(u) + T(du) - g) - cond
+ return -p1 + alpha1 / m1 * (T(tdata, u) + T(tdata, du) - g) - cond
end
@@ -80,13 +88,13 @@ function myrun()
end
end
- @inline function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2)
- m1 = max(gamma1, norm(T(u) - g))
- cond1 = norm(T(u) - g) > gamma1 ?
- dot(T(u) - g, T(du)) / norm(T(u) - g)^2 * p1 :
+ @inline 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 :
zeros(size(p1))
- a1 = alpha1 / m1 * dot(T(du), T(phi)) -
- dot(cond1, T(phi))
+ a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
+ dot(cond1, T(tdata, phi))
m2 = max(gamma2, norm(nablau))
cond2 = norm(nablau) > gamma2 ?
@@ -95,20 +103,20 @@ function myrun()
a2 = lambda / m2 * dot(nabladu, nablaphi) -
dot(cond2, nablaphi)
- aB = alpha2 * dot(T(du), T(phi)) +
+ aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
beta * dot(S(du), S(phi))
return a1 + a2 + aB
end
- @inline function du_l(x, phi, nablaphi; g, u, nablau)
- aB = alpha2 * dot(T(u), T(phi)) +
+ @inline function du_l(x, phi, nablaphi; g, u, nablau, tdata)
+ aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
beta * dot(S(u), S(phi))
- m1 = max(gamma1, norm(T(u) - g))
- p1part = alpha1 / m1 * dot(T(u) - g, T(phi))
+ 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(phi))
+ gpart = alpha2 * dot(g, T(tdata, phi))
return -aB - p1part - p2part + gpart
end
@@ -122,17 +130,17 @@ function myrun()
pvd = paraview_collection("$(name).pvd")
save_step!(pvd, 0)
- for i = 1:5
- print("newton step $i ...")
+ for i = 1:20
+ print("[$(lpad(i, 3, '0'))] ")
# solve du
- A = assemble(Vu, du_a; g, u, nablau, p1, p2)
- b = assemble_rhs(Vu, du_l; g, u, nablau)
- print(" assembled ...")
+ print("assemble ... ")
+ A = assemble(Vu, du_a; g, u, nablau, p1, p2, tdata)
+ b = assemble_rhs(Vu, du_l; g, u, nablau, tdata)
+ print("solve ... ")
du.data .= A \ b
- println(" solved ...")
# solve dp1, dp2
- interpolate!(dp1, dp1_update; g, u, p1, du)
+ interpolate!(dp1, dp1_update; g, u, p1, du, tdata)
interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu)
# newton update
@@ -144,7 +152,9 @@ function myrun()
p1_project!()
p2_project!()
+ print("save ... ")
save_step!(pvd, i)
+ println()
end
vtk_save(pvd)
return
--
GitLab