Skip to content

subdivide_quad

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:

Name Type Description Default
ind int

Index of cell to subdivide into input tree

required
C1 numpy double array

Matrix of cell centers of input tree

required
W1 numpy double array

Vector of cell half widths of input tree

required
CH1 numpy int array

Matrix of child indeces (-1 if leaf node) of input tree

required
PAR1 numpy int array

Vector of immediate parent indeces (to traverse upwards) of input tree

required
D1 numpy int array

Vector of tree depths of input tree

required
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).

required
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).

True

Returns:

Name Type Description
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:

TODO

Source code in src/gpytoolbox/subdivide_quad.py
  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
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
    --------
    TODO
    """

    # 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