add preliminary optflow changes

scripts/run_experiments.jl

@@ -851,17 +851,18 @@ function inpaint(img, imgmask; name, params...)
 end
 
 function optflow(ctx)
-    # coord flip
+    tol_newton = 1e-6 # cauchy criterion
+
+    size(ctx.params.imgf0) == size(ctx.params.imgf1) ||
+	throw(ArgumentError("non-matching image sizes"))
+
+    m = 2 # range dim of vector field
     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, (size(imgf0) .÷ 16)...)
     #mesh = init_grid(imgf0)
+    mesh = init_grid(imgf0, (size(imgf0) .÷ 2^3)...)
+    domain_factor = 1 / sqrt(area(mesh))
 
     # optflow specific stuff
     Vg = FeSpace(mesh, P1(), (1,))
@@ -869,7 +870,7 @@ function optflow(ctx)
     f1 = FeFunction(Vg, name="f1")
     fw = FeFunction(Vg, name="fw")
 
-    T(tdata, u) = tdata * u # tdata = nablafw
+    T(tdata, u) = tdata * u # here tdata will be nabla(fw)
     S(u, nablau) = nablau
     #S(u, nablau) = u
 
@@ -881,22 +882,19 @@ function optflow(ctx)
         imgfw = warp_backwards(imgf1, sample(st.u))
         project_l2_lagrange!(fw, imgfw)
 
-        # 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
             nablafw * u - (fw - f0)
         interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = st.tdata)
     end
 
-    function reproject!()
-        project_qi_lagrange!(f0, imgf0)
-        project_qi_lagrange!(f1, imgf1)
+    function interpolate_image_data!()
+        project_l2_pixel!(f0, imgf0)
+        project_l2_pixel!(f1, imgf1)
     end
-    reproject!()
-    warp!()
+
+    interpolate_image_data!()
+    warp!() # just to fill st.g
 
     save_step(i) =
         output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
@@ -904,31 +902,36 @@ function optflow(ctx)
 
     i = 0
     pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
+
         pvd[i] = save_step(i)
-        while true
-            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)
+
+        #norm_g_old = norm_l2(st.g)
+
+        for j in 1:3
+
+            norm_step_ = Inf
+            while norm_step_ > tol_newton
                 println()
-                end
-                warp!()
-                norm_g = norm_l2(st.g)
                 i += 1
+                step!(st)
+                #estimate!(st)
                 pvd[i] = save_step(i)
-                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))
+                norm_step_ = norm_step(st) * domain_factor
+                println("norm_step = $norm_step_")
 
-            println("refining ...")
+                # residual is basically unusable: takes long to compute and
+                # oscillates
+                #norm_residual_ = norm_residual(st) * domain_factor
+                #println("norm_residual = $norm_residual_")
+
+                display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
+            end
+            # poor man's hash
+            #println(integrate(st.mesh, (x; tdata) -> tdata; st.tdata))
 
+            println("refine ...")
+            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,
                 f0, f1, fw)
@@ -936,16 +939,46 @@ 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
+            #i += 1
             pvd[i] = save_step(i)
+
+            println("interpolate image data ...")
+            interpolate_image_data!()
+
+            println("warp ...")
+            warp!() # yay
         end
+
+
+        #warp!()
+        #norm_g = norm_l2(st.g)
+        #i += 1
+        #pvd[i] = save_step(i)
+        #i >= 50 && break
+        #continue
+
+        #marked_cells = mark(st; theta = 0.5)
+        #marked_cells = Set(axes(mesh.cells, 2))
+
+        #println("refining ...")
+
+        #mesh, fs = refine(mesh, marked_cells;
+        #    st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2,
+        #    f0, f1, fw)
+        #st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fs.fw), S,
+        #    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("interpolate image data ...")
+        #interpolate_image_data!()
+        #    i += 1
+        #    pvd[i] = save_step(i)
     end
-    display(plot(colorflow(to_img(sample(st.u)); maxflow)))
+    display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
 
     #CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
     #saveimg(joinpath(ctx.outdir, "output_glob.png"), fetch_u(states.glob))
@@ -957,8 +990,9 @@ end
 function experiment_optflow_middlebury(ctx)
     imgf0 = loadimg(joinpath(ctx.indir, "frame10.png"))
     imgf1 = loadimg(joinpath(ctx.indir, "frame11.png"))
+    maxflow = 4.6157 # Dimetrodon
 
-    ctx = Util.Context(ctx; imgf0, imgf1)
+    ctx = Util.Context(ctx; imgf0, imgf1, maxflow)
     return optflow(ctx)
 end