From 803a8ffee196a793f46f577de6000066eae05789 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Thu, 12 Aug 2021 13:16:05 +0200 Subject: [PATCH] implemented two methods for nearest triangle --- field.py | 89 ++++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 80 insertions(+), 9 deletions(-) diff --git a/field.py b/field.py index 07a0d84..ad3abdd 100644 --- a/field.py +++ b/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]