diff --git a/src/image.jl b/src/image.jl
index e9b315e8bebbf16deba8ca16bbdd8c6e67a1195b..609c9a7e258d63f04197cbfa8a6d641d7b4e0ca3 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