From 6e3564c3c157fed6c63527ea6828d37ddc56c58c Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 13 Aug 2021 02:04:15 +0200 Subject: [PATCH] implemented feature search. it should work now and perform well. --- field.py | 129 ++++++++++++++++++++++++++----------------------------- 1 file changed, 61 insertions(+), 68 deletions(-) diff --git a/field.py b/field.py index 45e11e1..024cc5d 100644 --- a/field.py +++ b/field.py @@ -678,6 +678,10 @@ class Features3d: self.triangulate(contour_method=contour_method, cellvol_normal_component=cellvol_normal_component, report=report) + # Set some state variable + self._kdtree = None + self._kdaxis = None + self._kdradius = None # Drop input if requested: this may free memory in case it was copied/temporary. if not keep_input: self._input = None return @@ -870,7 +874,37 @@ 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 build_kdtree(self,kdaxis=0,leafsize=100,compact_nodes=False,balanced_tree=False,report=False): + '''Builds a KD-tree for feature search.''' + from scipy import spatial + if kdaxis==0: + min1 = np.amin(self._vertices[:,:,1],axis=1) + max1 = np.amax(self._vertices[:,:,1],axis=1) + min2 = np.amin(self._vertices[:,:,2],axis=1) + max2 = np.amax(self._vertices[:,:,2],axis=1) + elif kdaxis==1: + min1 = np.amin(self._vertices[:,:,0],axis=1) + max1 = np.amax(self._vertices[:,:,0],axis=1) + min2 = np.amin(self._vertices[:,:,2],axis=1) + max2 = np.amax(self._vertices[:,:,2],axis=1) + elif kdaxis==2: + min1 = np.amin(self._vertices[:,:,0],axis=1) + max1 = np.amax(self._vertices[:,:,0],axis=1) + min2 = np.amin(self._vertices[:,:,1],axis=1) + max2 = np.amax(self._vertices[:,:,1],axis=1) + else: + raise ValueError("Invalid ray axis.") + center = np.stack((0.5*(max1+min1),0.5*(max2+min2)),axis=1) + radius = 0.5*np.amax(np.sqrt((max1-min1)**2+(max2-min2)**2)) + self._kdtree = spatial.KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes, + copy_data=False,balanced_tree=balanced_tree) + self._kdaxis = kdaxis + self._kdradius = radius + if report: + print('[Features3d.build_kdtree]',end=' ') + print('KD-Tree built for axis {}.'.format(axis)) + + def inside_feature(self,coords): '''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 @@ -889,80 +923,39 @@ 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 + coords = np.array(coords) + assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array." - # pts = np.random.random((10000,3)) + if not self._kdtree: + self.build_kdtree() + if self._kdaxis==0: query_axis = [1,2] + elif self._kdaxis==1: query_axis = [0,2] + elif self._kdaxis==2: query_axis = [0,1] - t = time() + cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius) - xmin = np.amin(self._vertices[:,:,0],axis=1) - xmax = np.amax(self._vertices[:,:,0],axis=1) - zmin = np.amin(self._vertices[:,:,2],axis=1) - zmax = np.amax(self._vertices[:,:,2],axis=1) - 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._vertices[bla[ii],0,:], - self._vertices[bla[ii],1,:], - self._vertices[bla[ii],2,:]) + N = coords.shape[0] + raydir = np.zeros((3,),dtype=self._vertices.dtype) + raydir[self._kdaxis] = 1.0 + in_feat = np.empty((N,),dtype=np.int) + for ii in range(N): + hit_idx,t,N = Features3d.ray_triangle_intersection(coords[ii],raydir, + self._vertices[cand[ii],0,:], + self._vertices[cand[ii],1,:], + self._vertices[cand[ii],2,:]) if hit_idx is None: - is_inside[ii] = False + in_feat[ii] = -1 else: idx = np.argmin(np.abs(t)) - is_inside[ii] = t[idx]*(N[idx,:]*dR).sum(axis=-1)>0 - print('inside check',time()-t__) + if t[idx]*(N[idx,:]*raydir).sum(axis=-1)>0: + in_feat[ii] = self.feature_by_cell(cand[ii][idx]) + else: + in_feat[ii] = -1 + return in_feat - # 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 + def feature_by_cell(self,idx_cell): + '''Gets feature ID for a given cell.''' + return np.argmax(idx_cell>self._offset) @staticmethod def ray_triangle_intersection(R,dR,A,B,C):