fix l2_pixel projection for non-square images

src/image.jl

@@ -310,6 +310,7 @@ Base.eltype(it::PixelIterator) = eltype(it.grid)
 
 Return an iterator over all integer cartesian indices (corresponding to pixel
 centers) intersecting a given `cell` or geometry `A`.
+The returned indices represent coordinates in the given mesh.
 """
 function PixelIterator(mesh::Mesh, cell::Int)
     A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
@@ -359,7 +360,7 @@ function project_l2_pixel!(u::FeFunction, img)
 
     # first loop to count number of triangles intersecting pixel evaluation
     # points
-    ncells = zeros(axes(img))
+    ncells = zeros(axes(img, 2), axes(img, 1))
     for cell in cells(mesh)
         pixels = PixelIterator(mesh, cell)
         for I in pixels