first implementation of find_features(). At the moment it only returns whether or not a point is inside

field.py

@@ -771,11 +771,11 @@ class Features3d:
         # of the normal has to be chosen which defaults to the z-component and is set by 
         # 'cellvol_normal_component'.
         if report: print('[Features3d.triangulate] calculating area and volume per cell...')
-        U = points[faces[:,1],:]
+        A = points[faces[:,1],:]
-        V = points[faces[:,2],:]
+        B = points[faces[:,2],:]
-        W = points[faces[:,3],:]
+        C = points[faces[:,3],:]
-        cn = np.cross(V-U,W-U)
+        cn = np.cross(B-A,C-A)
-        cc = (U+V+W)/3
+        cc = (A+B+C)/3
         area = 0.5*np.sqrt(np.square(cn).sum(axis=1))
         vol  = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component]
         # Now the label is known per cell. We only need to find all cells with the same label
@@ -880,99 +880,7 @@ class Features3d:
         '''Returns an array with volumes of all features.'''
         return np.add.reduceat(self._cell_volumes,self._offset[:-1])
 
-    def find_feature(self):
+    def find_feature(self,pts):
-        # coords = np.array(coords)
-        # assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array."
-        from time import time
-        from scipy import spatial
-
-        t = time()
-
-        xmin = np.amin(self._points[self._faces[:,1:],0],axis=1)
-        xmax = np.amax(self._points[self._faces[:,1:],0],axis=1)
-        zmin = np.amin(self._points[self._faces[:,1:],2],axis=1)
-        zmax = np.amax(self._points[self._faces[:,1:],2],axis=1)
-
-        center = np.stack((0.5*(xmax+xmin),0.5*(zmax+zmin)),axis=1)
-        radius = np.sqrt(2.0)*np.maximum(
-            np.amax(0.5*(xmax-xmin)),
-            np.amax(0.5*(zmax-zmin)))
-        print('Preparation for KD tree in',time()-t)
-
-        t = time()
-        kd = spatial.KDTree(center,leafsize=100,compact_nodes=False,copy_data=False,balanced_tree=False)
-        # kd = spatial.KDTree(center,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
-        # kd = spatial.KDTree(self.points[:,1:],leafsize=20,compact_nodes=False,copy_data=False,balanced_tree=False)
-        # kd = spatial.KDTree(self._points[:,1:],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
-        print('KD tree built in',time()-t)
-
-        query = np.random.random((10000,2))
-
-        # t = time()
-        # map_ = self._points.shape[0]*[[]]
-        # for ii,face in enumerate(self._faces[:,1:].ravel()):
-        #     map_[face].append(ii)
-        # faces_connected_to_vertex = {}
-        # for face_index,face in enumerate(self._faces[:,1:]):
-        #     for vertex_index in face:
-        #         faces_connected_to_vertex.setdefault(vertex_index,[]).append(face_index)
-        # print('Inverted cells-faces',time()-t)
-        # radius = np.sqrt(np.sum(self.spacing[1:]**2))
-        
-        t = time()
-        kd.query_ball_point(query,radius)
-        print('KD tree query',time()-t)
-
-        #    The provided vertices do not need to form a closed surface! This means they can be 
-        #    filtered before being passed to this function.
-
-        # print(radius)
-        # print(len(kd.query_ball_point((0.5,0.0),radius,return_sorted=True)),self._points.shape[0])
-
-
-        # yset = set(x.data for x in sorted(treey.at(0.0)))
-        # zset = set(x.data for x in sorted(treez.at(0.0)))
-        return
-
-    @staticmethod
-    def ray_triangle_intersection(R,dR,A,B,C):
-        '''Implements the Möller–Trumbore ray-triangle intersection algorithm. I modified the 
-           formulation of 
-           https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d
-           because it computes the cell normal on the way, which is needed to determine the 
-           direction of the hit, i.e. from the inside or outside.
-           Input:
-                R, dR: origin and direction of the ray as (3,) numpy arrays
-                A,B,C: vertices of N triangles as (N,3) numpy arrays
-           Returns:
-                is_hit:  did the ray hit any triangle? [bool]
-                t:       parameter to determine intersection point (x = R+t*dR) [float]
-                hit_dir: direction from which the triangle was hit, from inward/outward = +1,-1 [int]
-        '''
-        E1     = B-A
-        E2     = C-A
-        N      = np.cross(E1,E2)
-        det    = -(dR*N).sum(axis=-1) # dot product
-        det[np.abs(det)<1e-6] = np.nan # mask to avoid numpy runtime errors
-        invdet = 1.0/det
-        AO     = R-A
-        DAO    = np.cross(AO,dR)
-        # u,v,1-u-v are the barycentric coordinates of intersection
-        u =  (E2*DAO).sum(axis=-1)*invdet
-        v = -(E1*DAO).sum(axis=-1)*invdet
-        # Intersection point is R+t*dR
-        t = (AO*N).sum(axis=-1)*invdet
-        # Boolean array indicating hits
-        is_hit = np.logical_and(
-                   np.logical_and(
-                     np.logical_and(
-                       np.isfinite(det),u>=0.0),
-                     v>=0.0),
-                  (u+v)<=1.0)
-        return is_hit,t
-
-    @staticmethod
-    def nearest_triangle(R,dR,A,B,C):
         '''Find nearest triangle from point R in direction ±dR and its orientation w.r.t dR. 
            The algorithm determines the closest intersection of a ray cast from the origin of the 
            point in the specified direction (with positiv and negative sign). The ray is supposed to 
@@ -991,6 +899,102 @@ class Features3d:
                 t:         parameter to determine intersection point (x = R+t*dR) [float]
                 hit_dir:   direction from which the triangle was hit, from inward/outward = +1,-1 [int]
         '''
