From 1a8afb97d0048e4c1a51ab6dfea67bdefd64e8ff Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Thu, 12 Aug 2021 21:53:20 +0200 Subject: [PATCH] first implementation of find_features(). At the moment it only returns whether or not a point is inside --- field.py | 223 +++++++++++++++++++++++++++---------------------------- 1 file changed, 111 insertions(+), 112 deletions(-) diff --git a/field.py b/field.py index ad3abdd..8a0da86 100644 --- a/field.py +++ b/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],:] - V = points[faces[:,2],:] - W = points[faces[:,3],:] - cn = np.cross(V-U,W-U) - cc = (U+V+W)/3 + A = points[faces[:,1],:] + B = points[faces[:,2],:] + C = points[faces[:,3],:] + cn = np.cross(B-A,C-A) + 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,99 +880,7 @@ class Features3d: '''Returns an array with volumes of all features.''' return np.add.reduceat(self._cell_volumes,self._offset[:-1]) - def find_feature(self): - # 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 - - t = time() - - xmin = np.amin(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) - zmax = np.amax(self._points[self._faces[:,1:],2],axis=1) - - center = np.stack((0.5*(xmax+xmin),0.5*(zmax+zmin)),axis=1) - radius = np.sqrt(2.0)*np.maximum( - np.amax(0.5*(xmax-xmin)), - np.amax(0.5*(zmax-zmin))) - print('Preparation for KD tree in',time()-t) - - t = time() - kd = spatial.KDTree(center,leafsize=100,compact_nodes=False,copy_data=False,balanced_tree=False) - # kd = spatial.KDTree(center,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True) - # kd = spatial.KDTree(self.points[:,1:],leafsize=20,compact_nodes=False,copy_data=False,balanced_tree=False) - # 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()): - # map_[face].append(ii) - # faces_connected_to_vertex = {} - # for face_index,face in enumerate(self._faces[:,1:]): - # 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)) - - t = time() - 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]) - - - # 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 - - @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 = -(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 - # Intersection point is R+t*dR - 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): + 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 @@ -991,6 +899,102 @@ class Features3d: 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) + # xmax = np.amax(self._points[self._faces[:,1:],0],axis=1) + # zmin = np.amin(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)) + print(time()-t) + del xmin,xmax,zmin,zmax + print('Preparation for KD tree in',time()-t) + + t = time() + kd = spatial.KDTree(center,leafsize=100,compact_nodes=False,copy_data=False,balanced_tree=False) + # kd = spatial.KDTree(center,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True) + # kd = spatial.KDTree(self.points[:,1:],leafsize=20,compact_nodes=False,copy_data=False,balanced_tree=False) + # kd = spatial.KDTree(self._points[:,1:],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True) + print('KD tree built in',time()-t) + + # t = time() + # map_ = self._points.shape[0]*[[]] + # for ii,face in enumerate(self._faces[:,1:].ravel()): + # map_[face].append(ii) + # faces_connected_to_vertex = {} + # for face_index,face in enumerate(self._faces[:,1:]): + # 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)) + + t = time() + 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 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) @@ -1002,23 +1006,18 @@ class Features3d: # 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 + # 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) + 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 + 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] - # 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 + return hit_idx,t,N[hit_idx] def clean_points(self,report=False): nfaces_ = self._faces.shape[0]