From 928cc9c9d21b8c085be51cbb8089a58fe28e8fc0 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Thu, 19 Aug 2021 09:29:43 +0200 Subject: [PATCH] it works now --- field.py | 394 +++++++++++++++++++++++-------------------------------- 1 file changed, 166 insertions(+), 228 deletions(-) diff --git a/field.py b/field.py index 3c99349..5f03cfd 100644 --- a/field.py +++ b/field.py @@ -650,7 +650,7 @@ def gaussian_filter_umean_channel(array,spacing,sigma,truncate=4.0): class Features3d: def __init__(self,input,threshold,origin,spacing,periodicity, invert=False,has_ghost=False,keep_input=False, - contour_method='flying_edges',cellvol_normal_component=2, + contour_method='flying_edges', report=False): assert len(origin)==3, "'origin' must be of length 3" assert len(spacing)==3, "'spacing' must be of length 3" @@ -681,7 +681,6 @@ class Features3d: self._input[1:-1,1:-1,1:-1] = np.pad(sign_invert*input,pw,mode='wrap') # Triangulate self.triangulate(contour_method=contour_method, - cellvol_normal_component=cellvol_normal_component, report=report) # Set some state variable self._kdtree = None @@ -694,14 +693,13 @@ class Features3d: @classmethod def from_field(cls,fld,threshold,periodicity,invert=False,has_ghost=False, keep_input=False,contour_method='flying_edges', - cellvol_normal_component=2,report=False): + correct_normals=True,report=False): return cls(fld.data,threshold,fld.origin,fld.spacing,periodicity, invert=invert,has_ghost=has_ghost,keep_input=keep_input, contour_method=contour_method, - cellvol_normal_component=cellvol_normal_component, report=report) - def triangulate(self,contour_method='flying_edges',cellvol_normal_component=2,report=False): + def triangulate(self,contour_method='flying_edges',report=False): import pyvista as pv import vtk from scipy import ndimage, spatial @@ -723,8 +721,9 @@ class Features3d: if report: print('[Features3d.triangulate] computing isocontour using {}...'.format(contour_method)) contour = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False,compute_gradients=True) assert contour.is_all_triangles(), "Contouring produced non-triangle cells." - points = contour.points - faces = contour.faces.reshape(-1,4) + points = contour.points + faces = contour.faces.reshape(-1,4) + gradients = contour.point_arrays['Gradients'][faces[:,1],:] # Python for loops are horribly slow, but for loops are the easiest way to perform # some of the necessary check. Therefore let's define some functions to be jit compiled # by numba to speed things up. @@ -734,9 +733,9 @@ class Features3d: assert points.shape[1]==3 assert len(bounds)==6 assert len(tolerance)==3 - ncells = faces.shape[0] - output = np.full((ncells,),-1,dtype=np.int8) - for ii in range(ncells): + nfaces = faces.shape[0] + output = np.full((nfaces,),-1,dtype=np.int8) + for ii in range(nfaces): for jj in range(3): if (points[faces[ii,1],jj]=0) - self._offset = np.zeros((nregion+1,),dtype=np.int64) - self._offset[1:] = np.cumsum(np.bincount(regionid[idx],minlength=nregion)) - self._offsetbd = np.full((nregion+1,),self._offset[-1],dtype=np.int64) - self._offsetbd[1:] += np.cumsum(np.bincount(regionid[idxbd],minlength=nregion)) + self._offset = np.zeros((2*nfeatures+1,),dtype=np.int64) + self._offset[1:nfeatures+1] = np.cumsum(np.bincount(regionid[idx],minlength=nfeatures)) + self._offset[nfeatures+1:] = np.cumsum(np.bincount(regionid[idxbd],minlength=nfeatures)) + self._offset[nfeatures+1:] += self._offset[nfeatures] self._faces = np.concatenate(( faces[idx,:][np.argsort(regionid[idx]),:], faces[idxbd,:][np.argsort(regionid[idxbd]),:]),axis=0) - # # - # # self._gradient = (contour.point_arrays['Gradients'][self._faces[:,1],:] - - # # Compute the volume and area per cell. For the volume computation, an arbitrary component - # # 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 = self._points[self._faces[:,1],:] - # B = self._points[self._faces[:,2],:] - # C = self._points[self._faces[:,3],:] - # cn = np.cross(B-A,C-A) - # # Check if cell normal points in direction of gradient. If not, switch vertex order. - # idx = (contour.point_arrays['Gradients'][self._faces[:,1],:]*cn).sum(axis=-1)>0 - # # print(idx.shape,np.sum(idx),self._faces.shape,self._faces[idx,2:].shape,self._faces[idx,3:1:-1].shape) - # # self._faces[np.ix_(idx,[2,3])] = self._faces[np.ix_(idx,[3,2])] - # # cn[idx,:] = -cn[idx,:] - # # Compute area and signed volume per cell - # cc = (A+B+C)/3 - # self._cell_areas = 0.5*np.sqrt(np.square(cn).sum(axis=1)) - # self._cell_volumes = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component] + gradients = np.concatenate(( + gradients[idx,:][np.argsort(regionid[idx]),:], + gradients[idxbd,:][np.argsort(regionid[idxbd]),:]),axis=0) + # Compute the volume and area per cell. + if report: print('[Features3d.triangulate] calculating surface area and volume...') + A = self._points[self._faces[:,1],:] + B = self._points[self._faces[:,2],:] + C = self._points[self._faces[:,3],:] + cn = np.cross(B-A,C-A) + # Check if cell normal points in direction of gradient. If not, switch vertex order. + idx = (gradients*cn).sum(axis=-1)>0 + self._faces[np.ix_(idx,[2,3])] = self._faces[np.ix_(idx,[3,2])] + cn[idx,:] = -cn[idx,:] + # Compute area and signed volume per cell + cc = (A+B+C)/3 + cell_areas = 0.5*np.sqrt(np.square(cn).sum(axis=1)) + cell_volumes = 0.5*np.mean(cn*cc,axis=1) + @jit(nopython=True) + def __sum_up_cell_data(cell_data,offset,incl_boundary): + nfeat = (len(offset)-1)//2 + niter = 2*nfeat if incl_boundary else nfeat + output = np.zeros(nfeat,dtype=cell_data.dtype) + for ifeat in range(niter): + for ii in range(offset[ifeat],offset[ifeat+1]): + output[ifeat%nfeat] += cell_data[ii] + return output + self._areas = __sum_up_cell_data(cell_areas, self._offset,False) + self._volumes = __sum_up_cell_data(cell_volumes,self._offset,True ) return @property @@ -832,23 +841,19 @@ class Features3d: def nfeatures(self): return self._nfeatures @property - def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1]) + def areas(self): return self._areas @property - def cell_volumes(self): return np.split(self._cell_volumes,self._offset[1:-1]) - - def fill_holes(self,report=False): - '''Remove triangulation which is fully enclosed by another. These regions - will have a negative volume due to the direction of their normal vector.''' - self.discard_features(np.flatnonzero(self.volumes()<0),report=report) - return + def volumes(self): return self._volumes def reduce_noise(self,threshold=1e-5,is_relative=True,report=False): '''Discards all objects with smaller volume than a threshold.''' if is_relative: vol_domain = np.prod(self.spacing*self.dimensions) threshold = threshold*vol_domain - self.discard_features(np.flatnonzero(np.abs(self.volumes())0: + self.discard_features(features,report=report) return def discard_features(self,features,clean_points=False,report=False): @@ -867,9 +872,9 @@ class Features3d: # Save former number of faces for report nfaces = self._faces.shape[0] # Delete indexed elements from arrays - self._faces = np.delete(self._faces,idx,axis=0) - self._cell_areas = np.delete(self._cell_areas,idx,axis=0) - self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0) + self._faces = np.delete(self._faces,idx,axis=0) + self._areas = np.delete(self._areas,features,axis=0) + self._volumes = np.delete(self._volumes,features,axis=0) # Correct offset self._offset[1:] = self._offset[1:]-np.cumsum(gapsize) self._offset = np.delete(self._offset,features,axis=0) @@ -897,52 +902,43 @@ class Features3d: '''Returns the surface area of feature. If feature is None, total surface area of all features is returned.''' if feature is None: - return np.sum(self._cell_areas) + return np.sum(self.areas) else: - return np.sum(self._cell_areas[self._offset[feature:feature+2]]) - - def areas(self): - '''Returns an array with surface areas of all features.''' - return np.add.reduceat(self._cell_areas,self._offset[:-1]) + return self.areas[feature] def volume(self,feature): '''Returns volume enclosed by feature. If feature isNone, total volume of all features is returned.''' if feature is None: - return np.sum(self._cell_volumes) + return np.sum(self.volumes) else: - return np.sum(self._cell_volumes[self._offset[feature:feature+2]]) - - def volumes(self): - '''Returns an array with volumes of all features.''' - return np.add.reduceat(self._cell_volumes,self._offset[:-1]) - - + return self.volumes[feature] 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._points[self._faces[:,1:],1],axis=1) - max1 = np.amax(self._points[self._faces[:,1:],1],axis=1) - min2 = np.amin(self._points[self._faces[:,1:],2],axis=1) - max2 = np.amax(self._points[self._faces[:,1:],2],axis=1) - elif kdaxis==1: - min1 = np.amin(self._points[self._faces[:,1:],0],axis=1) - max1 = np.amax(self._points[self._faces[:,1:],0],axis=1) - min2 = np.amin(self._points[self._faces[:,1:],2],axis=1) - max2 = np.amax(self._points[self._faces[:,1:],2],axis=1) - elif kdaxis==2: - min1 = np.amin(self._points[self._faces[:,1:],0],axis=1) - max1 = np.amax(self._points[self._faces[:,1:],0],axis=1) - min2 = np.amin(self._points[self._faces[:,1:],1],axis=1) - max2 = np.amax(self._points[self._faces[:,1:],1],axis=1) - else: - raise ValueError("Invalid kdaxis.") - 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) + from scipy.spatial import KDTree + from numba import jit + @jit(nopython=True,cache=True) + def __get_center_and_search_radius(points,faces,kdaxis): + nfaces = faces.shape[0] + center = np.zeros((nfaces,2),dtype=points.dtype) + radius = 0.0 + if kdaxis==0: i1,i2 = 1,2 + elif kdaxis==1: i1,i2 = 0,2 + elif kdaxis==2: i1,i2 = 0,1 + for iface in range(nfaces): + min1 = min(min(points[faces[iface,1],i1],points[faces[iface,2],i1]),points[faces[iface,3],i1]) + max1 = max(max(points[faces[iface,1],i1],points[faces[iface,2],i1]),points[faces[iface,3],i1]) + min2 = min(min(points[faces[iface,1],i2],points[faces[iface,2],i2]),points[faces[iface,3],i2]) + max2 = max(max(points[faces[iface,1],i2],points[faces[iface,2],i2]),points[faces[iface,3],i2]) + center[iface,0] = 0.5*(max1+min1) + center[iface,1] = 0.5*(max2+min2) + radius = max(radius,(max1-min1)**2+(max2-min2)**2) + radius = 0.5*np.sqrt(radius) + return center,radius + center,radius = __get_center_and_search_radius(self._points,self._faces,kdaxis) + self._kdtree = KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes, + copy_data=False,balanced_tree=balanced_tree) self._kdaxis = kdaxis self._kdradius = radius if report: @@ -968,6 +964,9 @@ 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] ''' + from numba import jit + from time import time + coords = np.array(coords) assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as Nx3 array." @@ -977,44 +976,68 @@ class Features3d: elif self._kdaxis==1: query_axis = [0,2] elif self._kdaxis==2: query_axis = [0,1] - cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius) + t__ = time() + candidates = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius).tolist() + print('query',time()-t__) + + t__ = time() + cand_num = np.asarray(tuple(map(len,candidates))) + cand_arr = np.concatenate(candidates).astype(np.int) + print('remap',time()-t__) - Ncoord = coords.shape[0] raydir = np.zeros((3,),dtype=self._points.dtype) raydir[self._kdaxis] = 1.0 - in_feat = np.empty((Ncoord,),dtype=np.int) - is_hit = np.empty((Ncoord,),dtype=np.int) - hit_dir_ = np.empty((Ncoord,),dtype=np.int) - face_ = np.empty((Ncoord,),dtype=np.int) - print('point 12156 = ',coords[12156]) - - for ii in range(Ncoord): - hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir, - self._points[self._faces[cand[ii],1],:], - self._points[self._faces[cand[ii],2],:], - self._points[self._faces[cand[ii],3],:]) - # if ii==12279: print('DEBUG',hit_idx,t,N) - if ii==12156: print('DEBUG',hit_idx,t,hit_dir) - if hit_idx is None: - in_feat[ii] = -1 - else: - idx = np.argmin(np.abs(t)) - if ii==12156: - for kk in hit_idx: - print(cand[ii][kk]) - # print(self._points[self._faces[cand[ii][hit_idx[idx]],1:],:]) - # print(cand[ii][hit_idx[idx]]) - # ifeat = self.feature_from_face(cand[ii][hit_idx[idx]]) - # print(self._offset[ifeat],self._offset[ifeat+1]) - # stop - # if hit_dir[idx]>0: - # in_feat[ii] = self.feature_from_face(cand[ii][hit_idx[idx]]) - # else: - # in_feat[ii] = -1 + # + @jit(nopython=True,error_model='numpy',cache=True) + def __ray_triangle_intersect(R,dR,A,B,C): + '''Implements the Möller–Trumbore ray-triangle intersection algorithm. Taken from + https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d''' + E1 = B-A + E2 = C-A + N = np.cross(E1,E2) + norm = np.sqrt(np.sum(N*N)) + det = -np.sum(dR*N) + invdet = 1.0/det + AO = R-A + DAO = np.cross(AO,dR) + u = np.sum(E2*DAO)*invdet + v = -np.sum(E1*DAO)*invdet + t = np.sum(AO*N) *invdet + return (np.abs(det)>=1e-7*norm and t>=0.0 and u>=0 and v>=0 and (u+v)-1.0<=0), t + # @jit(nopython=True,cache=True) + def __get_nearest_face(coords,raydir,cand_arr,cand_num,faces,points): + N = coords.shape[0] + nearface = np.zeros(N,dtype=np.int64) + isinside = np.zeros(N,dtype=np.bool_) + jj = 0 + for ii in range(N): + t = np.inf + f = -1 + n = 0 + for kk in range(cand_num[ii]): + idx = cand_arr[jj] + hit_,t_ = __ray_triangle_intersect(coords[ii],raydir, + points[faces[idx,1],:], + points[faces[idx,2],:], + points[faces[idx,3],:]) + if hit_: + n += 1 + if t_=0),Ncoord)) - return in_feat + print('{} of {} points are located inside of features.'.format(sum(isinside),coords.shape[0])) + return output def feature_from_face(self,idx_cell): '''Gets feature ID for a given cell.''' @@ -1024,30 +1047,20 @@ class Features3d: # as soon as idx_cell is equal or larger than offset, the argmax function # short circuits and returns the index of the first occurence without # checking any value afterwards. - return np.searchsorted(self._offset,idx_cell,side='right')-1 + return (np.searchsorted(self._offset,idx_cell,side='right')-1)%self._nfeatures def faces_from_feature(self,features,incl_boundary=False): '''Returns indices of cells which belong to given features.''' - from collections.abc import Iterable + features = self.list_of_features(features) faces = [] - if features is None: - faces.append(self._faces) - if incl_boundary: - faces.append(self._faces) - elif not isinstance(features,Iterable): - idx = slice(self._offset[features],self._offset[features+1]) + for feature in features: + idx = np.arange(self._offset[feature],self._offset[feature+1]) faces.append(self._faces[idx,:]) - if incl_boundary: - idx = slice(self._offsetbd[features],self._offsetbd[features+1]) - faces.append(self._faces[idx,:]) - else: + if incl_boundary: for feature in features: - idx = np.arange(self._offset[feature],self._offset[feature+1]) + idx = np.arange(self._offset[self._nfeatures+feature], + self._offset[self._nfeatures+feature+1]) faces.append(self._faces[idx,:]) - if incl_boundary: - for feature in features: - idx = np.arange(self._offsetbd[feature],self._offsetbd[feature+1]) - faces.append(self._faces[idx,:]) return np.concatenate(faces,axis=0) def regionid_from_feature(self,features,incl_boundary=False): @@ -1059,18 +1072,23 @@ class Features3d: regionid.append(np.full(nfaces,feature,dtype=np.int64)) if incl_boundary: for feature in features: - nfaces = self._offsetbd[feature+1]-self._offsetbd[feature] + nfaces = (self._offset[self._nfeatures+feature+1]- + self._offset[self._nfeatures+feature]) regionid.append(np.full(nfaces,feature,dtype=np.int64)) return np.concatenate(regionid,axis=0) - def nfaces_of_feature(self,features): + def nfaces_of_feature(self,features,incl_boundary=False): '''Returns number of cells which belong to given features. Can be used to construct new offsets.''' features = self.list_of_features(features) - ncells = [] + nfaces = [] for feature in features: - ncells.append(self._offset[feature+1]-self._offset[feature]) - return np.array(ncells) + nfaces.append(self._offset[feature+1]-self._offset[feature]) + if incl_boundary: + for feature in features: + nfaces.append(self._offset[self._nfeatures+feature+1]- + self._offset[self._nfeatures+feature]) + return np.concatenate(nfaces,axis=0) def list_of_features(self,features): '''Ensures that 'features' is a list.''' @@ -1082,84 +1100,6 @@ class Features3d: else: return list(features) - # @staticmethod - # def ray_triangle_intersection(r0,dr,v0,v1,v2): - # '''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. - # ''' - # v0v1 = v1-v0 - # v0v2 = v2-v0 - # pvec = np.cross(dr,v0v2) - # det = (v0v1*pvec).sum(axis=-1) # det<0: from behind, det>0 from front ("culling") - # norm = 1e0 # we should normalize the unit vector? - # det[np.abs(det)<1e-6*norm] = np.nan # mask to avoid numpy runtime errors - # invDet = 1./det - # tvec = r0-v0 - # qvec = np.cross(tvec,v0v1) - # u = (tvec*pvec).sum(axis=-1)*invDet - # v = (dr*qvec).sum(axis=-1)*invDet - # is_hit = np.logical_and( - # np.logical_and( - # np.logical_and( - # np.isfinite(det),u>=-1e-6), - # v>=1e-6), - # (u+v)-1.0<=1e-6) - # hit_idx = np.flatnonzero(is_hit) - # if len(hit_idx)==0: return (None,None,None) - # t = (v0v2[hit_idx,:]*qvec[hit_idx,:]).sum(axis=-1)*invDet[hit_idx]; - # return hit_idx,t,np.sign(det[hit_idx]) - @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>=-1e-6), - v>=-1e-6), - (u+v)-1.0<=1e-6) - 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,np.sign(t)*np.sign(det[hit_idx]) - return hit_idx,t,N[hit_idx] - def cmap_features(self,nfeatures,name='tab10'): from matplotlib.colors import ListedColormap if name=='tab10': # seaborn/matplotlib tab10 @@ -1185,10 +1125,8 @@ class Features3d: def to_vtk(self,features,incl_regionid=False,incl_boundary=False): import pyvista as pv - faces = self.faces_from_feature(features,incl_boundary=incl_boundary) - print(faces.shape) + faces = self.faces_from_feature(features,incl_boundary=incl_boundary) output = pv.PolyData(self._points,faces) - print(output.n_faces) if incl_regionid: output.cell_arrays['RegionId'] = self.regionid_from_feature(features,incl_boundary=incl_boundary) return output