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(U-W,V-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,23 +880,51 @@ 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):
+        '''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 
+           be sent in a direction perpendicular to any boundary. Periodicity does not matter here.
+           The ray-triangle intersection is implemented by the Möller–Trumbore algorithm (see method
+           'ray_triangle_intersection()' for details). Whether a hit occured from the interior
+           or the exterior is determined by the direction of the surface normal of the triangle 
+           and the direction of the ray.
+           This method is very similar to 'ray_triangle_intersection()', but takes advantage of
+           the fact that only the nearest intersection is to be returned.
+           Input:
+                R, dR: the point to be checked and the direction of the ray as (3,) numpy arrays
+                A,B,C: vertices of N triangles as (N,3) numpy arrays
+           Returns:
+                is_inside: is point inside? [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]
+        '''
         # 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)
+        # xmin = np.amin(self._points[self._faces[:,1:],0],axis=1)
-        xmax = np.amax(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)
+        # zmin = np.amin(self._points[self._faces[:,1:],2],axis=1)
-        zmax = np.amax(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 = np.sqrt(2.0)*np.maximum(
+        radius = 0.5*np.amax(np.sqrt((xmax-xmin)**2+(zmax-zmin)**2))
-            np.amax(0.5*(xmax-xmin)),
+        print(time()-t)
-            np.amax(0.5*(zmax-zmin)))
+        del xmin,xmax,zmin,zmax
         print('Preparation for KD tree in',time()-t)
 
         t = time()
@@ -906,8 +934,6 @@ class Features3d:
         # 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()):
@@ -917,20 +943,81 @@ class Features3d:
         #     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))
-        radius = np.sqrt(np.sum(self.spacing[1:]**2))
         
         t = time()
-        kd.query_ball_point(query,radius)
+        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
+        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)
+        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
+        # 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)
+        hit_idx = np.flatnonzero(is_hit)
+        if len(hit_idx)==0: return (None,None,None)
+        # Intersection point is R+t*dR
+        t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx]
+        return hit_idx,t,N[hit_idx]
 
     def clean_points(self,report=False):
         nfaces_ = self._faces.shape[0]