1# distutils: include_dirs = LIB_DIR
2# distutils: libraries = STD_LIBS
3"""
4Geometry selection routines.
5
6
7
8
9"""
10
11
12import numpy as np
13
14cimport cython
15cimport numpy as np
16cimport oct_visitors
17from cython cimport floating
18from libc.math cimport sqrt
19from libc.stdlib cimport free, malloc
20
21from yt.utilities.lib.bitarray cimport ba_get_value, ba_set_value
22from yt.utilities.lib.fnv_hash cimport c_fnv_hash as fnv_hash
23from yt.utilities.lib.fp_utils cimport fclip, fmax, fmin, iclip, imax, imin
24from yt.utilities.lib.geometry_utils cimport (
25    bounded_morton_dds,
26    decode_morton_64bit,
27    encode_morton_64bit,
28    morton_neighbors_coarse,
29    morton_neighbors_refined,
30)
31from yt.utilities.lib.grid_traversal cimport sampler_function, walk_volume
32from yt.utilities.lib.volume_container cimport VolumeContainer
33
34from .oct_container cimport Oct, OctreeContainer
35from .oct_visitors cimport cind
36
37
38cdef extern from "math.h":
39    double exp(double x) nogil
40    float expf(float x) nogil
41    long double expl(long double x) nogil
42    double floor(double x) nogil
43    double ceil(double x) nogil
44    double fmod(double x, double y) nogil
45    double log2(double x) nogil
46    long int lrint(double x) nogil
47    double fabs(double x) nogil
48
49# use this as an epsilon test for grids aligned with selector
50# define here to avoid the gil later
51cdef np.float64_t grid_eps = np.finfo(np.float64).eps
52grid_eps = 0.0
53
54cdef inline np.float64_t dot(np.float64_t* v1,
55                             np.float64_t* v2) nogil:
56    return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
57
58cdef inline np.float64_t norm(np.float64_t* v) nogil:
59    return sqrt(dot(v, v))
60
61# These routines are separated into a couple different categories:
62#
63#   * Routines for identifying intersections of an object with a bounding box
64#   * Routines for identifying cells/points inside a bounding box that
65#     intersect with an object
66#   * Routines that speed up some type of geometric calculation
67
68# First, bounding box / object intersection routines.
69# These all respect the interface "dobj" and a set of left_edges, right_edges,
70# sometimes also accepting level and mask information.
71
72@cython.boundscheck(False)
73@cython.wraparound(False)
74@cython.cdivision(True)
75def convert_mask_to_indices(np.ndarray[np.uint8_t, ndim=3, cast=True] mask,
76            int count, int transpose = 0):
77    cdef int i, j, k, cpos
78    cdef np.ndarray[np.int64_t, ndim=2] indices
79    indices = np.zeros((count, 3), dtype='int64')
80    cpos = 0
81    for i in range(mask.shape[0]):
82        for j in range(mask.shape[1]):
83            for k in range(mask.shape[2]):
84                if mask[i, j, k] == 1:
85                    if transpose == 1:
86                        indices[cpos, 0] = k
87                        indices[cpos, 1] = j
88                        indices[cpos, 2] = i
89                    else:
90                        indices[cpos, 0] = i
91                        indices[cpos, 1] = j
92                        indices[cpos, 2] = k
93                    cpos += 1
94    return indices
95
96
97@cython.boundscheck(False)
98@cython.wraparound(False)
99@cython.cdivision(True)
100cdef _mask_fill(np.ndarray[np.float64_t, ndim=1] out,
101                np.int64_t offset,
102                np.ndarray[np.uint8_t, ndim=3, cast=True] mask,
103                np.ndarray[floating, ndim=3] vals):
104    cdef np.int64_t count = 0
105    cdef int i, j, k
106    for i in range(mask.shape[0]):
107        for j in range(mask.shape[1]):
108            for k in range(mask.shape[2]):
109                if mask[i, j, k] == 1:
110                    out[offset + count] = vals[i, j, k]
111                    count += 1
112    return count
113
114def mask_fill(np.ndarray[np.float64_t, ndim=1] out,
115              np.int64_t offset,
116              np.ndarray[np.uint8_t, ndim=3, cast=True] mask,
117              np.ndarray vals):
118    if vals.dtype == np.float32:
119        return _mask_fill[np.float32_t](out, offset, mask, vals)
120    elif vals.dtype == np.float64:
121        return _mask_fill[np.float64_t](out, offset, mask, vals)
122    else:
123        raise RuntimeError
124
125@cython.cdivision(True)
126@cython.boundscheck(False)
127@cython.wraparound(False)
128def points_in_cells(
129        np.float64_t[:] cx,
130        np.float64_t[:] cy,
131        np.float64_t[:] cz,
132        np.float64_t[:] dx,
133        np.float64_t[:] dy,
134        np.float64_t[:] dz,
135        np.float64_t[:] px,
136        np.float64_t[:] py,
137        np.float64_t[:] pz):
138    # Take a list of cells and particles and calculate which particles
139    # are enclosed within one of the cells.  This is used for querying
140    # particle fields on clump/contour objects.
141    # We use brute force since the cells are a relatively unordered collection.
142
143    cdef int p, c, n_p, n_c
144    cdef np.ndarray[np.uint8_t, ndim=1, cast=True] mask
145
146    n_p = px.size
147    n_c = cx.size
148    mask = np.zeros(n_p, dtype="bool")
149
150    for p in range(n_p):
151        for c in range(n_c):
152            if (fabs(px[p] - cx[c]) <= 0.5 * dx[c] and
153                fabs(py[p] - cy[c]) <= 0.5 * dy[c] and
154                fabs(pz[p] - cz[c]) <= 0.5 * dz[c]):
155                mask[p] = True
156                break
157
158    return mask
159
160include "_selection_routines/selector_object.pxi"
161include "_selection_routines/point_selector.pxi"
162include "_selection_routines/sphere_selector.pxi"
163include "_selection_routines/region_selector.pxi"
164include "_selection_routines/cut_region_selector.pxi"
165include "_selection_routines/disk_selector.pxi"
166include "_selection_routines/cutting_plane_selector.pxi"
167include "_selection_routines/slice_selector.pxi"
168include "_selection_routines/ortho_ray_selector.pxi"
169include "_selection_routines/ray_selector.pxi"
170include "_selection_routines/data_collection_selector.pxi"
171include "_selection_routines/ellipsoid_selector.pxi"
172include "_selection_routines/grid_selector.pxi"
173include "_selection_routines/octree_subset_selector.pxi"
174include "_selection_routines/indexed_octree_subset_selector.pxi"
175include "_selection_routines/always_selector.pxi"
176include "_selection_routines/compose_selector.pxi"
177include "_selection_routines/halo_particles_selector.pxi"
178include "_selection_routines/boolean_selectors.pxi"
179