diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 1dbe7ad6d03784bd0e7f315e9ba75fe58145c7a1..292b5fc485323806c0b01624dd85b76ee6284d56 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -6,6 +6,7 @@ import ColorTypes
 import CSV
 import DataFrames: DataFrame
 import FileIO
+using ImageQualityIndexes: assess_psnr, assess_ssim
 using OpticalFlowUtils
 using WriteVTK: paraview_collection
 using Plots
@@ -982,24 +983,45 @@ 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" => interpolate!,
-        "l2_lagrange" => project_l2_lagrange!,
-        #"clement" => projec_clement!,
-        "qi_lagrange" => project_qi_lagrange!,
-        #"qi_lagrange_avg" => project_qi_lagrange!,
-        #"l2_pixel" => project_l2_pixel!,
-    ]
-
-    map(methods) do (name, f)
+    df_psnr = DataFrame()
+    df_ssim = DataFrame()
+
+    for mesh_size in (32, 16, 13)
+        mesh = init_grid(imgf, mesh_size)
+        space = FeSpace(mesh, P1(), (1,))
         u = FeFunction(space)
-        f(u, imgf)
-        save_csv(joinpath(ctx.outdir, name * ".csv"), u)
+
+        # all methods use bilinear interpolation for image evaluations
+        methods = [
+            "nodal" => interpolate!,
+            "l2_lagrange" => project_l2_lagrange!,
+            #"clement" => projec_clement!,
+            "qi_lagrange" => project_qi_lagrange!,
+            #"qi_lagrange_avg" => project_qi_lagrange!,
+            "l2_pixel" => project_l2_pixel!,
+        ]
+
+        qualities = map(methods) do (method, f!)
+            u.data .= false
+
+            f!(u, imgf)
+            save_csv(joinpath(ctx.outdir, "$(mesh_size)_$(method).csv"), u)
+
+            imgu = sample(u)
+            saveimg(joinpath(ctx.outdir, "$(mesh_size)_$(method).png"), imgu)
+
+            return method => (
+                psnr = assess_psnr(imgu, imgf),
+                ssim = assess_ssim(imgu, imgf))
+        end
+        psnr = map(x -> Symbol(first(x)) => last(x).psnr, qualities)
+        ssim = map(x -> Symbol(first(x)) => last(x).ssim, qualities)
+        push!(df_psnr, (;mesh_size, psnr...))
+        push!(df_ssim, (;mesh_size, ssim...))
     end
+
+    CSV.write(joinpath(ctx.outdir, "psnr.csv"), df_psnr)
+    CSV.write(joinpath(ctx.outdir, "ssim.csv"), df_ssim)
+
+    #savedata(joinpath(ctx.outdir, "data.tex"); energy_min, algparams...)
 end