changed computation of triangle normal: should give the same result as before. Implemented Möller-Trumbore ray-triangle intersection calculation

field.py

@@ -774,7 +774,7 @@ class Features3d:
         U = points[faces[:,1],:]
         V = points[faces[:,2],:]
         W = points[faces[:,3],:]
-        cn = np.cross(U-W,V-U)
+        cn = np.cross(V-U,W-U)
         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]
@@ -917,8 +917,7 @@ 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)
@@ -932,6 +931,27 @@ 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(r0,dr,v0,v1,v2):
+        v0v1 = v1-v0
+        v0v2 = v2-v0
+        pvec = np.cross(dr,v0v2)
+        det = np.dot(v0v1,pvec)
+        if np.abs(det)<1e-6: 
+            return float('-inf')
+        invDet = 1.0 / det
+        tvec = r0-v0
+        u = np.dot(tvec,pvec)*invDet
+        if u<0 or u>1: 
+            return float('-inf')
+        qvec = np.cross(tvec,v0v1)
+        v = np.dot(dr,qvec)*invDet
+        if v<0 or u+v>1: 
+            return float('-inf')
+        return np.dot(v0v2,qvec)*invDet
+
     def clean_points(self,report=False):
         nfaces_ = self._faces.shape[0]
         ind,inv = np.unique(self._faces[:,1:],return_inverse=True)