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