use warping

there is still something wrong apparently since one observes smoothing
after warping?

scripts/run.jl

@@ -7,6 +7,7 @@ import DataFrames: DataFrame
 import FileIO
 using OpticalFlowUtils
 using WriteVTK: paraview_collection
+using Plots
 
 using SemiSmoothNewton
 using SemiSmoothNewton: HMesh, ncells, refine
@@ -691,12 +692,13 @@ function optflow(ctx)
     # coord flip
     imgf0 = from_img(ctx.params.imgf0)
     imgf1 = from_img(ctx.params.imgf1)
+    maxflow = 4.6157 # Dimetrodon
     size(imgf0) == size(imgf1) ||
 	throw(ArgumentError("non-matching dimensions"))
 
     m = 2
     #mesh = init_grid(imgf0; type=:vertex)
-    mesh = init_grid(imgf0, 1, 1)
+    mesh = init_grid(imgf0, (size(imgf0) .÷ 16)...)
     #mesh = init_grid(imgf0)
 
     # optflow specific stuff
@@ -706,25 +708,33 @@ function optflow(ctx)
     fw = FeFunction(Vg, name="fw")
 
     T(tdata, u) = tdata * u # tdata = nablafw
-    #S(u, nablau) = nablau
-    S(u, nablau) = u
+    S(u, nablau) = nablau
+    #S(u, nablau) = u
 
     st = L1L2TVContext("run", mesh, m; T, tdata = nabla(fw), S,
         ctx.params.alpha1, ctx.params.alpha2, ctx.params.lambda, ctx.params.beta,
         ctx.params.gamma1, ctx.params.gamma2)
 
-    function reproject!()
-        project_img!(f0, imgf0)
-        project_img!(f1, imgf1)
+    function warp!()
+        imgfw = warp_backwards(imgf1, sample(st.u))
+        project_img!(fw, imgfw)
 
-        # TODO: should be warping from imgfw
-        fw.data .= f1.data
+        # replace new tdata
+        st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fw), S,
+            st.alpha1, st.alpha2, st.beta, st.lambda, st.gamma1, st.gamma2,
+            st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
 
         g_optflow(x; u, f0, fw, nablafw) =
             nablafw * u - (fw - f0)
-        interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = nabla(fw))
+        interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = st.tdata)
+    end
+
+    function reproject!()
+        project_img2!(f0, imgf0)
+        project_img2!(f1, imgf1)
     end
     reproject!()
+    warp!()
 
     save_step(i) =
         output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
@@ -734,14 +744,26 @@ function optflow(ctx)
     pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
         pvd[i] = save_step(i)
         while true
-            for k in 1:10
+            for j in 1:4
+                norm_g_old = norm_l2(st.g)
+                for k in 1:5
+                    i += 1
+                    step!(st)
+                    estimate!(st)
+                    pvd[i] = save_step(i)
+                    println()
+                end
+                warp!()
+                norm_g = norm_l2(st.g)
                 i += 1
-                step!(st)
-                estimate!(st)
                 pvd[i] = save_step(i)
-                println()
+                display(plot(colorflow(to_img(sample(st.u)); maxflow)))
             end
+            i >= 50 && break
+            #continue
+
             marked_cells = mark(st; theta = 0.5)
+            #marked_cells = Set(axes(mesh.cells, 2))
 
             println("refining ...")
 
@@ -752,12 +774,16 @@ function optflow(ctx)
                 st.alpha1, st.alpha2, st.beta, st.lambda, st.gamma1, st.gamma2,
                 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)
+                i += 1
+                pvd[i] = save_step(i)
 
             println("reprojecting ...")
             reproject!()
+                i += 1
+                pvd[i] = save_step(i)
         end
     end
-    display(plot(colorflow(to_img(sample(st.u)))))
+    display(plot(colorflow(to_img(sample(st.u)); maxflow)))
 
     #CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
     #saveimg(joinpath(ctx.outdir, "output_glob.png"), fetch_u(states.glob))

src/image.jl

-export interpolate_bilinear, halve
+export interpolate_bilinear, halve, warp_backwards
 
 eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...]