refactor optflow examples

scripts/run_experiments.jl

@@ -1051,6 +1051,7 @@ function optflow(ctx)
     size(ctx.params.imgf0) == size(ctx.params.imgf1) ||
 	throw(ArgumentError("non-matching image sizes"))
 
+    project_image! = project_l2_lagrange!
     tol_newton = 1e-5 # cauchy criterion for inner newton loop
 
     m = 2 # range dim of vector field
@@ -1073,14 +1074,15 @@ function optflow(ctx)
         ctx.params.gamma1, ctx.params.gamma2)
 
     function warp!()
+        println("warp and reproject ...")
         # warp image into imgfw / fw
         imgfw = warp_backwards(imgf1, sample(st.u))
-        project_l2_pixel!(fw, imgfw)
+        project_image!(fw, imgfw)
 
         # recompute optflow operator T based on u0 and fw
-        # TODO: investigate what julia does for singular matrices here
         function tdata_optflow(x; u0_deriv, nablafw)
-            res = nablafw / (I + u0_deriv')
+            #res = nablafw / (I + u0_deriv')
+            res = nablafw
             all(isfinite, res) || throw(DivideError("singular optflow matrix"))
             return res
         end
@@ -1088,15 +1090,17 @@ function optflow(ctx)
             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) =
-            #-(fw - f0)
             st.T(tdata, u0) - (fw - f0)
         interpolate!(st.g, g_optflow; u0 = st.u, f0, fw, st.tdata)
     end
 
     function interpolate_image_data!()
-        project_l2_pixel!(f0, imgf0)
-        project_l2_pixel!(f1, imgf1)
+        println("interpolate image data ...")
+        project_image!(f0, imgf0)
+        project_image!(f1, imgf1)
     end
 
 
@@ -1104,55 +1108,59 @@ function optflow(ctx)
         output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
 	    st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw)
 
-    i = 0
     pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
 
-        pvd[i] = save_step(i)
+        interpolate_image_data!()
+        warp!() # in the first step just to fill st.g
+        pvd[0] = save_step(0)
 
-        #norm_g_old = norm_l2(st.g)
-
-        for j in 1:5
-            println("interpolate image data ...")
-            interpolate_image_data!()
-            println("warp ...")
-            warp!() # in the first step just to fill st.g
-            i += 1
+        i = 0
+        k = 0
+        while true
+            step!(st)
+            k += 1
+            #estimate_pd!(st)
 
-            # interior newton loop
-            norm_step_ = Inf
-            k = 0
-            while norm_step_ > tol_newton && k < 20
-                k += 1
-                println()
-                step!(st)
-                #estimate_pd!(st)
+            norm_step_ = norm_step(st) / sqrt(mesh_area)
+            println("norm_step = $norm_step_")
 
-                norm_step_ = norm_step(st) / sqrt(mesh_area)
-                println("norm_step = $norm_step_")
+            # residual is basically unusable: takes long to compute and
+            # oscillates
+            #norm_residual_ = norm_residual(st) / sqrt(mesh_area)
+            #println("norm_residual = $norm_residual_")
 
-                # residual is basically unusable: takes long to compute and
-                # oscillates
-                #norm_residual_ = norm_residual(st) / sqrt(mesh_area)
-                #println("norm_residual = $norm_residual_")
+            #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
 
-                #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
-            end
+            # interior newton stop criterion
+            norm_step_ > tol_newton && k < 10 && continue
+            k = 0
 
             display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
 
-            estimate_res!(st)
+            i += 1
             pvd[i] = save_step(i)
 
-            println("refine ...")
-            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(mesh, st.d, st.m, st.T, fs.tdata, st.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)
+            # warping stop criterion
+            warp!()
+            if i > 10
+                break
+            end
+
+            # refinement stop criterion
+            if false
+                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(mesh, st.d, st.m, st.T, fs.tdata, st.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)
+                interpolate_image_data!()
+            end
         end
     end
     #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
@@ -1167,9 +1175,11 @@ end
 function experiment_optflow_middlebury(ctx)
     imgf0 = loadimg(joinpath(ctx.indir, "frame10.png"))
     imgf1 = loadimg(joinpath(ctx.indir, "frame11.png"))
-    maxflow = 4.6157 # Dimetrodon
+    gtflow = FileIO.load(joinpath(ctx.indir, "flow10.flo"))
+    maxflow = OpticalFlowUtils.maxflow(gtflow)
 
     ctx = Util.Context(ctx; imgf0, imgf1, maxflow)
+    mkpath(ctx.outdir)
     return optflow(ctx)
 end