From c5908e00f5ad23473da1bf27340c22fea9ec1e0f Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Mon, 16 Aug 2021 12:58:23 +0200 Subject: [PATCH] update: inside/outside detection does not work properly --- field.py | 248 +++++++++++++++++++++++-------------------------------- 1 file changed, 105 insertions(+), 143 deletions(-) diff --git a/field.py b/field.py index 2d42186..4c261c0 100644 --- a/field.py +++ b/field.py @@ -714,12 +714,6 @@ 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." - - # print(contour.point_arrays['Gradients'].shape) - # print('gradient',time()-t) - # print(np.sum(gn>0),'of',len(gn)) - - # Compute the connectivity of the triangulated surface: first we run an ordinary # connectivity filter neglecting periodic wrapping. if report: print('[Features3d.triangulate] computing connectivity...') @@ -767,70 +761,38 @@ class Features3d: while target_ != map_[target_]: # map it recursively target_ = map_[target_] map_[source_] = target_ - # - map_ = np.unique(map_,return_inverse=True)[1] + map_ = np.unique(map_,return_inverse=True)[1] point_labels = map_[point_labels] self._nfeatures = np.max(map_)+1 # starts with zero # Labels are now stored as point data. To efficiently convert it to cell data, the first # point of each cell determines the value for this cell. - if report: print('[Features3d.triangulate] identifying cell based labels...') - cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]] - cell_gradient = contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1]] - print(cell_gradient.shape) - # We are done with VTK for now. Since we are only dealing with triangles, let us - # store them in a structure which is quicker to access and already sorted based - # on labels. - # _vertices has dim(ncells,nvert,ndim), where nvert=3 (three vertices per triangle) - # and ndim=3 (three dimensional coordinates) - # _offset is the cumulative number of cells for each feature and the offset to - # access the _vertices array by feature. - if report: print('[Features3d.triangulate] storing vertices...') - ind = np.argsort(cell_labels) - self._offset = np.zeros((self._nfeatures+1,),dtype=np.int64) - self._offset[1:] = np.cumsum(np.bincount(cell_labels)) - self._vertices = contour.points[contour.faces.reshape(contour.n_faces,4)[ind,1:]] + if report: print('[Features3d.triangulate] sorting faces by labels...') + cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]] + offset = np.zeros((self._nfeatures+1,),dtype=np.int64) + offset[1:] = np.cumsum(np.bincount(cell_labels)) + ind = np.argsort(cell_labels) + self._offset = offset + self._faces = contour.faces.reshape(contour.n_faces,4)[ind,:] + self._points = contour.points + # self._offset = offset + # self._faces = contour.faces.reshape(contour.n_faces,4) + # self._points = contour.points # 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._vertices[:,0,:] - B = self._vertices[:,1,:] - C = self._vertices[:,2,:] + 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 + 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 - # print(contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]].shape) - cg = (contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,0,:]+ - contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,1,:]+ - contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,2,:])/3 - ng = (cg*cn).sum(axis=-1) - print('interpolated:') - print(cg) - print('component 0:') - print(contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1:]][:,0,:]) - print('projection') - print(sum(ng>0),len(ng)) - - # self.cg = cg - # self.cc = cc - self.cg = contour.point_arrays['Gradients'] - self.cc = contour.points - - # .sum(axis=1))/3)*cn).sum(axis=-1) 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] - - # # Compute gradient of field at cell center to get correct cell normal direction - # t = time() - # cn = cn/np.sqrt(np.square(cn).sum()) - # ccpix = ((cc-self.origin)/self.spacing).transpose() - # ccnpix = (((cc+1e-4*cn)-self.origin)/self.spacing).transpose() - # val_cc = ndimage.map_coordinates(self._input,ccpix,order=1,prefilter=False) - # val_ccn = ndimage.map_coordinates(self._input,ccnpix,order=1,prefilter=False) - # gn = (val_ccn-val_cc) - # print('gradient',time()-t) - # print(cg) - # print(np.sum(cg>0),'of',cg.shape) - # stop return @property @@ -859,7 +821,7 @@ class Features3d: self.discard_features(np.flatnonzero(np.abs(self.volumes())12279: - print('Point coordinates:',coords[12279,:]) - Ncoord = coords.shape[0] - raydir = np.zeros((3,),dtype=self._vertices.dtype) + 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._vertices[cand[ii],0,:], - self._vertices[cand[ii],1,:], - self._vertices[cand[ii],2,:]) + 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 - is_hit[ii] = -1 - hit_dir_[ii]= 0 - face_[ii] = -1 else: idx = np.argmin(np.abs(t)) - # print(cand[ii][hit_idx[idx]]) - if ii==12279: - # print(cand[ii][idx],N[idx,:],N[idx,:]/np.sqrt(np.square(N[idx,:]).sum()),(N[idx,:]*raydir).sum(axis=-1),t[idx]*(N[idx,:]*raydir).sum(axis=-1)) - print(t[idx],hit_dir[idx],raydir) - print(self._vertices[cand[ii][hit_idx[idx]],:,:]) - hit_dir_[ii]= hit_dir[idx] - face_[ii] = cand[ii][hit_idx[idx]] - if hit_dir[idx]<0: - in_feat[ii] = self.feature_from_cellidx(cand[ii][hit_idx[idx]]) - is_hit[ii] = 1 - else: - in_feat[ii] = -1 - is_hit[ii] = 0 + 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 if report: print('[Features3d.inside_feature]',end=' ') print('{} of {} points are located inside of features.'.format(sum(in_feat>=0),Ncoord)) - return in_feat #,is_hit,hit_dir_,face_ + return in_feat - def feature_from_cellidx(self,idx_cell): + def feature_from_face(self,idx_cell): '''Gets feature ID for a given cell.''' # This is an optimized solution for monotonically increasing arrays: # as long as offset is smaller than the cell index, the boolean operation @@ -1020,13 +991,13 @@ 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.argmax(self._offset>=idx_cell) + return np.searchsorted(self._offset,idx_cell,side='right')-1 - def cellidx_from_feature(self,features,concat=True): + def faces_from_feature(self,features,concat=True): '''Returns indices of cells which belong to given features.''' from collections.abc import Iterable if features is None: - idx = None + idx = slice(None) elif not isinstance(features,Iterable): idx = slice(self._offset[features],self._offset[features+1]) elif concat: @@ -1034,7 +1005,10 @@ class Features3d: for feature in features: rng_ = np.arange(self._offset[feature],self._offset[feature+1]) idx.append(rng_) - idx = np.concatenate(idx,axis=0) + if len(idx)==1: + idx = np.array(idx) + else: + idx = np.concatenate(idx,axis=0) else: idx = [] for feature in features: @@ -1042,7 +1016,7 @@ class Features3d: idx.append(rng_) return idx - def ncells_from_feature(self,features): + def nfaces_of_feature(self,features): '''Returns number of cells which belong to given features. Can be used to construct new offsets.''' features = self.list_of_features(features) @@ -1136,27 +1110,8 @@ class Features3d: 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]) - - def faces_points_representation(self,features,reduce_points=False): - from collections.abc import Iterable - # Select relevant cells - idx = self.cellidx_from_feature(features) - if reduce_points: - points,faces_ = np.unique(self._vertices[idx,:,:].reshape(-1,3),return_inverse=True,axis=0) - ncells = faces_.shape[0]//3 - faces = np.empty((ncells,4),dtype=np.int) - faces[:,0] = 3 - faces[:,1:] = faces_.reshape(ncells,3) - print('reduced:',faces.shape) - else: - points = self._vertices[idx,:,:].reshape((-1,3)) - ncells = points.shape[0]//3 - faces = np.empty((ncells,4),dtype=np.int) - faces[:,0] = 3 - faces[:,1:] = np.arange(0,3*ncells).reshape(ncells,3) - print('full:',faces.shape) - return faces,points + # 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 @@ -1181,24 +1136,31 @@ class Features3d: icolor = (icolor+1)%ncolor return ListedColormap(output) - def to_vtk(self,features,reduce_points=False,store_labels=False,remap_labels=False): - '''Returns a VTK/pyvista object for the isocontour of a single feature, - list of features, or the full isocontour (when None is passed).''' + def to_vtk(self,features): import pyvista as pv - faces,points = self.faces_points_representation(features,reduce_points=reduce_points) - output = pv.PolyData(points,faces) - if store_labels: - output.cell_arrays['label'] = np.empty(faces.shape[0],dtype=np.int) - offset = 0 - for ii,feature in enumerate(self.list_of_features(features)): - ncells = self.ncells_from_feature(feature) - if remap_labels: - output.cell_arrays['label'][offset:offset+ncells] = ii - else: - output.cell_arrays['label'][offset:offset+ncells] = feature - offset += ncells + idx = self.faces_from_feature(features) + return pv.PolyData(self._points,self._faces[idx,:]) + + def to_vtk2(self,features): + import pyvista as pv + idx = self.faces_from_feature(features) + faces_ = self._faces[idx,:] + A = self._points[faces_[:,1],:] + B = self._points[faces_[:,2],:] + C = self._points[faces_[:,3],:] + cn = np.cross(B-A,C-A) + cn = cn/np.sqrt(np.square(cn).sum(axis=-1,keepdims=True)) + cc = (A+B+C)/3 + # print(A[0,:],B[0,:],C[0,:],cn[0,:]) + # print(A.shape,B.shape,C.shape,cn.shape,faces_.shape,self._faces.shape) + output = pv.PolyData(cc) + output['cellnormal'] = cn return output + def vtk_faces(self,face_idx): + import pyvista as pv + return pv.PolyData(self._points,self._faces[face_idx,:]) + class BinaryFieldNd: def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False): assert isinstance(input,np.ndarray) and input.dtype==bool,\