-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathmaketree.m
313 lines (256 loc) · 9.87 KB
/
maketree.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
function [tr] = maketree(rp,varargin)
%MAKETREE assemble an AABB search tree for a collection of
%(hyper-)rectangles.
% [TR] = MAKETREE(RP) returns an axis-aligned bounding-box
% (AABB) tree for a collection of d-rectangles embedded in
% R^d. The rectangles are defined as an NR-by-NDIM*2 array
% of min/max coordinates RP = [PMIN,PMAX], where PMIN =
% [A1,A2,...,ANDIM] and PMAX = [B1,B2,...,BNDIM] are the
% minimum/maximum coordinates associated with each rectan-
% gle.
%
% The resulting tree is returned as a structure containing
% an NT-by-NDIM*2 array of tree-node coordinates TR.XX =
% [PMIN,PMAX], an NT-by-2 array of tree-node indexing
% TR.II = [PI,CI], and an NT-by-1 cell array of node lists
% TR.LL. PI,CI are the parent/child pointers associated
% with each node in the tree, and TR.LL{II} is the list of
% rectangles associated with the II-th node.
%
% MAKETREE forms a binary search tree, such that each
% "leaf" node in the tree encloses a maximum number of re-
% ctangles. The tree is formed by recursively subdividing
% the bounding-box of the collection. At each division, a
% simple heuristic is used to determine a splitting axis
% and to position an axis-aligned splitting (hyper-)plane.
% The associated collection of rectangles is partitioned
% between two new child nodes. The dimensions of each node
% in the tree are selected to provide a minimal enclosure
% of the rectangles in its associated sub-tree. Tree nodes
% may overlap as a result.
%
% [...] = MAKETREE(RP,OP) also passes an additional user-
% defined options structure. OP.NOBJ = {32} is the maximum
% allowable number of rectangles per tree-node. OP.LONG =
% {.75} is a relative length tolerance for "long" rectang-
% les, such that any rectangles with RMAX(IX)-RMIN(IX) >=
% OP.LONG * (NMAX(IX)-NMIN(IX)) remain in the parent node.
% Here RMIN,RMAX are the coordinates of the rectangle,
% NMIN,NMAX are the coordinates of the enclosing node in
% the tree, and IX is the splitting axis. Nodes that beco-
% me "full" of "long" items may exceed their maximum rect-
% angle capacity. OP.VTOL = {.55} is a "volume" splitting
% criteria, designed to continue subdivision while the net
% node volume is reducing. Specifically, a node is split
% if V1+V2 <= OP.VTOL*VP, where VP is the d-dim. "volume"
% of the parent node, and V1,V2 are the volumes associated
% with its children.
%
% See also DRAWTREE, QUERYSET, MAPVERT, MAPRECT
% 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 : de2363@columbia.edu
% Last updated : 08/10/2017
tr.xx = []; tr.ii = []; tr.ll = {}; op = [];
%------------------------------ quick return on empty inputs
if (isempty(rp)), return; end
%---------------------------------------------- basic checks
if (~isnumeric(rp))
error('maketree:incorrectInputClass', ...
'Incorrect input class.');
end
%---------------------------------------------- basic checks
if (nargin < +1 || nargin > +2)
error('maketree:incorrectNoOfInputs', ...
'Incorrect number of inputs.');
end
if (ndims(rp) ~= +2 || ...
mod(size(rp,2),2)~= +0)
error('maketree:incorrectDimensions', ...
'Incorrect input dimensions.');
end
%------------------------------- extract user-defined inputs
if (nargin >= +2), op = varargin{1}; end
isoctave = exist( ...
'OCTAVE_VERSION','builtin') > 0;
if (isoctave)
%-- faster for OCTAVE with large tree block size; slower
%-- loop execution...
NOBJ = +1024 ;
else
%-- faster for MATLAB with small tree block size; better
%-- loop execution...
NOBJ = + 32 ;
end
%--------------------------------------- user-defined inputs
if (~isstruct(op))
op.nobj = NOBJ ;
op.long = .75;
op.vtol = .55;
else
%-------------------------------- bound population count
if (isfield(op,'nobj'))
if (op.nobj <= +0 )
error('Invalid options.') ;
end
else
op.nobj = NOBJ ;
end
%-------------------------------- bound "long" tolerance
if (isfield(op,'long'))
if (op.long < +.0 || op.long > +1.)
error('Invalid options.') ;
end
else
op.long = .75;
end
%-------------------------------- bound "long" tolerance
if (isfield(op,'vtol'))
if (op.vtol < +.0 || op.vtol > +1.)
error('Invalid options.') ;
end
else
op.vtol = .55;
end
end
%---------------------------------- dimensions of rectangles
nd = size(rp,2) / +2 ;
ni = size(rp,1) ;
%------------------------------------------ alloc. workspace
xl = zeros(ni*1,1*nd);
xr = zeros(ni*1,1*nd);
ii = zeros(ni*1,2);
ll = cell (ni*1,1);
ss = zeros(ni*1,1);
%------------------------------------ min & max coord. masks
lv = false(size(rp,2),1);
rv = false(size(rp,2),1);
lv((1:nd)+nd*+0) = true ;
rv((1:nd)+nd*+1) = true ;
%----------------------------------------- inflate rectangle
r0 = min(rp(:,lv),[],1) ;
r1 = max(rp(:,rv),[],1) ;
rd = repmat(r1-r0,ni,1) ;
rp(:,lv) = ...
rp(:,lv) - rd * eps^.8;
rp(:,rv) = ...
rp(:,rv) + rd * eps^.8;
%----------------------------------------- rectangle centres
rc = rp(:,lv)+rp(:,rv);
rc = rc * .5 ;
%----------------------------------------- rectangle lengths
rd = rp(:,rv)-rp(:,lv);
%------------------------------ root contains all rectangles
ll{1} = (+1:ni)' ;
%------------------------------------ indexing for root node
ii(1,1) = +0 ;
ii(1,2) = +0 ;
%------------------------------ root contains all rectangles
xl(1,:) = min(rp(:,lv),[],1);
xr(1,:) = max(rp(:,rv),[],1);
%-- main loop : divide nodes until all constraints satisfied
ss(+1) = +1; ns = +1; nn = +1;
while (ns ~= +0)
%----------------------------------- pop node from stack
ni = ss(ns) ;
ns = ns - 1 ;
%----------------------------------- push child indexing
n1 = nn + 1 ;
n2 = nn + 2 ;
%--------------------------- set of rectangles in parent
li = ll{ni} ;
%--------------------------- split plane on longest axis
dd = xr(ni,:) ...
- xl(ni,:) ;
[dd,ia] = sort(dd);
for id = nd : -1 : +1
%--------------------------- push rectangles to children
ax = ia (id) ;
mx = dd (id) ;
il = rd(li,ax) > ...
op.long * mx ;
lp = li( il) ; % "long" rectangles
ls = li(~il) ; % "short" rectangles
if (length(lp) < ...
0.5*length(ls)&& ...
length(lp) < ...
0.5 * op.nobj)
break ;
end
end
if (isempty(ls) )
%-------------------------------- partition empty, done!
continue ;
end
% select the split position: take the mean of the set of
% (non-"long") rectangle centres along axis AX
%-------------------------------------------------------
sp = sum(rc(ls,ax))/length(ls);
%---------------------------- partition based on centres
i2 = rc(ls,ax)>sp ;
l1 = ls(~i2) ; % "left" rectangles
l2 = ls( i2) ; % "right" rectangles
if (isempty(l1) || ...
isempty(l2) )
%-------------------------------- partition empty, done!
continue ;
end
%-------------------------------- finalise node position
xl(n1,:) = ...
min(rp(l1,lv),[],1) ;
xr(n1,:) = ...
max(rp(l1,rv),[],1) ;
xl(n2,:) = ...
min(rp(l2,lv),[],1) ;
xr(n2,:) = ...
max(rp(l2,rv),[],1) ;
%--------------------------- push child nodes onto stack
if (length(li) <= op.nobj )
vi = prod(xr(ni,:) ... % upper d-dim "vol."
- xl(ni,:) ) ;
v1 = prod(xr(n1,:) ... % lower d-dim "vol."
- xl(n1,:) ) ;
v2 = prod(xr(n2,:) ...
- xl(n2,:) ) ;
if (v1+v2 < op.vtol*vi)
%-------------------------------- parent--child indexing
ii(n1,1) = ni ;
ii(n2,1) = ni ;
ii(ni,2) = n1 ;
%-------------------------------- finalise list indexing
ll{ni,1} = lp ;
ll{n1,1} = l1 ;
ll{n2,1} = l2 ;
ss(ns+1) = n1 ;
ss(ns+2) = n2 ;
ns = ns+2; nn = nn+2;
end
else
%-------------------------------- parent--child indexing
ii(n1,1) = ni ;
ii(n2,1) = ni ;
ii(ni,2) = n1 ;
%-------------------------------- finalise list indexing
ll{ni,1} = lp ;
ll{n1,1} = l1 ;
ll{n2,1} = l2 ;
ss(ns+1) = n1 ;
ss(ns+2) = n2 ;
ns = ns+2; nn = nn+2;
end
end
%----------------------------------------------- trim alloc.
xl = xl(1:nn,:) ;
xr = xr(1:nn,:) ;
ii = ii(1:nn,:) ;
ll = ll(1:nn,:) ;
%----------------------------------------------- pack struct
tr.xx =[xl, xr] ;
tr.ii = ii ;
tr.ll = ll ;
end