-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfindtria.m
320 lines (265 loc) · 9.9 KB
/
findtria.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
function [tp,tj,tr] = findtria(pp,tt,pj,varargin)
%FINDTRIA spatial queries for collections of d-simplexes.
% [TP,TI] = FINDTRIA(PP,TT,PJ) finds the set of simple-
% xes that intersect with a given spatial query. Simplexes
% are specified via the vertex array PP = [X1,X2,...,XN]
% and the indexing array TT = [T1,T2,...,TM], such that
% the vertex positions for the II-th simplex are the poin-
% ts [PP(TT(II,1),:),PP(TT(II,2),:),...,PP(TT(II,M),:)].
%
% Simplexes are NOT required to form a conforming triangu-
% lation. Specifically, non-delaunay, non-convex and even
% overlapping configurations are supported. Multiple matc-
% hes may be returned if the collection is overlapping.
%
% A set of intersecting simplexes is returned for each
% query point in PI, such that the II-th point is associa-
% ted with the simplexes TI(TP(II,1):TP(II,2)). Unenclosed
% points have TP(II,1)==+0.
%
% In general, query points may be matched to multiple sim-
% plexes, but in cases when single matches are guaranteed,
% or if only a single match is desired, the following ret-
% urns a singly-matched indexing array that is consistent
% with MATLAB's existing point-location routines:
%
% [tp,tj] = findtria(pp,tt,pj) ;
% ti = nan(size(tp,1),1);
% in = tp(:,1) > +0;
% ti(in) = tj(tp(in,+1));
%
% [TP,TI,TR] = FINDTRIA(PP,TT,PI) additionally returns the
% supporting aabb-tree used internally to compute the que-
% ry. If the underlying collection [PP,TT] is static, the
% tree TR may be passed to subsequent calls, via [...] =
% FINDTRIA(PP,TT,PJ,TR). This syntax may lead to improved
% performance, especially when the number of simplexes
% is large w.r.t. the number of query points. Note that in
% such cases the underlying simplexes are NOT permitted to
% change between calls, or erroneous results may be retur-
% ned. Additional parameters used to control the creation
% of the underlying aabb-tree may also be passed via [...]
% = FINDTRIA(PP,TT,PI,TR,OP). See MAKETREE for additional
% information.
%
% See also MAKETREE, EXCHANGE, QUERYSET
% Please see the following for additional information:
%
% Darren Engwirda, "Locally-optimal Delaunay-refinement &
% optimisation-based mesh generation". Ph.D. Thesis, Scho-
% ol of Mathematics and Statistics, Univ. of Sydney, 2014:
% http://hdl.handle.net/2123/13148
% Darren Engwirda : 2014 --
% Email : [email protected]
% Last updated : 10/03/2018
tp = []; tj = []; tr = []; op = [];
filename = mfilename('fullpath');
filepath = fileparts( filename );
addpath([filepath,'/aabb-tree']);
%---------------------------------------------- basic checks
if (nargin < +3 || nargin > +6)
error('queryset:incorrectNumInputs', ...
'Incorrect number of inputs.');
end
%------------------------------- fast return on empty inputs
if (isempty(pj)), return; end
%------------------------------- extract user-defined inputs
if (nargin >= +4), tr = varargin{1}; end
if (nargin >= +5), op = varargin{2}; end
%---------------------------------------------- basic checks
if (~isnumeric(pp) || ~isnumeric(tt) || ...
~isnumeric(pj) )
error('findtria:incorrectInputClass', ...
'Incorrect input class.') ;
end
%---------------------------------------------- basic checks
if (ndims(pp) ~= +2 || size(pp,2) < +2 ...
|| size(pp,2) > size(tt,2))
error('findtria:incorrectDimensions', ...
'Incorrect input dimensions.');
end
if (ndims(tt) ~= +2 || size(tt,2) < +3)
error('findtria:incorrectDimensions', ...
'Incorrect input dimensions.');
end
%------------------------------- test query array for VERT's
if (ndims(pj) ~= +2 || ...
size(pj,2) ~= size(pp,2) )
error('findtria:incorrectDimensions', ...
'Incorrect input dimensions.');
end
%---------------------------------------------- basic checks
if (~isempty(tr) && ~isstruct(tr) )
error('findtria:incorrectInputClass', ...
'Incorrect input class.') ;
end
if (~isempty(op) && ~isstruct(op) )
error('findtria:incorrectInputClass', ...
'Incorrect input class.') ;
end
if (isempty(tr))
%------------------------------ compute aabb's for triangles
bi = pp(tt(:,1),:); bj = pp(tt(:,1),:);
for ii = 2 : size(tt,2)
bi = min(bi,pp(tt(:,ii),:)) ;
bj = max(bj,pp(tt(:,ii),:)) ;
end
bb = [bi,bj];
%------------------------------ compute aabb-tree for aabb's
if (isempty(op)) % scale against |pj|
op.nobj = ceil(size(tt,1) / ...
size(pj,1)) * +4 ;
op.nobj = max( +32,op.nobj) ;
op.nobj = min(+256,op.nobj) ;
end
tr = maketree(bb,op); % compute aabb-tree
end
%------------------------------ compute tree-to-vert mapping
tm = mapvert (tr,pj);
%------------------------------ compute vert-to-tria queries
x0 = min(pp,[],1);
x1 = max(pp,[],1);
rt = prod(x1 - x0) * eps^.8 ;
[ti,ip,tj] = ...
queryset(tr,tm,@triakern,pj,pp,tt,rt) ;
%------------------------------ re-index onto full obj. list
tp = zeros(size(pj,1),2);
tp( :,2) = -1 ;
if (isempty(ti)), return; end
tp(ti,:) = ip ;
end
function [ip,it] = triakern(pk,tk,pi,pp,tt,rt)
%TESTPTS compute the "point-tria" matches within a tile.
mp = length(pk); mt = length(tk);
switch (size(tt,2))
case 3
%------------------------------------ pts in 2-simplexes
pk = pk.' ;
pk = pk(ones(mt,1),:);
pk = pk(:);
tk = tk(:,ones(1,mp));
tk = tk(:);
in = intria2( ...
pp,tt(tk,:),pi(pk,:),rt);
ip = pk(in) ;
it = tk(in) ;
case 4
%------------------------------------ pts in 3-simplexes
pk = pk.' ;
pk = pk(ones(mt,1),:);
pk = pk(:);
tk = tk(:,ones(1,mp));
tk = tk(:);
in = intria3( ...
pp,tt(tk,:),pi(pk,:),rt);
ip = pk(in) ;
it = tk(in) ;
otherwise
%------------------------------------ pts in d-simplexes
[il,jl] = intrian( ...
pp,tt(tk,:),pi(pk,:));
ip = pk(il(:));
it = tk(jl(:));
end
end
function [in] = intria2(pp,tt,pi,rt)
%INTRIA2 returns TRUE for points enclosed by 2-simplexes.
t1 = tt(:,1); t2 = tt(:,2) ;
t3 = tt(:,3);
vi = pp(t1,:) - pi ;
vj = pp(t2,:) - pi ;
vk = pp(t3,:) - pi ;
%------------------------------- compute sub-volume about PI
aa = zeros(size(tt,1),3) ;
aa(:,1) =(vi(:,1).*vj(:,2) - ...
vj(:,1).*vi(:,2) ) ;
aa(:,2) =(vj(:,1).*vk(:,2) - ...
vk(:,1).*vj(:,2) ) ;
aa(:,3) =(vk(:,1).*vi(:,2) - ...
vi(:,1).*vk(:,2) ) ;
%------------------------------- PI is internal if same sign
rt = rt ^ 2 ;
in = aa(:,1).*aa(:,2) >= -rt ...
& aa(:,2).*aa(:,3) >= -rt ...
& aa(:,3).*aa(:,1) >= -rt ;
end
function [in] = intria3(pp,tt,pi,rt)
%INTRIA3 returns TRUE for points enclosed by 3-simplexes.
t1 = tt(:,1); t2 = tt(:,2) ;
t3 = tt(:,3); t4 = tt(:,4) ;
v1 = pi - pp(t1,:) ;
v2 = pi - pp(t2,:) ;
v3 = pi - pp(t3,:) ;
v4 = pi - pp(t4,:) ;
%------------------------------- compute sub-volume about PI
aa = zeros(size(tt,1),4) ;
aa(:,1) = ...
+v1(:,1).*(v2(:,2).*v3(:,3) - ...
v2(:,3).*v3(:,2) ) ...
-v1(:,2).*(v2(:,1).*v3(:,3) - ...
v2(:,3).*v3(:,1) ) ...
+v1(:,3).*(v2(:,1).*v3(:,2) - ...
v2(:,2).*v3(:,1) ) ;
aa(:,2) = ...
+v1(:,1).*(v4(:,2).*v2(:,3) - ...
v4(:,3).*v2(:,2) ) ...
-v1(:,2).*(v4(:,1).*v2(:,3) - ...
v4(:,3).*v2(:,1) ) ...
+v1(:,3).*(v4(:,1).*v2(:,2) - ...
v4(:,2).*v2(:,1) ) ;
aa(:,3) = ...
+v2(:,1).*(v4(:,2).*v3(:,3) - ...
v4(:,3).*v3(:,2) ) ...
-v2(:,2).*(v4(:,1).*v3(:,3) - ...
v4(:,3).*v3(:,1) ) ...
+v2(:,3).*(v4(:,1).*v3(:,2) - ...
v4(:,2).*v3(:,1) ) ;
aa(:,4) = ...
+v3(:,1).*(v4(:,2).*v1(:,3) - ...
v4(:,3).*v1(:,2) ) ...
-v3(:,2).*(v4(:,1).*v1(:,3) - ...
v4(:,3).*v1(:,1) ) ...
+v3(:,3).*(v4(:,1).*v1(:,2) - ...
v4(:,2).*v1(:,1) ) ;
%------------------------------- PI is internal if same sign
rt = rt ^ 2 ;
in = aa(:,1).*aa(:,2) >= -rt ...
& aa(:,1).*aa(:,3) >= -rt ...
& aa(:,1).*aa(:,4) >= -rt ...
& aa(:,2).*aa(:,3) >= -rt ...
& aa(:,2).*aa(:,4) >= -rt ...
& aa(:,3).*aa(:,4) >= -rt ;
end
function [ii,jj] = intrian(pp,tt,pi)
%INTRIAN return a list of points and enclosing "n"-simplexes.
[np,pd] = size(pi);
[nt,td] = size(tt);
%---------------- coefficient matrices for barycentric coord.
mm = zeros(pd,pd,nt);
for id = +1 : pd
for jd = +1 : pd
mm(id,jd,:) = pp(tt(:,jd),id) - ...
pp(tt(:,td),id) ;
end
end
%---------------- solve linear systems for barycentric coord.
xx = zeros(pd,np,nt);
vp = zeros(pd,np,+1);
for ti = +1 : nt
%---------------------------------------- form rhs coeff.
for id = +1 : pd
vp(id,:) = ...
pi(:,id) - pp(tt(ti,td),id) ;
end
%---------------------------- actually faster to call LU-
[ll,uu] = lu(mm(:,:,ti));
%---------------------------- and then forward/back solve
xx(:,:,ti) = uu\(ll\vp);
end
%-------------------- PI is internal if coord. have same sign
in = all(xx >= +.0-(eps^.8),1) & ...
sum(xx,1) <= +1.+(eps^.8) ;
%-------------------- find lists of matching points/simplexes
[ii,jj] = ...
find(reshape(in, [np,nt])) ;
end