From e40cc4c052ae6fa863c5e00800058374a3bbb745 Mon Sep 17 00:00:00 2001 From: Stephan Hilb <stephan@ecshi.net> Date: Sat, 27 Nov 2021 18:52:00 +0100 Subject: [PATCH] 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 | 16 ++++++++++++++-- src/mesh.jl | 22 ++++++++++++++-------- src/operator.jl | 4 +++- 3 files changed, 31 insertions(+), 11 deletions(-) diff --git a/src/function.jl b/src/function.jl index 1d00318..dfc7cfc 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 39ff191..a8c5c2e 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 1a21b05..87c3710 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) -- GitLab