+        # coords = np.array(coords)
+        # assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array."
+        from time import time
+        from scipy import spatial
+        from .jit import minmax
+
+        # pts = np.random.random((10000,3))
+
+        print(self._points[self._faces[:,1:],0].shape)
+        t = time()
+
+        # xmin = np.amin(self._points[self._faces[:,1:],0],axis=1)
+        # xmax = np.amax(self._points[self._faces[:,1:],0],axis=1)
+        # zmin = np.amin(self._points[self._faces[:,1:],2],axis=1)
+        # zmax = np.amax(self._points[self._faces[:,1:],2],axis=1)
+        # print(time()-t)
+        xmin,xmax = minmax(self._points[self._faces[:,1:],0])
+        zmin,zmax = minmax(self._points[self._faces[:,1:],2])
+        print(time()-t)
+        # print(xmin.shape,xmin2.shape)
+        # print(np.all(np.isclose(xmin,xmin2)))
+
+        center = np.stack((0.5*(xmax+xmin),0.5*(zmax+zmin)),axis=1)
+        radius = 0.5*np.amax(np.sqrt((xmax-xmin)**2+(zmax-zmin)**2))
+        print(time()-t)
+        del xmin,xmax,zmin,zmax
+        print('Preparation for KD tree in',time()-t)
+
+        t = time()
+        kd = spatial.KDTree(center,leafsize=100,compact_nodes=False,copy_data=False,balanced_tree=False)
+        # kd = spatial.KDTree(center,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
+        # kd = spatial.KDTree(self.points[:,1:],leafsize=20,compact_nodes=False,copy_data=False,balanced_tree=False)
+        # kd = spatial.KDTree(self._points[:,1:],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
+        print('KD tree built in',time()-t)
+
+        # t = time()
+        # map_ = self._points.shape[0]*[[]]
+        # for ii,face in enumerate(self._faces[:,1:].ravel()):
+        #     map_[face].append(ii)
+        # faces_connected_to_vertex = {}
+        # for face_index,face in enumerate(self._faces[:,1:]):
+        #     for vertex_index in face:
+        #         faces_connected_to_vertex.setdefault(vertex_index,[]).append(face_index)
+        # print('Inverted cells-faces',time()-t)
+        # radius = np.sqrt(np.sum(self.spacing[1:]**2))
+        
+        t = time()
+        bla = kd.query_ball_point(pts[:,[0,2]],radius)
+        print('KD tree query',time()-t)
+
+        t__ = time()
+        Npts = pts.shape[0]
+        print(pts.shape,Npts)
+        dR = (0,1,0)
+        is_inside = np.empty((Npts,),dtype=bool)
+        print(is_inside.shape)
+        for ii in range(Npts):
+            hit_idx,t,N = Features3d.ray_triangle_intersection(pts[ii],dR,
+                                    self._points[self._faces[bla[ii],1]],
+                                    self._points[self._faces[bla[ii],2]],
+                                    self._points[self._faces[bla[ii],3]])
+            if hit_idx is None:
+                is_inside[ii] = False
+            else:
+                idx = np.argmin(np.abs(t))
+                is_inside[ii] = t[idx]*(N[idx,:]*dR).sum(axis=-1)>0
+        print('inside check',time()-t__)
+
+        # print(is_inside)
+        #    The provided vertices do not need to form a closed surface! This means they can be 
+        #    filtered before being passed to this function.
+
+        # print(radius)
+        # print(len(kd.query_ball_point((0.5,0.0),radius,return_sorted=True)),self._points.shape[0])
+
+
+        # yset = set(x.data for x in sorted(treey.at(0.0)))
+        # zset = set(x.data for x in sorted(treez.at(0.0)))
+        return is_inside
+
+    @staticmethod
+    def ray_triangle_intersection(R,dR,A,B,C):
+        '''Implements the Möller–Trumbore ray-triangle intersection algorithm. I modified the 
+           formulation of 
+           https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d
+           because it computes the cell normal on the way, which is needed to determine the 
+           direction of the hit, i.e. from the inside or outside.
+           Input:
+                R, dR: origin and direction of the ray as (3,) numpy arrays
+                A,B,C: vertices of N triangles as (N,3) numpy arrays
+           Returns:
+                hit_idx: index of the input triangles which returned a hit. [(Nhit,) int]
+                t:       parameter to determine intersection point (x = R+t*dR) [(Nhit,) float]
+                N:       normal vector of triangles which were hit [(Nhit,3) float]
+           All returned values are None if not hit occured.
+        '''
         E1     = B-A
         E2     = C-A
         N      = np.cross(E1,E2)
