diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index f12db0fb14b6bd6081e5dfba137ea7dc374c3491..6b3fdf0a259a69239786bb5c4a652853374eb25f 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -959,3 +959,47 @@ function experiment_optflow_middlebury(ctx)
     ctx = Util.Context(ctx; imgf0, imgf1)
     return optflow(ctx)
 end
+
+
+function test_image(n = 2^6; supersample_factor = 16)
+    q = supersample_factor
+    imgs = zeros(n, n)
+
+    f((a, b)) = (sin(0.5 * 6 / (a^2 + b^2 + 1e-1)) + 1) / 2
+
+    for I in CartesianIndices((n, n))
+        for J in CartesianIndices((q, q))
+            imgs[I] += f((Tuple(I * q - J) .+ 0.5) ./ (n * q))
+        end
+        imgs[I] /= q^2
+    end
+
+    return imgs
+end
+
+
+function experiment_image_mesh_interpolation(ctx)
+    imgf = loadimg(joinpath(ctx.indir, "input.png"))
+
+    # TODO: test other meshes
+    mesh = init_grid(imgf)
+    space = FeSpace(mesh, P1(), (1,))
+    u = FeFunction(space)
+
+    # all methods use bilinear interpolation for image evaluations
+    methods = [
+        #"nodal_interpolation" => interpolate!,
+        ## techniques that use lagrange-lattice quadratures
+        "l2_projection" => project_img!,
+        #"clement" => projec_clement!,
+        "l1_stable_quasi_interpolation" => project_img2!,
+        #"l1_stable_quasi_interpolation_avg" => project_img2!,
+        #"l2_regression" => project_img_regression!,
+    ]
+
+    map(methods) do (name, f)
+        u = FeFunction(space)
+        f(u, imgf)
+        save_csv(joinpath(ctx.outdir, name * ".csv"), u)
+    end
+end