compute error measures

scripts/run_experiments.jl

 using LinearAlgebra: I, det, dot, norm, normalize
 using SparseArrays: sparse, ishermitian
-using Statistics: mean
+using Statistics: mean, stdm
 using Random: MersenneTwister
 
 using Colors: Gray
@@ -82,6 +82,26 @@ function to_img(arr::AbstractArray{<:Any,3})
     return out
 end
 
+# Optical Flow measures
+
+function map_flow(f, xs...)
+    ax = axes(first(xs))[2:end]
+    return map(CartesianIndices(ax)) do i
+        vs = (view(x, :, i) for x in xs)
+        any(v -> any(ismissing, v), vs) && return missing
+        return f(vs...)
+    end
+end
+
+function angular_error(x, y)
+    z = (1. + dot(x, y)) / sqrt((1. + dot(x, x)) * (1. + dot(y, y)))
+    # gracefully handle z > 1 in case of rounding errors
+    return acos(clamp(z, 0., 1.))
+end
+
+function endpoint_error(x, y)
+    return norm(Iterators.map((x, y) -> x - y, x, y))
+end
 
 """
     logfilter(dt; a=20)
@@ -881,6 +901,7 @@ function denoise(ctx)
 
             # interior newton stop criterion
             norm_step_ > eps_newton && k_newton < 10 && continue
+            k_newton >= 10 && @warn "Newton reached maximum number of iterations"
             k_newton = 0
 
             # plot
@@ -1360,6 +1381,7 @@ function inpaint(ctx)
 
             # interior newton stop criterion
             norm_step_ > ctx.params.eps_newton && k_newton < 30 && continue
+            k_newton >= 30 && @warn "Newton reached maximum number of iterations"
             k_newton = 0
 
             # plot
@@ -1550,6 +1572,7 @@ function optflow(ctx)
 
             # interior newton stop criterion
             norm_step_ > eps_newton && k_newton < 30 && continue
+            k_newton >= 30 && @warn "Newton reached maximum number of iterations"
             ctx.params.warp || break
             k_newton = 0
 
@@ -1594,13 +1617,25 @@ function optflow(ctx)
     saveimg(joinpath(ctx.outdir, "f0.png"), to_img(imgf0))
     saveimg(joinpath(ctx.outdir, "f1.png"), to_img(imgf1))
     saveimgdiff(joinpath(ctx.outdir, "g.png"), to_img(imgf0), to_img(imgf1))
-    saveimg(joinpath(ctx.outdir, "output.png"), colorflow(to_img(u_sampled); ctx.params.maxflow))
+    u_flow = to_img(u_sampled)
+    saveimg(joinpath(ctx.outdir, "output.png"), colorflow(u_flow; ctx.params.maxflow))
+
+    endpoint_errors = collect(skipmissing(map_flow(endpoint_error, u_flow, ctx.params.gtflow)))
+    angular_errors = collect(skipmissing(map_flow(angular_error, u_flow, ctx.params.gtflow)))
+    endpoint_error_mean = mean(endpoint_errors)
+    endpoint_error_stddev = stdm(endpoint_errors, endpoint_error_mean)
+    angular_error_mean = mean(angular_errors)
+    angular_error_stddev = stdm(angular_errors, angular_error_mean)
+
     imgfw = warp_backwards(imgf1, u_sampled)
     saveimg(joinpath(ctx.outdir, "fw.png"), to_img(imgfw))
     savedata(joinpath(ctx.outdir, "data.tex");
         eps_newton, eps_warp, ctx.params.n_refine,
         st.alpha1, st.alpha2, st.lambda, st.beta, st.gamma1, st.gamma2,
-        width=size(u_sampled, 1), height=size(u_sampled, 2))
+        width=size(u_sampled, 2), height=size(u_sampled, 3),
+        endpoint_error_mean, endpoint_error_stddev,
+        angular_error_mean, angular_error_stddev,
+    )
     return st
 end
 
@@ -1610,7 +1645,7 @@ function experiment_optflow_middlebury(ctx)
     gtflow = FileIO.load(joinpath(ctx.indir, "flow10.flo"))
     maxflow = OpticalFlowUtils._maxflow(gtflow)
 
-    ctx = Util.Context(ctx; imgf0, imgf1, maxflow)
+    ctx = Util.Context(ctx; imgf0, imgf1, maxflow, gtflow)
     saveimg(joinpath(ctx.outdir, "ground_truth.png"), colorflow(gtflow; maxflow))
     return optflow(ctx)
 end
@@ -1618,6 +1653,7 @@ end
 function experiment_optflow_middlebury_all_benchmarks(ctx)
     for example in ["Dimetrodon", "Grove2", "Grove3", "Hydrangea",
             "RubberWhale", "Urban2", "Urban3", "Venus"]
+        #example == "Dimetrodon" && continue
         ctx(experiment_optflow_middlebury, example;
             alpha1 = 10., alpha2 = 0., lambda = 1., beta = 1e-5,
             gamma1 = 1e-4, gamma2 = 1e-4)