don't use matrix coordinates within the framework

Julia matrix/image coordinates (i.e. first index downwards, second
rightwards) should not spill into the framework.
Instead the conversion should happen before projecting and after
sampling. Meshes should be created accordingly.

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
 
 

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)

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)