diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 6b3fdf0a259a69239786bb5c4a652853374eb25f..a266a0c68e74aeca24f56264208e5269f06face1 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -12,7 +12,7 @@ using Plots
 
 using SemiSmoothNewton
 using SemiSmoothNewton: HMesh, ncells, refine, area, mesh_size, ndofs
-using SemiSmoothNewton: project_img!, project_img2!, project!
+using SemiSmoothNewton: project_l2_lagrange!, project_qi_lagrange!, project!
 using SemiSmoothNewton: vtk_mesh, vtk_append!, vtk_save
 
 include("util.jl")
@@ -501,7 +501,7 @@ function denoise(img; name, params...)
 
     ctx = L1L2TVContext(name, mesh, m; T, tdata = nothing, S, params...)
 
-    project_img!(ctx.g, img)
+    project_l2_lagrange!(ctx.g, img)
     #interpolate!(ctx.g, x -> evaluate_bilinear(img, x))
     #m = (size(img) .- 1) ./ 2 .+ 1
     #interpolate!(ctx.g, x -> norm(x .- m) < norm(m .- 1) / 3)
@@ -538,7 +538,7 @@ function denoise(img; name, params...)
 	ctx, _ = refine(ctx, marked_cells)
 	test_mesh(ctx.mesh)
 
-	project_img!(ctx.g, img)
+	project_l2_lagrange!(ctx.g, img)
 
 	k >= 100 && break
     end
@@ -571,7 +571,7 @@ function denoise_pd(st, img; df=nothing, name, algorithm, params_...)
     sigma = inv(tau * L^2)
     theta = 1.
 
-    #project_img!(st.g, img)
+    #project_l2_lagrange!(st.g, img)
     interpolate!(st.g, x -> evaluate_bilinear(img, x))
     st.u.data .= st.g.data
 
@@ -687,7 +687,7 @@ function denoise_approximation(ctx)
         ctx.params.alpha1, ctx.params.alpha2, ctx.params.lambda, ctx.params.beta,
         ctx.params.gamma1, ctx.params.gamma2)
 
-    #project_img!(st.g, img)
+    #project_l2_lagrange!(st.g, img)
     interpolate!(st.g, x -> evaluate_bilinear(ctx.params.img, x))
     st.u.data .= st.g.data
     #st.u.data .= rand(size(st.u.data))
@@ -877,7 +877,7 @@ function optflow(ctx)
 
     function warp!()
         imgfw = warp_backwards(imgf1, sample(st.u))
-        project_img!(fw, imgfw)
+        project_l2_lagrange!(fw, imgfw)
 
         # replace new tdata
         st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fw), S,
@@ -890,8 +890,8 @@ function optflow(ctx)
     end
 
     function reproject!()
-        project_img2!(f0, imgf0)
-        project_img2!(f1, imgf1)
+        project_qi_lagrange!(f0, imgf0)
+        project_qi_lagrange!(f1, imgf1)
     end
     reproject!()
     warp!()
@@ -988,13 +988,12 @@ function experiment_image_mesh_interpolation(ctx)
 
     # all methods use bilinear interpolation for image evaluations
     methods = [
-        #"nodal_interpolation" => interpolate!,
-        ## techniques that use lagrange-lattice quadratures
-        "l2_projection" => project_img!,
+        "nodal" => interpolate!,
+        "l2_lagrange" => project_l2_lagrange!,
         #"clement" => projec_clement!,
-        "l1_stable_quasi_interpolation" => project_img2!,
-        #"l1_stable_quasi_interpolation_avg" => project_img2!,
-        #"l2_regression" => project_img_regression!,
+        "qi_lagrange" => project_qi_lagrange!,
+        #"qi_lagrange_avg" => project_qi_lagrange!,
+        #"l2_pixel" => project_l2_pixel!,
     ]
 
     map(methods) do (name, f)
diff --git a/src/image.jl b/src/image.jl
index d0d1b332befd15c636890fee41ccdc94332b0be7..1a291fd5c06735ace13a3fafb35f43815568e699 100644
--- a/src/image.jl
+++ b/src/image.jl
@@ -135,10 +135,10 @@ end
 # TODO: refine interface: projection vs interpolation, unify different
 # algorithms
 
-project_img(space::FeSpace, img) =
-    (u = FeFunction(space); project_img!(u, img))
+project_l2_lagrange(space::FeSpace, img) =
+    (u = FeFunction(space); project_l2_lagrange!(u, img))
 
-function project_img!(u::FeFunction, img::AbstractArray)
+function project_l2_lagrange!(u::FeFunction, img::AbstractArray)
     d = 2 # domain dimension
     space = u.space
     mesh = space.mesh
@@ -219,8 +219,8 @@ function project_img!(u::FeFunction, img::AbstractArray)
     return u
 end
 
-project_img2(space::FeSpace, img) =
-    (u = FeFunction(space); project_img2!(u, img))
+project_qi_lagrange(space::FeSpace, img) =
+    (u = FeFunction(space); project_qi_lagrange!(u, img))
 
 # L1-stable quasi-interpolation operator
 # (https://arxiv.org/pdf/1505.06931.pdf)
@@ -228,7 +228,7 @@ project_img2(space::FeSpace, img) =
 # FIXME: currently only approximate quadrature by sampling on lagrange lattice
 # based on element size. Exact evaluation is tricky to implement (lots of
 # intersection handling).
-function project_img2!(u::FeFunction, img)
+function project_qi_lagrange!(u::FeFunction, img)
     u.space.element == P1() ||
         throw(ArgumentError("element unsupported"))
     d = 2 # domain dimension
@@ -296,7 +296,7 @@ function quadrature_cell_pixels(mesh, cell; m0)
 end
 
 
-function interpolate_img_l2pixel!(u::FeFunction, img)
+function project_l2_pixel!(u::FeFunction, img)
     d = 2 # domain dimension
     space = u.space
     mesh = space.mesh