Skip to content
Snippets Groups Projects
Commit 83c0116d authored by Stephan Hilb's avatar Stephan Hilb
Browse files

refactor

parent c63a8eb3
No related branches found
No related tags found
No related merge requests found
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))
......@@ -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
......@@ -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
......@@ -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).
......@@ -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
......@@ -362,10 +386,6 @@ function project_l2_pixel!(u::FeFunction, img)
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
# pixels as quadrature points
pixels = PixelIterator(mesh, cell)
for P in pixels
......@@ -418,5 +438,4 @@ function project_l2_pixel!(u::FeFunction, img)
u.data .= A \ b
return u
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment