implement optical flow

src/run.jl

@@ -8,15 +8,22 @@ function myrun()
     name = "test"
     mesh = init_grid(50)
     d = 2
-    m = 1
+    m = 2
 
-    Vg = FeSpace(mesh, DP0(), (1,))
+    Vg = FeSpace(mesh, P1(), (1,))
     Vu = FeSpace(mesh, P1(), (m,))
     Vp1 = FeSpace(mesh, DP0(), (1,))
     Vp2 = FeSpace(mesh, DP1(), (m, d))
 
-    g = FeFunction(Vg, name="g")
+    # inpainting
     mask = FeFunction(Vg, name="mask")
+    # optflow
+    f0 = FeFunction(Vg, name="f0")
+    f1 = FeFunction(Vg, name="f1")
+    fw = FeFunction(Vg, name="fw")
+    nablafw = nabla(fw)
+
+    g = FeFunction(Vg, name="g")
     u = FeFunction(Vu, name="u")
     p1 = FeFunction(Vp1, name="p1")
     p2 = FeFunction(Vp2, name="p2")
@@ -36,7 +43,7 @@ function myrun()
 
     alpha1 = 0.
     alpha2 = 10.
-    beta = 0.
+    beta = 1e-3
     lambda = 0.01
     gamma1 = 1e-3
     gamma2 = 1e-3
@@ -44,16 +51,24 @@ function myrun()
     interpolate!(g, x -> norm(x - SA[0.5, 0.5]) < 0.3)
     interpolate!(mask, x -> abs(x[2] - 0.5) > 0.1)
 
+    interpolate!(f0, x -> x[1])
+    interpolate!(f1, x -> x[1] - 0.01)
+    fw.data .= f1.data
+    g_optflow(x; u, f0, fw, nablafw) =
+	nablafw * u - (fw - f0)
+    interpolate!(g, g_optflow; u, f0, fw, nablafw)
+
     #T(data, u) = u
     #tdata = u # FIXME
 
-    @inline T(data, u) = iszero(data) ? zero(u) : u
-    tdata = mask
+    #T(data, u) = iszero(data) ? zero(u) : u
+    #tdata = mask
 
-    #T(data, u) = dot(data, u)
-    #tdata = nabla(fw)
+    T(data, u) = data * u
+    tdata = nabla(fw)
 
-    S(u) = u
+    #S(u, nablau) = u
+    S(u, nablau) = nablau
 
     function dp1_update(x; g, u, p1, du, tdata)
 	m1 = max(gamma1, norm(T(tdata, u) - g))
@@ -88,7 +103,7 @@ function myrun()
 	end
     end
 
-    @inline function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
+    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 :
@@ -104,14 +119,14 @@ function myrun()
 	    dot(cond2, nablaphi)
 
 	aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
-	    beta * dot(S(du), S(phi))
+	    beta * dot(S(du, nabladu), S(phi, nablaphi))
 
 	return a1 + a2 + aB
     end
 
-    @inline function du_l(x, phi, nablaphi; g, u, nablau, tdata)
+    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))
+	    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))