From d95453ad45408c85cdc14dba27f29187efc493bb Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 13 Aug 2021 22:05:35 +0200 Subject: [PATCH] WIP: completely messed up. trying to fix surface normal directions --- field.py | 256 +++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 213 insertions(+), 43 deletions(-) diff --git a/field.py b/field.py index 024cc5d..2d42186 100644 --- a/field.py +++ b/field.py @@ -659,7 +659,7 @@ class Features3d: # Assign basic properties to class variables self.origin = np.array(origin,dtype=np.float) self.spacing = np.array(spacing,dtype=np.float) - self.dimensions = input.shape + self.dimensions = np.array(input.shape,dtype=np.int) self.periodicity = tuple(bool(x) for x in periodicity) # If regions are supposed to be inverted, i.e. the interior consists of values # smaller than the threshold instead of larger, change the sign of the array. @@ -712,8 +712,14 @@ class Features3d: # Compute full contour: ensure we only have triangles to efficiently extract points # later on. if report: print('[Features3d.triangulate] computing isocontour using {}...'.format(contour_method)) - contour = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False) + 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...') @@ -769,6 +775,8 @@ class Features3d: # 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. @@ -780,8 +788,7 @@ class Features3d: 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 = np.ascontiguousarray( - contour.points[contour.faces.reshape(contour.n_faces,4)[ind,1:]]) + self._vertices = contour.points[contour.faces.reshape(contour.n_faces,4)[ind,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'. @@ -791,8 +798,39 @@ class Features3d: C = self._vertices[:,2,:] cn = np.cross(B-A,C-A) 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 @@ -821,19 +859,17 @@ 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[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, + 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) + 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,:]) + # if ii==12279: print('DEBUG',hit_idx,t,N) 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)) - if t[idx]*(N[idx,:]*raydir).sum(axis=-1)>0: - in_feat[ii] = self.feature_by_cell(cand[ii][idx]) + # 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 - return in_feat + is_hit[ii] = 0 + 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_ - def feature_by_cell(self,idx_cell): + def feature_from_cellidx(self,idx_cell): '''Gets feature ID for a given cell.''' - return np.argmax(idx_cell>self._offset) + # This is an optimized solution for monotonically increasing arrays: + # as long as offset is smaller than the cell index, the boolean operation + # returns False. As soon as the first True is returned, i.e. + # 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) + def cellidx_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 + elif not isinstance(features,Iterable): + idx = slice(self._offset[features],self._offset[features+1]) + elif concat: + idx = [] + for feature in features: + rng_ = np.arange(self._offset[feature],self._offset[feature+1]) + idx.append(rng_) + idx = np.concatenate(idx,axis=0) + else: + idx = [] + for feature in features: + rng_ = slice(self._offset[feature],self._offset[feature+1]) + idx.append(rng_) + return idx + + def ncells_from_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) + ncells = [] + for feature in features: + ncells.append(self._offset[feature+1]-self._offset[feature]) + return np.array(ncells) + + def list_of_features(self,features): + '''Ensures that 'features' is a list.''' + from collections.abc import Iterable + if features is None: + return np.arange(0,self._nfeatures) + if not isinstance(features,Iterable): + return [features] + 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 @@ -988,46 +1129,75 @@ class Features3d: 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.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,N[hit_idx] + return hit_idx,t,np.sign(t)*np.sign(det[hit_idx]) - def faces_points_representation(self,labels,reduce_points=False): + def faces_points_representation(self,features,reduce_points=False): from collections.abc import Iterable # Select relevant cells - if labels is None: - idx = None - else: - if not isinstance(labels,Iterable): labels = [labels] - idx = [] - for label in labels: - rng_ = np.arange(self._offset[label],self._offset[label+1]) - idx.append(rng_) - idx = np.concatenate(idx,axis=0) + 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 - def to_vtk(self,labels,reduce_points=False): + def cmap_features(self,nfeatures,name='tab10'): + from matplotlib.colors import ListedColormap + if name=='tab10': # seaborn/matplotlib tab10 + cmap = [(0.12156862745098039, 0.4666666666666667, 0.7058823529411765), + (1.0, 0.4980392156862745, 0.054901960784313725), + (0.17254901960784313, 0.6274509803921569, 0.17254901960784313), + (0.8392156862745098, 0.15294117647058825, 0.1568627450980392), + (0.5803921568627451, 0.403921568627451, 0.7411764705882353), + (0.5490196078431373, 0.33725490196078434, 0.29411764705882354), + (0.8901960784313725, 0.4666666666666667, 0.7607843137254902), + (0.4980392156862745, 0.4980392156862745, 0.4980392156862745), + (0.7372549019607844, 0.7411764705882353, 0.13333333333333333), + (0.09019607843137255, 0.7450980392156863, 0.8117647058823529)] + else: + raise NotImplementedError('Colormap not implemented.') + ncolor = len(cmap) + icolor = 0 + output = np.empty((nfeatures,3)) + for ii in range(nfeatures): + output[ii] = cmap[icolor] + 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).''' import pyvista as pv - from collections.abc import Iterable - faces,points = self.faces_points_representation(labels,reduce_points=reduce_points) - return pv.PolyData(points,faces) + 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 + return output class BinaryFieldNd: def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False):