implemented two methods for nearest triangle

field.py

@@ -923,6 +923,9 @@ class Features3d:
         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])
 
@@ -931,23 +934,91 @@ class Features3d:
         # zset = set(x.data for x in sorted(treez.at(0.0)))
         return
 
-    # @staticmethod
-    # def ray_triangle_intersect(self):
-
-    def ray_triangle_intersect(R,dR,A,B,C): 
+    @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    = -np.dot(dR,N)
+        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      =  np.dot(E2,DAO)*invdet
-        v      = -np.dot(E1,DAO)*invdet
+        u =  (E2*DAO).sum(axis=-1)*invdet
+        v = -(E1*DAO).sum(axis=-1)*invdet
         # Intersection point is R+t*dR
-        t      =  np.dot(AO,N)*invdet
-        return (t,t*np.dot(N,dR)) if (np.abs(det) >= 1e-6 and u >= 0.0 and v >= 0.0 and (u+v) <= 1.0) else None
+        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 
+           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]
+        '''
+        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
+        # Now we computed all the criteria to determine a hit
+        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 (False,np.nan,None)
+        # Compute the intersection distance: intersection occurs at R+t*dR
+        t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx]
+        # Get closest hit
+        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]