From 83c0116d82510a77ab7cc13373e43ae9ec251f43 Mon Sep 17 00:00:00 2001 From: Stephan Hilb <stephan@ecshi.net> Date: Thu, 16 Dec 2021 10:47:14 +0100 Subject: [PATCH] refactor --- src/image.jl | 273 +++++++++++++++++++++++++++------------------------ 1 file changed, 146 insertions(+), 127 deletions(-) diff --git a/src/image.jl b/src/image.jl index e9b315e..609c9a7 100644 --- a/src/image.jl +++ b/src/image.jl @@ -1,23 +1,19 @@ export evaluate_bilinear, halve, warp_backwards -# ImageFunction +# General Utility Methods on Images -# TODO: inherit from some abstract function type -struct ImageFunction{Img} - mesh::Mesh - img::Img - cell::Base.RefValue{Int} +function mean_squared_error(img, img_ref) + axes(img) == axes(img_ref) || + throw(ArgumentError("images need matching axes")) + n = length(img) + return sum((img .- img_ref) .^ 2 ./ n) end -ImageFunction(mesh, img) = - ImageFunction(mesh, img, Ref(1)) - -bind!(f::ImageFunction, cell) = f.cell[] = cell -# transform coordinates to image/matrix indexing space -img_coord(img, x) = (size(img, 1) - x[2] + 1, x[1]) -evaluate(f::ImageFunction, xloc) = - evaluate_bilinear(f.img, - img_coord(f.img, elmap(f.mesh, f.cell[])(xloc))) +function peak_signal_to_noise_ratio(img, img_ref) + mse = mean_squared_error(img, img_ref) + imax = 1.0 # maximum intensity + return 10 * log10(imax / mse) +end # FIXME: unused? from_img(img) = permutedims(reverse(arr; dims = 1)) @@ -31,10 +27,10 @@ function evaluate_bilinear(img, x) val = zero(eltype(img)) for idx in CartesianIndices(ntuple(_->0:1, ndims(img))) - cornerbool = Bool.(Tuple(idx)) - λ = ifelse.(cornerbool, x .- x0, x1 .- x) - corner = ifelse.(cornerbool, x1, x0) - val += prod(λ) * eval_neumann(img, corner) + cornerbool = Bool.(Tuple(idx)) + λ = ifelse.(cornerbool, x .- x0, x1 .- x) + corner = ifelse.(cornerbool, x1, x0) + val += prod(λ) * eval_neumann(img, corner) end return val end @@ -42,12 +38,12 @@ end function warp_backwards(img, u) d = ndims(img) axes(img) == axes(u)[2:end] && 1:d == axes(u, 1) || - throw(ArgumentError("invalid dimensions")) + throw(ArgumentError("invalid dimensions")) res = similar(img) for I in CartesianIndices(img) - x = Tuple(I) .+ ntuple(i -> u[i, I], d) - res[I] = evaluate_bilinear(img, x) + x = Tuple(I) .+ ntuple(i -> u[i, I], d) + res[I] = evaluate_bilinear(img, x) end return res end @@ -67,6 +63,24 @@ function halve(img) return res end +# ImageFunction mesh function wrapper + +# TODO: inherit from some abstract mesh function type +# TODO: make mesh type a parameter for performance +struct ImageFunction{Img} + mesh::Mesh + img::Img + cell::Base.RefValue{Int} +end + +ImageFunction(mesh, img) = ImageFunction(mesh, img, Ref(0)) + +bind!(f::ImageFunction, cell) = f.cell[] = cell +# transform coordinates to image/matrix indexing space +img_coord(img, x) = (size(img, 1) - x[2] + 1, x[1]) +evaluate(f::ImageFunction, xloc) = evaluate_bilinear(f.img, + img_coord(f.img, elmap(f.mesh, f.cell[])(xloc))) + # Sampling # evaluate the function on a rectangular grid of coordinates with integer @@ -93,15 +107,15 @@ function _sample(f::FeFunction) out = similar(f.data, outsize) for cell in cells(mesh) bind!(f, cell) - A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}( - view(mesh.vertices, :, view(mesh.cells, :, cell))) + A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}( + view(mesh.vertices, :, view(mesh.cells, :, cell))) # TODO: cleanup when minimum/maximum with dims on StaticArrays becomes # type-stable - I0 = round.(Int, SVector(minimum(A[1, :]), minimum(A[2, :])) .- m0) .+ 1 - I1 = round.(Int, SVector(maximum(A[1, :]), maximum(A[2, :])) .- m0) .+ 1 + I0 = round.(Int, SVector(minimum(A[1, :]), minimum(A[2, :])) .- m0) .+ 1 + I1 = round.(Int, SVector(maximum(A[1, :]), maximum(A[2, :])) .- m0) .+ 1 - for I in CartesianIndex(I0[1], I0[2]):CartesianIndex(I1[1], I1[2]) + for I in CartesianIndex(I0[1], I0[2]):CartesianIndex(I1[1], I1[2]) x = SVector(Tuple(I) .- 1 .+ m0) # TODO: move to global-to-local function or inside-triangle test xloc = (A[:, SUnitRange(2, 3)] .- A[:, 1])::SArray \ @@ -110,7 +124,7 @@ function _sample(f::FeFunction) if all(xloc .>= -ε) && sum(xloc) <= 1 + ε out[fcolon..., I] .= evaluate(f, xloc) end - end + end end # convert to image indexing @@ -121,23 +135,30 @@ function _sample(f::FeFunction) return out2 end -# Interpolation + +# Mesh Function Projection Methods + +# TODO: refine interface: projection vs interpolation, unify different +# algorithms """ interprets `img` as a bilinearly interpolated continuous function and applies the default interpolation operator for the discrete function u. """ +# TODO: should be called "interpolate_nodal" function interpolate!(u::FeFunction, img::AbstractArray) f = ImageFunction(u.space.mesh, img) interpolate!(u, @inline (x; f) -> f; f) end -# TODO: refine interface: projection vs interpolation, unify different -# algorithms - project_l2_lagrange(space::FeSpace, img) = (u = FeFunction(space); project_l2_lagrange!(u, img)) +""" +standard l2 projection but using an adaptive quadrature rule based on cell size +""" +# TODO: this could be a general l2 projection, accepting a mesh function and +# some custom quadrature rule? function project_l2_lagrange!(u::FeFunction, img::AbstractArray) d = 2 # domain dimension space = u.space @@ -159,57 +180,57 @@ function project_l2_lagrange!(u::FeFunction, img::AbstractArray) # mesh cells for cell in cells(mesh) - foreach(f -> bind!(f, cell), opparams) - # cell map is assumed to be constant per cell - delmap = jacobian(elmap(mesh, cell), SA[0., 0.]) - delmapinv = inv(delmap) - intel = abs(det(delmap)) - - p = ceil(Int, diam(mesh, cell)) - qw, qx = quadrature_composite_lagrange_midpoint(p) - nqpts = length(qw) # number of quadrature points - - qphi = zeros(nrdims, nrdims, nldofs, nqpts) - dqphi = zeros(nrdims, d, nrdims, nldofs, nqpts) - for k in axes(qx, 2) - for r in 1:nrdims - qphi[r, r, :, k] .= evaluate_basis(space.element, SVector{d}(view(qx, :, k))) - dqphi[r, :, r, :, k] .= transpose(jacobian(x -> evaluate_basis(space.element, x), SVector{d}(view(qx, :, k)))) - end - end - - # quadrature points - for k in axes(qx, 2) - xhat = SVector{d}(view(qx, :, k)) - x = elmap(mesh, cell)(xhat) - opvalues = map(f -> evaluate(f, xhat), opparams) - - # local test-function dofs - for jdim in 1:nrdims, ldofj in 1:nldofs - gdofj = space.dofmap[jdim, ldofj, cell] - - phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj, k)) - dphij = SArray{Tuple{space.size..., d}}( - SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj, k)) * delmapinv) - - lv = qw[k] * l(x, phij, dphij; opvalues...) * intel - b[gdofj] += lv - - # local trial-function dofs - for idim in 1:nrdims, ldofi in 1:nldofs - gdofi = space.dofmap[idim, ldofi, cell] - - phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi, k)) - dphii = SArray{Tuple{space.size..., d}}( - SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi, k)) * delmapinv) - - av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel - push!(I, gdofi) - push!(J, gdofj) - push!(V, av) - end - end - end + foreach(f -> bind!(f, cell), opparams) + # cell map is assumed to be constant per cell + delmap = jacobian(elmap(mesh, cell), SA[0., 0.]) + delmapinv = inv(delmap) + intel = abs(det(delmap)) + + p = ceil(Int, diam(mesh, cell)) + qw, qx = quadrature_composite_lagrange_midpoint(p) + nqpts = length(qw) # number of quadrature points + + qphi = zeros(nrdims, nrdims, nldofs, nqpts) + dqphi = zeros(nrdims, d, nrdims, nldofs, nqpts) + for k in axes(qx, 2) + for r in 1:nrdims + qphi[r, r, :, k] .= evaluate_basis(space.element, SVector{d}(view(qx, :, k))) + dqphi[r, :, r, :, k] .= transpose(jacobian(x -> evaluate_basis(space.element, x), SVector{d}(view(qx, :, k)))) + end + end + + # quadrature points + for k in axes(qx, 2) + xhat = SVector{d}(view(qx, :, k)) + x = elmap(mesh, cell)(xhat) + opvalues = map(f -> evaluate(f, xhat), opparams) + + # local test-function dofs + for jdim in 1:nrdims, ldofj in 1:nldofs + gdofj = space.dofmap[jdim, ldofj, cell] + + phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj, k)) + dphij = SArray{Tuple{space.size..., d}}( + SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj, k)) * delmapinv) + + lv = qw[k] * l(x, phij, dphij; opvalues...) * intel + b[gdofj] += lv + + # local trial-function dofs + for idim in 1:nrdims, ldofi in 1:nldofs + gdofi = space.dofmap[idim, ldofi, cell] + + phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi, k)) + dphii = SArray{Tuple{space.size..., d}}( + SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi, k)) * delmapinv) + + av = qw[k] * a(x, phii, dphii, phij, dphij; opvalues...) * intel + push!(I, gdofi) + push!(J, gdofj) + push!(V, av) + end + end + end end ngdofs = ndofs(space) @@ -222,9 +243,10 @@ end project_qi_lagrange(space::FeSpace, img) = (u = FeFunction(space); project_qi_lagrange!(u, img)) -# L1-stable quasi-interpolation operator -# (https://arxiv.org/pdf/1505.06931.pdf) -# +""" +L1-stable quasi-interpolation operator +(https://arxiv.org/pdf/1505.06931.pdf) +""" # FIXME: currently only approximate quadrature by sampling on lagrange lattice # based on element size. Exact evaluation is tricky to implement (lots of # intersection handling). @@ -249,9 +271,9 @@ function project_qi_lagrange!(u::FeFunction, img) p = ceil(Int, diam(mesh, cell)) qw, qx = quadrature_composite_lagrange_midpoint(p) - for k in axes(qx, 2) - xhat = SVector{d}(view(qx, :, k)) - x = elmap(mesh, cell)(xhat) + for k in axes(qx, 2) + xhat = SVector{d}(view(qx, :, k)) + x = elmap(mesh, cell)(xhat) # size: nrdims value = evaluate(f, xhat) @@ -261,7 +283,7 @@ function project_qi_lagrange!(u::FeFunction, img) # no integration element used #u.data[gdofs] .+= qw[k] ./ area_refel .* 3 .* value * phix' u.data[gdofs] .+= qw[k] ./ area_refel * value * (12 .* phix' .- 3) - end + end gdofcount[gdofs] .+= 1 end u.data ./= gdofcount @@ -311,6 +333,8 @@ end function Base.iterate(it::PixelIterator, state...) y = iterate(it.grid, state...) isnothing(y) && return nothing + # TODO: do something more sophisticated than walking the whole cartesian + # bounding box while !_intersects(it, first(y)) y = iterate(it.grid, last(y)) isnothing(y) && return nothing @@ -356,15 +380,11 @@ function project_l2_pixel!(u::FeFunction, img) dqphi = zero(MArray{Tuple{nrdims, d, nrdims, nldofs}}) for cell in cells(mesh) - foreach(f -> bind!(f, cell), opparams) - # cell map is assumed to be constant per cell - delmap = jacobian(elmap(mesh, cell), SA[0., 0.]) - delmapinv = inv(delmap) - intel = abs(det(delmap)) - - #p = ceil(Int, diam(mesh, cell)) - #qw, qx = quadrature_composite_lagrange_midpoint(p) - #nqpts = length(qw) # number of quadrature points + foreach(f -> bind!(f, cell), opparams) + # cell map is assumed to be constant per cell + delmap = jacobian(elmap(mesh, cell), SA[0., 0.]) + delmapinv = inv(delmap) + intel = abs(det(delmap)) # pixels as quadrature points pixels = PixelIterator(mesh, cell) @@ -373,44 +393,44 @@ function project_l2_pixel!(u::FeFunction, img) # TODO: move to global-to-local function or inside-triangle test xhat = (pixels.A[:, SUnitRange(2, 3)] .- pixels.A[:, 1])::SArray \ (x .- pixels.A[:, 1]) - opvalues = map(f -> evaluate(f, xhat), opparams) + opvalues = map(f -> evaluate(f, xhat), opparams) # fill basis value/derivative buffer - for r in 1:nrdims - qphi[r, r, :] .= + for r in 1:nrdims + qphi[r, r, :] .= evaluate_basis(space.element, SVector{d}(xhat)) dqphi[r, :, r, :]::MArray{Tuple{d, nldofs}, Float64, 2, d * nldofs} .= transpose(jacobian( x -> evaluate_basis(space.element, x), SVector{d}(xhat))) - end + end - # local test-function dofs - for jdim in 1:nrdims, ldofj in 1:nldofs - gdofj = space.dofmap[jdim, ldofj, cell] + # local test-function dofs + for jdim in 1:nrdims, ldofj in 1:nldofs + gdofj = space.dofmap[jdim, ldofj, cell] - phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj)) - dphij = SArray{Tuple{space.size..., d}}( - SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj)) * delmapinv) + phij = SArray{Tuple{space.size...}}(view(qphi, :, jdim, ldofj)) + dphij = SArray{Tuple{space.size..., d}}( + SArray{Tuple{nrdims, d}}(view(dqphi, :, :, jdim, ldofj)) * delmapinv) lv = pixelweights[P] * l(x, phij, dphij; opvalues...) * intel - b[gdofj] += lv - - # local trial-function dofs - for idim in 1:nrdims, ldofi in 1:nldofs - gdofi = space.dofmap[idim, ldofi, cell] - - phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi)) - dphii = SArray{Tuple{space.size..., d}}( - SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi)) * delmapinv) - - av = pixelweights[P] * a(x, phii, dphii, phij, dphij; opvalues...) * intel - push!(I, gdofi) - push!(J, gdofj) - push!(V, av) - end - end - end + b[gdofj] += lv + + # local trial-function dofs + for idim in 1:nrdims, ldofi in 1:nldofs + gdofi = space.dofmap[idim, ldofi, cell] + + phii = SArray{Tuple{space.size...}}(view(qphi, :, idim, ldofi)) + dphii = SArray{Tuple{space.size..., d}}( + SArray{Tuple{nrdims, d}}(view(dqphi, :, :, idim, ldofi)) * delmapinv) + + av = pixelweights[P] * a(x, phii, dphii, phij, dphij; opvalues...) * intel + push!(I, gdofi) + push!(J, gdofj) + push!(V, av) + end + end + end end ngdofs = ndofs(space) @@ -418,5 +438,4 @@ function project_l2_pixel!(u::FeFunction, img) u.data .= A \ b return u - end -- GitLab