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 | def subdivide_quad(ind,C1,W1,CH1,PAR1,D1,A1,graded=True):
""""Turn one cell of a quadtree or octree into four or eight, respectively
Subdivides the ind-th cell of a given octree, maintaining all the adjacency information
Parameters
----------
ind : int
Index of cell to subdivide into input tree
C1 : numpy double array
Matrix of cell centers of input tree
W1 : numpy double array
Vector of cell half widths of input tree
CH1 : numpy int array
Matrix of child indeces (-1 if leaf node) of input tree
PAR1 : numpy int array
Vector of immediate parent indeces (to traverse upwards) of input tree
D1 : numpy int array
Vector of tree depths of input tree
A1 : scipy sparse.csr_matrix
Sparse node adjacency matrix of input tree, where a value of a in the (i,j) entry means that node j is to the a-th direction of i (a=1: left; a=2: right; a=3: bottom; a=4: top).
graded : bool, optional (default True)
Whether to ensure that adjacent quads only differ by one in depth or not (this is useful for numerical applications, not so much for others like position queries).
Returns
-------
C2 : numpy double array
Matrix of cell centers of output tree
W2 : numpy double array
Vector of cell half widths of output tree
CH2 : numpy int array
Matrix of child indeces (-1 if leaf node) of output tree
PAR2 : numpy int array
Vector of immediate parent indeces (to traverse upwards) of output tree
D2 : numpy int array
Vector of tree depths of output tree
A2 : scipy sparse.csr_matrix
Sparse node adjacency matrix of output tree, where a value of a in the (i,j) entry means that node j is to the a-th direction of i (a=1: left; a=2: right; a=3: bottom; a=4: top).
See Also
--------
initialize_quadtree.
Examples
--------
```python
th = 2*np.pi*np.random.rand(200,1)
# Circle
P = 0.5*np.concatenate((np.cos(th),np.sin(th)),axis=1)
# Initialize quadtree
C,W,CH,PAR,D,A = gpytoolbox.initialize_quadtree(P,graded=False,max_depth=7,min_depth=4,vmin=np.array([-1,-1]),vmax=np.array([1,1]))
# Get leaf indices
leaf_ind = gpytoolbox.quadtree_children(CH)
# Let's choose an index:
ind = leaf_ind[20]
# And let's subdivide it
C1,W1,CH1,PAR1,D1,A1 = gpytoolbox.subdivide_quad(ind,C,W,CH,PAR,D,A,graded=False)
```
"""
# If quad A lays to the <lookup[:,1]> of quad B, then children
# <lookup[:,2:end]> resulting from subdividing A will still be
# neighbors of quad B
# legend for <lookup[:,1]> is left, right, bottom, top, (back, front)
# children legend is 1: back bottom left, 2: back bottom right, 3: back top
# left, 4: back top right, 5: front bottom left, 6: front bottom right, 7: front
# top left, 8: front top right
lookup_a = np.array([[1,2,4,8,6], # A to the left of B, so we add all right children of A
[2,1,3,5,7],
[3,3,4,7,8],
[4,1,2,5,6],
[5,5,6,7,8],
[6,1,2,3,4]],dtype=int)
# (reminder 2:end columns are 1-indexed)
lookup_a[:,1:5] = lookup_a[:,1:5] - 1
# If quad A lays to the <lookup_b[:,1]> of quad B, and we subdivide both,
# then child <lookup_b[:,2]> of quad B will have child <lookup_b[:,3]> of
# quad A as a new neighbor:
lookup_b = np.array([[1,1,2],
[1,3,4],
[1,5,6],
[1,7,8],
[2,2,1],
[2,4,3],
[2,6,5],
[2,8,7],
[3,1,3],
[3,2,4],
[3,5,7],
[3,6,8],
[4,3,1],
[4,4,2],
[4,7,5],
[4,8,6],
[5,3,7],
[5,4,8],
[5,1,5],
[5,2,6],
[6,7,3],
[6,8,4],
[6,5,1],
[6,6,2]],dtype=int)
# (reminder 2:end columns are 1-indexed)
lookup_b[:,1:3] = lookup_b[:,1:3] - 1
assert(CH1[ind,1]==-1) # can't subdivide if not a leaf node
# For simplicity:
dim1 = C1.shape[1]
ncp = 2**dim1 #number of children per parent for dim-agnostic
w = W1[ind]
c = C1[ind,:]
d = D1[ind]
p = PAR1[ind]
a_ind = A1[:,ind]
num_quads = C1.shape[0]
# Easy part: add four cells new cell centers order: bottom-left,
# bottom-right, top-left, top-right
if dim1==2:
C2 = np.vstack(
(
C1,
c[None,:] + 0.25*w*np.array([[-1,-1]]),
c[None,:] + 0.25*w*np.array([[1,-1]]),
c[None,:] + 0.25*w*np.array([[-1,1]]),
c[None,:] + 0.25*w*np.array([[1,1]])
)
)
else:
C2 = np.vstack(
(
C1,
c[None,:] + 0.25*w*np.array([[-1,-1,-1]]),
c[None,:] + 0.25*w*np.array([[1,-1,-1]]),
c[None,:] + 0.25*w*np.array([[-1,1,-1]]),
c[None,:] + 0.25*w*np.array([[1,1,-1]]),
c[None,:] + 0.25*w*np.array([[-1,-1,1]]),
c[None,:] + 0.25*w*np.array([[1,-1,1]]),
c[None,:] + 0.25*w*np.array([[-1,1,1]]),
c[None,:] + 0.25*w*np.array([[1,1,1]])
)
)
# New widths
W2 = np.concatenate((W1,np.tile(np.array([0.5*w]),ncp)))
# New depths
D2 = np.concatenate((D1,np.tile(np.array([d+1],dtype=int),ncp)))
# Keep track of child indeces
CH2 = np.vstack((
CH1,
np.tile(np.array([[-1]]),(ncp,ncp))
))
CH2[ind,:] = num_quads + np.linspace(0,ncp-1,ncp,dtype=int)
#np.array([0,1,2,3],dtype=int)
# And parent indeces
PAR2 = np.concatenate((PAR1,np.tile(np.array([ind],dtype=int),ncp)))
# Now the hard part, which is the adjacency Reminder:
# (left-right-bottom-top) Effectively we are concatenating [A , B; "-B", C]
# C is always the same square 4 by 4 matrix
square_mat = csr_matrix(np.array([
[0,2,4,0],
[1,0,0,4],
[3,0,0,2],
[0,3,1,0]
]))
if dim1==3:
square_mat = vstack((
hstack((square_mat,6*eye(4))),
hstack((5*eye(4),square_mat))
))
rect_mat = csr_matrix((num_quads,ncp))
# We will loop over all neighbors of ind print("A1") print(A1.toarray())
neighbors_ind = a_ind.nonzero()[0]
# print("a_ind") print(a_ind) print("neighbors_ind") print(neighbors_ind)
where_neighbors = a_ind[neighbors_ind].toarray().squeeze()
# print("where") print(where_neighbors)
for i in range(len(neighbors_ind)):
neighbor_ind = neighbors_ind[i]
neighbor_where = where_neighbors[i]
neighbor_depth = D1[neighbor_ind]
# This will help us populate rect_mat(:,neighbor_ind) We'll build I, J
# and val. There are two options here: if the neighbor has the same or
# higher depth, it will gain two new neighbors. If not, it will gain
# only one. Maybe there's a way to combine both options but for now
# let's do it one at a time.
# Let's start with the easy case: neighbor depth is low
if neighbor_depth<=d:
# J will always be the same (neighbor_ind will gain two neighbors)
num_new_nb = 2**(dim1-1)
J = np.tile(np.array([neighbor_ind]),num_new_nb)
# Orientation will also be the same that it was
vals = np.tile(np.array([neighbor_where]),num_new_nb)
I = lookup_a[lookup_a[:,0]==neighbor_where,1:(num_new_nb+1)].squeeze()
else:
# neighbor depth is high. We need to traverse the tree to find out
# which of the four d + 1 depth children this neighbor comes from
# originally. This, combined with the neighbor_where information,
# will tell us which one new neighbor this cell has. Should use a
# look-up table with 16 possibilities J will always be the same
# (neighbor_ind will gain 1 neighbor)
J = np.array([neighbor_ind])
vals = np.array([neighbor_where])
#
# The key question is which of the four depth d+1 children do we
# originally come from. TRAVERSE!
n_ind = neighbor_ind
n_depth = neighbor_depth
n_par = PAR1[neighbor_ind]
while n_depth>(d+1):
n_ind = n_par
assert(n_depth==(D1[n_ind]+1))
n_depth = D1[n_ind]
n_par = PAR1[n_ind]
which_child = np.nonzero(CH1[n_par,:]==n_ind)[0]
# print("which_child") print(which_child) print(neighbor_where)
# print(n_ind) print(CH1[n_par,:])
assert(d==D1[n_par])
# Reminder: bottom-left, bottom-right, top-left, top-right
lookup_row = ((lookup_b[:,0]==neighbor_where) & (lookup_b[:,1]==which_child))
I = lookup_b[lookup_row,2]
rect_mat = rect_mat + csr_matrix((vals,(J,I)),shape=(num_quads,ncp))
A2 = csr_matrix(vstack((
hstack((A1,rect_mat)),
hstack((transpose_orientation(rect_mat),square_mat))
)))
#print("A2") print(A2.toarray())
if graded:
a2_ind = A2[:,ind] # from the perspective of the neighbors
# Go over all neighbors of the original quad we subdivided
neighbors_ind = a2_ind.nonzero()[0]
for i in range(len(neighbors_ind)):
is_child_2 = ((CH2[neighbors_ind[i],1]==-1))
neighbor_depth = D2[neighbors_ind[i]]
if (is_child_2 and neighbor_depth<d):
C2,W2,CH2,PAR2,D2,A2 = subdivide_quad(neighbors_ind[i],C2,W2,CH2,PAR2,D2,A2,True)
return C2,W2,CH2,PAR2,D2,A2
|