Line data Source code
1 : """
2 : Context for supporting point location queryies on [`triangulated`](@ref)
3 : [`geometry`](@ref).
4 : """
5 : struct PointLocation{V, F, H, E, N, T}
6 :
7 : mesh::Mesh{V, F, H, E}
8 : vertices::Vector{VHnd}
9 : faces::Vector{FHnd}
10 : tree::DistanceQueries.SimplexTree{N, 3, T}
11 :
12 0 : PointLocation{V, F, H, E, N, T}(
13 : mesh::Mesh{V, F, H, E},
14 : vertices::Vector{VHnd},
15 : faces::Vector{FHnd},
16 : tree::DistanceQueries.SimplexTree{N, 3, T}) where {V, F, H, E, N, T} =
17 : new(mesh, vertices, faces, tree)
18 : end
19 :
20 : """
21 : pl = pointlocation(g[; maxpoints, maxdepth])
22 :
23 : d, f, λ = nearestface(pl, y)
24 : d, v = nearestvertex(pl, y)
25 :
26 : py = projection(pl, y)
27 :
28 : The implementation uses
29 : [`DistanceQueries`](https://pages.vc.cs.ovgu.de/distancequeries.jl/).
30 :
31 : This method takes a "snapshot" of the geometry that is stored
32 : independently of the `mesh` in a kd-tree.
33 :
34 : !!! warning
35 : Handles returned by [`nearestvertex`](@ref) and [`nearestface`](@ref)
36 : become *invalid* if the mesh changes after calling [`pointlocation`](@ref).
37 :
38 : !!! note
39 : Consider using [`compact`](@ref), i.e., ensure [`iscontiguous`](@ref)
40 : before calling [`pointlocation`](@ref)!
41 :
42 : See also [`geometry`](@ref), [`triangulated`](@ref), [`nearestface`](@ref),
43 : [`nearestvertex`](@ref), [`nearest`](@ref), [`projection`](@ref)
44 : """
45 0 : pointlocation(g::AbstractGeometry;
46 : maxpoints::Int = 10, maxdepth::Int =32) =
47 : pointlocation(g, _istriangulation(g); maxpoints, maxdepth)
48 :
49 : # get mapping of VHnd or [] if identity
50 0 : function _pl_vertices(mesh::Mesh) :: Vector{VHnd}
51 0 : if iscontiguous(mesh, VHnd)
52 0 : VHnd[]
53 : else
54 0 : collect(PMesh.vertices(mesh))
55 : end
56 : end
57 :
58 : # get mapping of FHnd or [] if identity
59 0 : function _pl_faces(mesh::Mesh) :: Vector{FHnd}
60 0 : if iscontiguous(mesh, FHnd)
61 0 : FHnd[]
62 : else
63 0 : collect(PMesh.faces(mesh))
64 : end
65 : end
66 :
67 : # map vertex positions from _pl_vertices
68 0 : function _pl_vpos(g::AbstractGeometry, vertices::Vector{VHnd})
69 0 : mesh = _mesh(g)
70 0 : @massert iscontiguous(mesh, VHnd) == isempty(vertices)
71 :
72 0 : x = vec(_vpos(g))
73 0 : VT = typeof(x)
74 0 : T = eltype(VT)
75 :
76 0 : if iscontiguous(mesh, VHnd)
77 0 : copy(x) :: VT
78 : else
79 : # TODO: make this a method
80 0 : xs = Vector{T}(undef, nv(mesh))
81 0 : for (i, v) in enumerate(vertices)
82 0 : xs[i] = x[v]
83 0 : end
84 0 : xs ::VT
85 : end
86 : end
87 :
88 : # map vertex indices in triangle table ts as specified by vertices
89 0 : function _pl_map_vidx!(mesh::Mesh, ts::Vector{SVector{3, Int}}, vertices::Vector{VHnd})
90 0 : @massert iscontiguous(mesh, VHnd) == isempty(vertices)
91 :
92 0 : if !iscontiguous(mesh, VHnd)
93 0 : rvidx = zeros(Int, nv(mesh))
94 :
95 0 : for (i, j) in enumerate(vertices)
96 0 : rvidx[Int(j)] = i
97 0 : end
98 :
99 0 : ts = map(ts) do ijk
100 0 : map(vi -> rvidx[vi], ijk)
101 : end
102 : end
103 : end
104 :
105 0 : function pointlocation(g::AbstractGeometry, ::IsTriangulation;
106 : maxpoints::Int, maxdepth::Int)
107 0 : mesh = _mesh(g)
108 :
109 : # get handles of triangular faces and vertices
110 : # (possibly not a contiguous sequence of handles)
111 0 : faces = _pl_faces(mesh)
112 0 : vertices = _pl_vertices(mesh)
113 :
114 : # get triangulation
115 0 : ts = triangletable(mesh)
116 :
117 : # get positions for used vertices (map by `vertices` if not contiguous)
118 0 : x = _pl_vpos(g, vertices)
119 :
120 : # update vertex indices in `ts` by `vertices` (if not contiguous)
121 0 : _pl_map_vidx!(mesh, ts, vertices)
122 :
123 0 : pointstree = kdtree(x; maxpoints, maxdepth)
124 0 : tree = kdtree(pointstree, scmatrix(ts))
125 :
126 0 : pointlocation(mesh, vertices, faces, tree)
127 : end
128 :
129 0 : function pointlocation(mesh::Mesh{V, F, H, E},
130 : vertices::Vector{VHnd},
131 : faces::Vector{FHnd},
132 : tree::DistanceQueries.SimplexTree{N, 3, T}) where {V, F, H, E, N, T}
133 0 : PointLocation{V, F, H, E, N, T}(mesh, vertices, faces, tree)
134 : end
135 :
136 : """
137 : g = triangulated(geometry(mesh, x))
138 : pl = pointlocation(g)
139 :
140 : d, f, λ = nearestface(pl, y)
141 :
142 : Given querypoint `y`, compute distance `d` to `mesh`, with the
143 : [`projection`](@ref) point in face `f` at barycentric cooordinates
144 : `λ`.
145 :
146 : For an *planar* mesh in 2d, `d` is zero if `y` is within the
147 : triangulation. If `y` is outside, `λ` specifies a point on a boundary
148 : edge.
149 :
150 : !!! warning "requires"
151 : The mesh must not have changed sicne calling [`pointlocation`](@ref)
152 :
153 : See also [`geometry`](@ref), [`triangulated`](@ref),
154 : [`nearestvertex`](@ref), [`nearest`](@ref), [`projection`](@ref)
155 : """
156 0 : function nearestface(pl::PointLocation, x)
157 0 : d, i, λ = DistanceQueries.nearest(pl.tree, x)
158 0 : f = isempty(pl.faces) ? FHnd(i) : pl.faces[i]
159 0 : d, f, λ
160 : end
161 :
162 : """
163 : g = triangulated(geometry(mesh, x))
164 : pl = pointlocation(g)
165 :
166 : d, v = nearestvertex(pl, y)
167 :
168 : Given querypoint `y`, compute the nearest vertex `v` and the
169 : distance `d` to this vertex.
170 :
171 : !!! warning "requires"
172 : The mesh must not have changed since calling [`pointlocation`](@ref)
173 :
174 : See also [`geometry`](@ref), [`triangulated`](@ref),
175 : [`nearestface`](@ref), [`nearest`](@ref), [`projection`](@ref)
176 : """
177 0 : function nearestvertex(pl::PointLocation, x)
178 0 : d, i = DistanceQueries.nearest(parent(pl.tree), x)
179 0 : d, isempty(pl.vertices) ? VHnd(i) : pl.vertices[i]
180 : end
181 :
182 : """
183 : Synomym for [`nearestface`](@ref)
184 : """
185 0 : DistanceQueries.nearest(pl::PointLocation, x) = nearestface(pl, x)
186 :
187 : """
188 : g = triangulated(geometry(mesh, x))
189 : pl = pointlocation(g)
190 :
191 : py = projection(pl, y)
192 :
193 : Given a point `y`, compute its projection onto the `mesh`.
194 :
195 : !!! note
196 : This method does not require any data from the `mesh`!
197 : It relies only on the "snapshot" taken by
198 : [`pointlocation`](@ref).
199 :
200 : See also [`geometry`](@ref), [`triangulated`](@ref),
201 : [`nearestface`](@ref)
202 : """
203 0 : function projection(pl::PointLocation, x)
204 0 : _, i, λ = DistanceQueries.nearest(pl.tree, x)
205 0 : point(pl.tree, i, λ)
206 : end
|
