scripts/run.jl

 include("run_experiments.jl")
 
-const datapath = joinpath(@__DIR__, "..", "data")
-const outpath = joinpath(@__DIR__, "..", "..", "data")
+const datapath = joinpath(@__DIR__, "../../data")
+const outpath = joinpath(@__DIR__, "../../data")
 
 ctx = Util.Context(datapath, outpath)
 
-# convergence
-
-# algorithm comparison
-#ctx(experiment_convergence_rate, "fem/convergence/rate")
-
-
-# adaptivity
-
-# residual estimator, discrete g, good
-# pd est, L2-TV, discrete g, good?
-# pd est, zero
-# pd est, L1-TV, bad
-
-
-# misc experiments
-
-# inpainting (OLD)
-#ctx(experiment_inpaint, "fem/inpaint")
-# adaptive optical flow
-#ctx(experiment_optflow_middlebury_all, "fem/optflow/middlebury")
+# NOTE: make sure middlebury data is available in each of the directories!
 
 function paper1(ctx)
     # newton vs other algs
@@ -33,12 +14,13 @@ function paper1(ctx)
     ctx(experiment_denoise, "fem/denoise")
     # combined inpainting + denoising with mixed noise
     ctx(experiment_inpaint_denoise, "fem/inpaint_denoise")
-
-    # comparison: warping and adaptivity
+    # comparison: vanilla / warping
     ctx(experiment_optflow_middlebury_warping_comparison, "fem/opticalflow/middlebury_warping_comparison")
 end
 
 function paper2(ctx)
-    # comparison: warping and adaptivity
-    ctx(experiment_optflow_middlebury_warping_comparison_adaptive, "fem/opticalflow/middlebury_warping_comparison")
+    # comparison: image mesh interpolation methods
+    ctx(experiment_image_mesh_interpolation, "image-mesh-interpolation")
+    # comparison: vanilla / warping / adaptive
+    ctx(experiment_optflow_middlebury_warping_comparison_adaptive, "fem/opticalflow/middlebury_warping_comparison_adaptive")
 end

scripts/run_experiments.jl

@@ -1512,7 +1512,8 @@ function optflow(ctx)
     project_image!(f, img) =
         ctx.params.n_refine == 0 ?
             interpolate!(f, x -> evaluate_bilinear(img, x)) : # fine mesh
-            project_l2_lagrange! # coarse adaptive mesh
+            project_l2_lagrange!(f, img) # coarse adaptive mesh
+
 
     eps_newton = 1e-3 # cauchy criterion for inner newton loop
     eps_warp = 0.05
@@ -1529,6 +1530,10 @@ function optflow(ctx)
     f0 = FeFunction(Vg, name = "f0")
     f1 = FeFunction(Vg, name = "f1")
     fw = FeFunction(Vg, name = "fw")
+    # initialized later
+    f0.data .= NaN
+    f1.data .= NaN
+    fw.data .= NaN
 
     st = OptFlowState(mesh;
         ctx.params.alpha1, ctx.params.alpha2,
@@ -1537,6 +1542,7 @@ function optflow(ctx)
 
     function warp!()
         println("warp and reproject ...")
+
         # warp image into imgfw / fw
         imgfw = warp_backwards(imgf1, sample(st.u))
         project_image!(fw, imgfw)
@@ -1671,7 +1677,6 @@ 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, Schoice = :nabla, newton_max_iters = 30)
@@ -1781,7 +1786,8 @@ function experiment_image_mesh_interpolation(ctx)
             save_csv(joinpath(ctx.outdir, "$(mesh_size)_$(method).csv"), u)
 
             imgu = sample(u)
-            saveimg(joinpath(ctx.outdir, "$(mesh_size)_$(method).png"), to_img(imgu))
+            clamp01 = x -> clamp(x, 0., 1.)
+            saveimg(joinpath(ctx.outdir, "$(mesh_size)_$(method).png"), map(clamp01, to_img(imgu)))
 
             return method => (
                 psnr = assess_psnr(imgu, imgf),

src/SemiSmoothNewton.jl

@@ -5,56 +5,4 @@ include("function.jl")
 include("operator.jl")
 include("image.jl")
 
-function clip_line(r0, r1, l0, l1)
-    d_l = -(l1[1] - l0[1])
-    d_r = +(l1[1] - l0[1])
-    d_b = -(l1[2] - l0[2])
-    d_t = +(l1[2] - l0[2])
-
-    t_l = -(r0[1] - l0[1]) / d_l
-    t_r = +(r1[1] - l0[1]) / d_r
-    t_b = -(r0[2] - l0[2]) / d_b
-    t_t = +(r1[2] - l0[2]) / d_t
-
-    t0 = max(
-	d_l <= 0 ? t_l : 0.,
-	d_r <= 0 ? t_r : 0.,
-	d_b <= 0 ? t_b : 0.,
-	d_t <= 0 ? t_t : 0.)
-
-    t1 = min(
-	d_l >= 0 ? t_l : 1.,
-	d_r >= 0 ? t_r : 1.,
-	d_b >= 0 ? t_b : 1.,
-	d_t >= 0 ? t_t : 1.)
-
-    # tolerance on non-intersection
-    t1 < t0 + eps() && return false
-
-    c0 = l0 + t0 * (l1 - l0)
-    c1 = l0 + t1 * (l1 - l0)
-
-    return c0, c1
-end
-
-"project array data onto a grid function"
-function project_array!(u, img::Array)
-    vertices = u.mesh.vertices
-    mesh_box = (
-	reduce((a, b) -> min.(a, b), vertices),
-	reduce((a, b) -> max.(a, b), vertices))
-
-
-    for cell in u.mesh.cells
-	cell_box = (
-	    reduce((a, b) -> min.(a, b), cell.vertices),
-	    reduce((a, b) -> max.(a, b), cell.vertices))
-
-	u[cell]
-
-    end
-end
-
-
-
 end # module

src/image.jl

@@ -15,9 +15,6 @@ function peak_signal_to_noise_ratio(img, img_ref)
     return 10 * log10(imax / mse)
 end
 
-# FIXME: unused?
-from_img(img) = permutedims(reverse(arr; dims = 1))
-
 eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...]
 
 # uses native array indexing, i.e. 1:n

