diff --git a/src/function.jl b/src/function.jl
index 1d00318c2ca9e073688e86fefa4de6207c728f07..dfc7cfcd5f8cb98140cabe711f118d90b83d6272 100644
--- a/src/function.jl
+++ b/src/function.jl
@@ -268,9 +268,11 @@ ImageFunction(mesh, img) =
     ImageFunction(mesh, img, Ref(1))
 
 bind!(f::ImageFunction, cell) = f.cell[] = cell
+img_coord(img, x) = (size(img, 1) - x[2] + 1, x[1])
 # TODO: precompute the offset from minimum mesh vertex
 evaluate(f::ImageFunction, xloc) =
-    interpolate_bilinear(f.img, elmap(f.mesh, f.cell[])(xloc) .+ (0.5, 0.5))
+    interpolate_bilinear(f.img,
+        img_coord(f.img, elmap(f.mesh, f.cell[])(xloc) .+ (0.5, 0.5)))
 
 
 struct FacetDivergence{F}
@@ -322,6 +324,10 @@ end
 bind!(f::FacetDivergence, cell) = f.cell[] = cell
 
 
+# evaluate the function on a rectangular grid of coordinates starting at (0.5,
+# 0.5) with integer spacing, i.e. pixel center coordinates of a bounding image.
+# indexing as for images: first index is y downwards, second index is x
+# rightwards
 # TODO: get rid of this hack as soon as we use size=() for scalar functions
 function sample(f::FeFunction)
     out = _sample(f)
@@ -361,7 +367,13 @@ function _sample(f::FeFunction)
             end
 	end
     end
-    return out
+
+    # convert to image indexing
+    d = length(f.space.size)
+    reverse!(out, dims = d + 2)
+    out2 = permutedims(out, (ntuple(identity, d)..., d + 2, d + 1))
+
+    return out2
 end
 
 
diff --git a/src/mesh.jl b/src/mesh.jl
index 39ff191477a5a265fd07c051463afc7127ff3157..a8c5c2e09056dee70516ff9d3836062e6cbc1e33 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -44,10 +44,12 @@ function init_hgrid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
     return HMesh(vertices, cells, zeros(Int, axes(cells)))
 end
 
-init_hgrid(img::Array{<:Any, 2}; type=:vertex) =
+function init_hgrid(img::Array{<:Any, 2}; type=:vertex)
+    s = (size(img, 2), size(img, 1))
     type == :vertex ?
-	init_grid(size(img, 1) - 1, size(img, 2) - 1, (1.0, 1.0), size(img)) :
-	init_grid(size(img, 1), size(img, 2), (0.5, 0.5), size(img) .- (0.5, 0.5))
+        init_hgrid((s .- 1)..., (1.0, 1.0), s) :
+	init_hgrid(s..., (0.5, 0.5), s .- (0.5, 0.5))
+end
 
 function sub_mesh(hmesh::HMesh, subcells = cells(hmesh))
     cells = collect(reshape(reinterpret(Int, hmesh.cells[subcells]), 3, :))
@@ -161,13 +163,17 @@ function init_grid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
     return Mesh(vertices, cells)
 end
 
-init_grid(img::Array{<:Any, 2}; type=:vertex) =
+function init_grid(img::Array{<:Any, 2}; type=:vertex)
+    s = (size(img, 2), size(img, 1))
     type == :vertex ?
-	init_grid(size(img, 1) - 1, size(img, 2) - 1, (1.0, 1.0), size(img)) :
-	init_grid(size(img, 1), size(img, 2), (0.5, 0.5), size(img) .- (0.5, 0.5))
+        init_grid((s .- 1)..., (1.0, 1.0), s) :
+	init_grid(s..., (0.5, 0.5), s .- (0.5, 0.5))
+end
 
-init_grid(img::Array{<:Any, 2}, m::Int, n::Int = m) =
-    init_grid(m, n, (0.5, 0.5), size(img) .- (0.5, 0.5))
+function init_grid(img::Array{<:Any, 2}, m::Int, n::Int = m)
+    s = (size(img, 2), size(img, 1))
+    init_grid(m, n, (0.5, 0.5), s .- (0.5, 0.5))
+end
 
 # horribly implemented, please don't curse me
 function bisect!(mesh::HMesh, marked_cells::Set)
diff --git a/src/operator.jl b/src/operator.jl
index 1a21b05ad99d82580ef8a2fabc49eaf92f31b04c..87c3710161f0f14c0e2786b8932f30a35eb500b9 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -116,10 +116,12 @@ function assemble(space::FeSpace, a, l; params...)
     return A, b
 end
 
+from_img(arr) = permutedims(reverse(arr; dims = 1))
+
 project_img(space::FeSpace, img) =
     (u = FeFunction(space); project_img!(u, img))
 
-# composite midpoint quadrature on lagrange point lattice
+# composite midpoint quadrature on 2d-lagrange point lattice
 function quadrature_composite_lagrange_midpoint(p)
     d_ = 2
     n = binomial(p + 2, 2)