@@ -1002,23 +1006,18 @@ class Features3d:
         # u,v,1-u-v are the barycentric coordinates of intersection
         u      =  (E2*DAO).sum(axis=-1)*invdet
         v      = -(E1*DAO).sum(axis=-1)*invdet
-        # Now we computed all the criteria to determine a hit
+        # Boolean array indicating hits
         is_hit = np.logical_and(
-                   np.logical_and(
+                np.logical_and(
-                     np.logical_and(
+                    np.logical_and(
-                       np.isfinite(det),u>=0.0),
+                    np.isfinite(det),u>=0.0),
-                     v>=0.0),
+                    v>=0.0),
-                   (u+v)<=1.0)
+                (u+v)<=1.0)
         hit_idx = np.flatnonzero(is_hit)
-        if len(hit_idx)==0: return (False,np.nan,None)
+        if len(hit_idx)==0: return (None,None,None)
-        # Compute the intersection distance: intersection occurs at R+t*dR
+        # Intersection point is R+t*dR
         t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx]
-        # Get closest hit
+        return hit_idx,t,N[hit_idx]
-        subidx = np.argmin(np.abs(t))
-        idx = hit_idx[subidx]
-        # Compute the direction of this hit
-        hit_dir = np.sign(t[subidx]*(N[idx,:]*dR).sum(axis=-1))
-        return idx,t[subidx],hit_dir
 
     def clean_points(self,report=False):
         nfaces_ = self._faces.shape[0]