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...')
-        A = points[faces[:,1],:]
+        U = points[faces[:,1],:]
-        B = points[faces[:,2],:]
+        V = points[faces[:,2],:]
-        C = points[faces[:,3],:]
+        W = points[faces[:,3],:]
-        cn = np.cross(B-A,C-A)
+        cn = np.cross(U-W,V-U)
-        cc = (A+B+C)/3
+        cc = (U+V+W)/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,51 +880,23 @@ class Features3d:
         '''Returns an array with volumes of all features.'''
         return np.add.reduceat(self._cell_volumes,self._offset[:-1])
 
-    def find_feature(self,pts):
+    def find_feature(self):
-        '''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 = 0.5*np.amax(np.sqrt((xmax-xmin)**2+(zmax-zmin)**2))
+        radius = np.sqrt(2.0)*np.maximum(
-        print(time()-t)
+            np.amax(0.5*(xmax-xmin)),
-        del xmin,xmax,zmin,zmax
+            np.amax(0.5*(zmax-zmin)))
         print('Preparation for KD tree in',time()-t)
 
         t = time()
@@ -934,6 +906,8 @@ 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()):
@@ -943,81 +917,20 @@ 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()
-        bla = kd.query_ball_point(pts[:,[0,2]],radius)
+        kd.query_ball_point(query,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
+        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:
-                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]