1cdef class EllipsoidSelector(SelectorObject):
2    cdef public np.float64_t vec[3][3]
3    cdef public np.float64_t mag[3]
4    cdef public np.float64_t center[3]
5
6    def __init__(self, dobj):
7        cdef int i
8        _ensure_code(dobj.center)
9        _ensure_code(dobj._e0)
10        _ensure_code(dobj._e1)
11        _ensure_code(dobj._e2)
12        _ensure_code(dobj._A)
13        _ensure_code(dobj._B)
14        _ensure_code(dobj._C)
15        for i in range(3):
16            self.center[i] = dobj.center[i]
17            self.vec[0][i] = dobj._e0[i]
18            self.vec[1][i] = dobj._e1[i]
19            self.vec[2][i] = dobj._e2[i]
20        self.mag[0] = dobj._A
21        self.mag[1] = dobj._B
22        self.mag[2] = dobj._C
23
24    @cython.boundscheck(False)
25    @cython.wraparound(False)
26    @cython.cdivision(True)
27    cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3]) nogil:
28        return self.select_point(pos)
29
30    @cython.boundscheck(False)
31    @cython.wraparound(False)
32    @cython.cdivision(True)
33    cdef int select_point(self, np.float64_t pos[3]) nogil:
34        cdef np.float64_t dot_evec[3]
35        cdef np.float64_t dist
36        cdef int i, j
37        dot_evec[0] = dot_evec[1] = dot_evec[2] = 0
38        # Calculate the rotated dot product
39        for i in range(3): # axis
40            dist = self.periodic_difference(pos[i], self.center[i], i)
41            for j in range(3):
42                dot_evec[j] += dist * self.vec[j][i]
43        dist = 0.0
44        for i in range(3):
45            dist += (dot_evec[i] * dot_evec[i])/(self.mag[i] * self.mag[i])
46        if dist <= 1.0: return 1
47        return 0
48
49    @cython.boundscheck(False)
50    @cython.wraparound(False)
51    @cython.cdivision(True)
52    cdef int select_sphere(self, np.float64_t pos[3], np.float64_t radius) nogil:
53        # this is the sphere selection
54        cdef int i
55        cdef np.float64_t dist, dist2_max, dist2 = 0
56        for i in range(3):
57            dist = self.periodic_difference(pos[i], self.center[i], i)
58            dist2 += dist * dist
59        dist2_max = (self.mag[0] + radius) * (self.mag[0] + radius)
60        if dist2 <= dist2_max:
61            return 1
62        return 0
63
64    @cython.boundscheck(False)
65    @cython.wraparound(False)
66    @cython.cdivision(True)
67    cdef int select_bbox(self, np.float64_t left_edge[3],
68                               np.float64_t right_edge[3]) nogil:
69        # This is the sphere selection
70        cdef int i
71        cdef np.float64_t box_center, relcenter, closest, dist, edge, dist_max
72        if left_edge[0] <= self.center[0] <= right_edge[0] and \
73           left_edge[1] <= self.center[1] <= right_edge[1] and \
74           left_edge[2] <= self.center[2] <= right_edge[2]:
75            return 1
76        # http://www.gamedev.net/topic/335465-is-this-the-simplest-sphere-aabb-collision-test/
77        dist = 0
78        for i in range(3):
79            box_center = (right_edge[i] + left_edge[i])/2.0
80            relcenter = self.periodic_difference(box_center, self.center[i], i)
81            edge = right_edge[i] - left_edge[i]
82            closest = relcenter - fclip(relcenter, -edge/2.0, edge/2.0)
83            dist += closest * closest
84        dist_max = self.mag[0] * self.mag[0]
85        if dist <= dist_max:
86            return 1
87        return 0
88
89    @cython.boundscheck(False)
90    @cython.wraparound(False)
91    @cython.cdivision(True)
92    cdef int select_bbox_edge(self, np.float64_t left_edge[3],
93                               np.float64_t right_edge[3]) nogil:
94        # This is the sphere selection
95        cdef int i
96        cdef np.float64_t box_center, relcenter, closest, farthest, cdist, fdist, edge
97        if left_edge[0] <= self.center[0] <= right_edge[0] and \
98           left_edge[1] <= self.center[1] <= right_edge[1] and \
99           left_edge[2] <= self.center[2] <= right_edge[2]:
100            fdist = 0
101            for i in range(3):
102                edge = right_edge[i] - left_edge[i]
103                box_center = (right_edge[i] + left_edge[i])/2.0
104                relcenter = self.periodic_difference(
105                    box_center, self.center[i], i)
106                farthest = relcenter + fclip(relcenter, -edge/2.0, edge/2.0)
107                fdist += farthest*farthest
108                if fdist >= self.mag[0]**2: return 2
109            return 1
110        # http://www.gamedev.net/topic/335465-is-this-the-simplest-sphere-aabb-collision-test/
111        cdist = 0
112        fdist = 0
113        for i in range(3):
114            box_center = (right_edge[i] + left_edge[i])/2.0
115            relcenter = self.periodic_difference(box_center, self.center[i], i)
116            edge = right_edge[i] - left_edge[i]
117            closest = relcenter - fclip(relcenter, -edge/2.0, edge/2.0)
118            farthest = relcenter + fclip(relcenter, -edge/2.0, edge/2.0)
119            cdist += closest * closest
120            fdist += farthest * farthest
121            if cdist > self.mag[0]**2: return 0
122        if fdist < self.mag[0]**2:
123            return 1
124        else:
125            return 2
126
127    def _hash_vals(self):
128        return (("vec[0][0]", self.vec[0][0]),
129                ("vec[0][1]", self.vec[0][1]),
130                ("vec[0][2]", self.vec[0][2]),
131                ("vec[1][0]", self.vec[1][0]),
132                ("vec[1][1]", self.vec[1][1]),
133                ("vec[1][2]", self.vec[1][2]),
134                ("vec[2][0]", self.vec[2][0]),
135                ("vec[2][1]", self.vec[2][1]),
136                ("vec[2][2]", self.vec[2][2]),
137                ("mag[0]", self.mag[0]),
138                ("mag[1]", self.mag[1]),
139                ("mag[2]", self.mag[2]),
140                ("center[0]", self.center[0]),
141                ("center[1]", self.center[1]),
142                ("center[2]", self.center[2]))
143
144    def _get_state_attnames(self):
145        return ("mag", "center", "vec")
146
147
148ellipsoid_selector = EllipsoidSelector
149