Line data Source code
1 : module BarycentricCoordinates
2 :
3 : using LinearAlgebra
4 : using StaticArrays
5 :
6 : #
7 : # Define simplex for barycentric interpolation
8 : #
9 :
10 : """
11 : Simplex for barycentric interpolation as a callable object.
12 :
13 : f = baryint(...) :: Simplex
14 : y = f(λ)
15 : y = f(λ₁, λ₂, ...)
16 :
17 : See also [`baryint`](@ref), [`MapBaryCoords`](@ref), [`mapbarycoords`](@ref)
18 : """
19 : struct Simplex{M, T, N, L} # NOTE: Must specify L for (or dynamic dispatch!)!
20 139 : X::SMatrix{M, N, T, L}
21 : end
22 :
23 13 : baryint(a::V, b::V) where {V} =
24 : Simplex{length(a), eltype(V), 2, 2length(a)}(SMatrix{length(a), 2}(a..., b...))
25 :
26 126 : baryint(a::V, b::V, c::V) where {V} =
27 : Simplex{length(a), eltype(V), 3, 3length(a)}(SMatrix{length(a), 3}(a..., b..., c...))
28 :
29 : """
30 : f = baryint(x1, x2, ...)
31 : f = baryint((x1, x2, ...))
32 :
33 : y = f(λ)
34 : y = f(λ₁, λ₂, ...)
35 :
36 : Get callable object `f` that applied barycentric interpolation
37 : `[x1 x2 ...] * λ` with `sum(λ) == 1`.
38 :
39 : !!! warning "precondition"
40 : Requires `sum(λ) == 1`. This is *not* checked!
41 :
42 : !!! note
43 : [`baryint`](@ref) supports only line segments and triangles,
44 : both in 2d and 3d.
45 :
46 : See also [`Simplex`](@ref), [`mapbarycoords`](@ref), [`barycoords`](@ref)
47 : """ baryint
48 :
49 3677 : (b::Simplex{M, T, N})(λ::SVector{N, T}) where {M, T, N} = b.X * λ
50 :
51 7 : (b::Simplex{M, T, N})(λ::AbstractVector) where {M, T, N} =
52 : b(SVector{N, T}(λ...))
53 :
54 0 : (b::Simplex)(λ) = b(SVector(λ...))
55 :
56 0 : (b::Simplex)(λ::R, λs...) where {R<:Number} = b.X * SVector(λ, λs...)
57 :
58 : #
59 : # Get affine map that computes barycentric coordinates
60 : #
61 :
62 : """
63 : Affine map that computes barycentric coordinates.
64 :
65 : b = mapbarycoords(x1, x2, ...)
66 :
67 : f = baryint(x1, x2, ...) :: Simplex
68 : b = mapbarycoords(f) :: MapBaryCoords
69 : λ = b(x)
70 : y = f(λ)
71 :
72 : λ = barycoords(f, x)
73 :
74 : @assert sum(λ) ≈ 1
75 :
76 : In case of a `n`-simplex `f` in `m>n` dimensions, the map "assumes"
77 : that is `x` is in the affine subspace defined by `f`. If this is not
78 : the case, there exist no barycentric coordinates, and the map computes
79 : barycentric coordinates for a *projection* into the affine
80 : subspace. Note that in this case, we have `f(λ) != x`.
81 :
82 : In case the simplex `f` is *degenerated*, there exist no barycentric
83 : coordinates. In this case, the map computes a "best fit":
84 :
85 : !!! note
86 : Mapping `x` to barycentric coordinates `λ` *never fails*, i.e.,
87 : yields *always finite* barycentric coordinates `λ`. Consequently,
88 : one may have `f(λ) != x`.
89 :
90 : See also [`mapbarycoords`](@ref), [`baryint`](@ref), [`Simplex`](@ref),
91 : [`barycoords`](@ref)
92 : """
93 : struct MapBaryCoords{M, N, D, T, L}
94 2252 : A::SMatrix{M, N, T, L}
95 : x0::SVector{D, T}
96 : end
97 :
98 3793 : function (bc::MapBaryCoords)(x)
99 3793 : λ = bc.A * (x - bc.x0)
100 3793 : SVector(λ..., 1-sum(λ))
101 : end
102 :
103 : """
104 : f = baryint(x1, x2, ...) :: Simplex
105 : λ = barycoords(f, x)
106 :
107 : @assert sum(λ) ≈ 1
108 :
109 : Evaluate barycentric coordinates `λ` for `x` in `f`.
110 :
111 : !!! note
112 : [`mapbarycoords`](@ref) constructs an affine map from `x` to `λ`.
113 : This is more efficient, if multiple `x` are mapped for the same
114 : simplex.
115 :
116 : See also [`mapbarycoords`](@ref)
117 : """
118 2170 : barycoords(s::Simplex, x) = mapbarycoords(s)(x)
119 :
120 0 : mapbarycoords(x::V, y::V, z...) where {V<:AbstractVector} =
121 : mapbarycoords(baryint(x, y, z...))
122 :
123 110 : function mapbarycoords(s::Simplex{N, T, 2}) where {N, T}
124 110 : x0 = s.X[:, 2]
125 110 : a = s.X[:, 1] - x0
126 :
127 110 : len2 = dot(a, a)
128 :
129 110 : (len2 > 16eps()) || (len2 = one(T))
130 :
131 110 : MapBaryCoords{1, N, N, T, N}(a / len2, x0)
132 : end
133 :
134 1480 : function mapbarycoords(s::Simplex{2, T, 3}) where {T}
135 1480 : x0 = s.X[:, 3]
136 1480 : A = hcat(s.X[:, 1]-x0, s.X[:, 2]-x0)
137 :
138 1480 : invA = inv(A)
139 1500 : all(isfinite, invA) || (invA = pinv(A))
140 :
141 1480 : MapBaryCoords{2, 2, 2, T, 4}(invA, x0)
142 : end
143 :
144 662 : function mapbarycoords(s::Simplex{3, T, 3}) where {T}
145 662 : x0 = s.X[:, 3]
146 662 : u = s.X[:, 1] - x0
147 662 : v = s.X[:, 2] - x0
148 :
149 662 : w = abs.(u × v)
150 1986 : i = argmax(w)
151 :
152 662 : if w[i] < eps()
153 60 : i = argmin(abs(u[i]) + abs(v[i]) for i in 1:length(u))
154 : end
155 :
156 662 : A, P =
157 : if i == 1
158 199 : @SMatrix([ u[2] v[2]; u[3] v[3] ]), SMatrix{2, 3, T}(0, 0, 1, 0, 0, 1)
159 463 : elseif i == 2
160 190 : @SMatrix([ u[1] v[1]; u[3] v[3] ]), SMatrix{2, 3, T}(1, 0, 0, 0, 0, 1)
161 : else
162 1125 : @SMatrix([ u[1] v[1]; u[2] v[2] ]), SMatrix{2, 3, T}(1, 0, 0, 1, 0, 0)
163 : end
164 :
165 : # TODO: test projections (add tests below)
166 : # TODO: update/assemble matrix inv(A) in-place (instead of inv(A)*P)
167 :
168 662 : invA =
169 : if w[i] < eps()
170 40 : pinv(A)
171 : else
172 1304 : inv(A)
173 : end
174 :
175 662 : MapBaryCoords{2, 3, 3, T, 6}(invA * P, x0)
176 : end
177 :
178 : """
179 : b = mapbarycoords(x1, x2, ...)
180 :
181 : f = baryint(x1, x2, ...) :: Simplex
182 : b = mapbarycoords(f) :: MapBaryCoords
183 : λ = b(x)
184 : y = f(λ)
185 :
186 : λ = barycoords(f, x)
187 :
188 : Get affine map [`MapBaryCoords`](@ref) that maps `x` to barycentric
189 : coordinates `λ`.
190 :
191 : !!! note
192 : For computing barycentric coordinates for *one single* point `x`,
193 : [`barycoords`](@ref) may be more convenient.
194 :
195 : In case of a `n`-simplex `f` in `m>n` dimensions, the map "assumes"
196 : that is `x` is in the affine subspace defined by `f`. If this is not
197 : the case, there exist no barycentric coordinates, and the map computes
198 : barycentric coordinates for a *projection* into the affine
199 : subspace. Note that in this case, we have `f(λ) != x`.
200 :
201 : In case the simplex `f` is *degenerated*, there exist no barycentric
202 : coordinates. In this case, the map computes a "best fit":
203 :
204 : !!! note
205 : Mapping `x` to barycentric coordinates `λ` *never fails*, i.e.,
206 : yields *always finite* barycentric coordinates `λ`. Consequently,
207 : one may have `f(λ) != x`.
208 :
209 : See [`MapBaryCoords`](@ref), [`barycoords`](@ref), [`baryint`](@ref)
210 : """ mapbarycoords
211 :
212 : end
|