src/mesh.jl

@@ -97,36 +97,6 @@ cells(mesh::Mesh) = axes(mesh.cells, 2)
 vertices(mesh::Mesh) = axes(mesh.vertices, 2)
 vertices(mesh::Mesh, cell) = ntuple(i -> mesh.cells[i, cell], nvertices_cell(mesh))
 
-#function Base.getproperty(obj::Mesh, sym::Symbol)
-#    if sym === :vertices
-#	return reinterpret(SVector{ndims_space(obj), Float64}, vec(obj.vertices))
-#    else
-#	return getfield(obj, sym)
-#    end
-#end
-
-# old grid initialization
-# produces regular grid but not suitable for newest vertex bisection
-function init_grid_old(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
-    r1 = LinRange(v0[1], v1[1], m + 1)
-    r2 = LinRange(v0[2], v1[2], n + 1)
-    coords = collect(Iterators.product(r1, r2))
-
-    vertices = [x[i] for i in 1:2, x in vec(coords)]
-    cells = Array{Int, 2}(undef, 3, 2 * m * n)
-
-    vidx = LinearIndices((m + 1, n + 1))
-    e1 = CartesianIndex(1, 0)
-    e2 = CartesianIndex(0, 1)
-    k = 0
-    for I in CartesianIndices((m, n))
-	cells[:, k += 1] .= (vidx[I], vidx[I + e1], vidx[I + e1 + e2])
-	cells[:, k += 1] .= (vidx[I], vidx[I + e1 + e2], vidx[I + e2])
-    end
-
-    return Mesh(vertices, cells)
-end
-
 # new grid initialization
 # produces regular grid suitable for newest vertex bisection
 function init_grid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
@@ -355,13 +325,6 @@ function geo_contains(A, v)
     return all(λ .>= 0) && sum(λ) <= 1
 end
 
-#function cell_contains(mesh, cell, v)
-#    geo = mesh.vertices[:, mesh.cells[:, cell]]
-#    J = jacobian(x -> geo * [1 - x[1] - x[2], x[1], x[2]], [0., 0.])
-#    λ = J \ (v - geo[:, 1])
-#    return all(λ .>= 0) && sum(λ) <= 1
-#end
-
 function vtk_mesh(filename, mesh::Mesh)
     cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, 3*(i-1)+1:3*(i-1)+3)
 	for i in axes(mesh.cells, 2)]
@@ -424,27 +387,3 @@ function elmap(mesh, cell)
 	@inbounds view(mesh.vertices, :, view(mesh.cells, :, cell)))
     return x -> A * SA[1 - x[1] - x[2], x[1], x[2]]
 end
-
-function test_mesh(mesh)
-    # assemble edge -> cells map
-    edgemap = Dict{NTuple{2, Int}, Vector{Int}}()
-    for cell in cells(mesh)
-	vs = sort(SVector(vertices(mesh, cell)))
-
-	e1 = (vs[1], vs[2])
-	e2 = (vs[1], vs[3])
-	e3 = (vs[2], vs[3])
-
-	edgemap[e1] = push!(get!(edgemap, e1, []), cell)
-	edgemap[e2] = push!(get!(edgemap, e2, []), cell)
-	edgemap[e3] = push!(get!(edgemap, e3, []), cell)
-    end
-    for (edge, cells) in edgemap
-	@assert(1 <= length(cells) <= 2)
-    end
-    # are cells positively oriented?
-    for cell in cells(mesh)
-	delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
-	@assert(det(delmap) > 0)
-    end
-end

test/runtests.jl

@@ -5,6 +5,30 @@ using SemiSmoothNewton
 
 Random.seed!(0)
 
+function test_mesh(mesh)
+    # assemble edge -> cells map
+    edgemap = Dict{NTuple{2, Int}, Vector{Int}}()
+    for cell in cells(mesh)
+	vs = sort(SVector(vertices(mesh, cell)))
+
+	e1 = (vs[1], vs[2])
+	e2 = (vs[1], vs[3])
+	e3 = (vs[2], vs[3])
+
+	edgemap[e1] = push!(get!(edgemap, e1, []), cell)
+	edgemap[e2] = push!(get!(edgemap, e2, []), cell)
+	edgemap[e3] = push!(get!(edgemap, e3, []), cell)
+    end
+    for (edge, cells) in edgemap
+	@assert(1 <= length(cells) <= 2)
+    end
+    # are cells positively oriented?
+    for cell in cells(mesh)
+	delmap = jacobian(elmap(mesh, cell), SA[0., 0.])
+	@assert(det(delmap) > 0)
+    end
+end
+
 @testset "imaging" begin
     img = rand(3, 3)
     mesh = init_grid(img)