From e95b7dc74af1c485c051c2f4509f55f3eec52ddb Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 18 Aug 2021 02:24:20 +0200 Subject: [PATCH] its getting better --- field.py | 317 +++++++++++++++++++++++++------------------------------ 1 file changed, 144 insertions(+), 173 deletions(-) diff --git a/field.py b/field.py index 8099cce..3c99349 100644 --- a/field.py +++ b/field.py @@ -657,8 +657,8 @@ class Features3d: assert len(periodicity)==3, "'periodicity' must be of length 3" assert isinstance(input,np.ndarray), "'input' must be numpy.ndarray." # Assign basic properties to class variables - self.origin = np.array(origin,dtype=np.float) self.spacing = np.array(spacing,dtype=np.float) + self.origin = np.array(origin,dtype=np.float)-self.spacing 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. @@ -667,13 +667,18 @@ class Features3d: # '_input' is the array which is to be triangulated. Here one trailing ghost cell is # required to get closed surfaces. The data will always be copied since it probably # needs to be copied for vtk anyway (needs FORTRAN contiguous array) and this way - # there will never be side effects on the input data. + # there will never be side effects on the input data. The array is padded with + # a very high value (cannot use NaN or Inf for vtk) in order to close the contour at + # the boundaries. if has_ghost: - self._input = sign_invert*input + self.dimensions = input.shape+np.array((2,2,2)) + self._input = np.full(self.dimensions,-1e30,dtype=input.dtype,order='F') + self._input[1:-1,1:-1,1:-1] = sign_invert*input else: pw = tuple((0,1) if x else (0,0) for x in periodicity) - self._input = np.pad(sign_invert*input,pw,mode='wrap') - self.dimensions = np.array(self._input.shape,dtype=np.int) + self.dimensions = input.shape+np.array((2,2,2))+np.array([sum(x) for x in pw]) + self._input = np.full(self.dimensions,-1e30,dtype=input.dtype,order='F') + 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, @@ -701,6 +706,7 @@ class Features3d: import vtk from scipy import ndimage, spatial from scipy.spatial import KDTree + from numba import jit from time import time # Check if '_input' is available: might have been dropped after initialization assert self._input is not None, "'_input' not available. Initialize object with keep_input=True flag." @@ -710,169 +716,113 @@ class Features3d: datavtk.origin = self.origin datavtk.spacing = self.spacing datavtk.point_arrays['data'] = self._input.ravel('F') - # Compute full contour: ensure we only have triangles to efficiently extract points - # later on. + # Compute full contour: due to the padding, the contour will result in closed object. + # This is necessary for proper inside/outside checks and volume computation. However, + # for surface area computation, the boundary faces have to be removed. Therefore, we + # will isolate them and store them in an extra array. 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." - # Compute contour on boundaries: this is necessary for inside/outside checks and proper - # volume computation for overlapping objects - t__ = time() - self._faces_bd = [] - self._points_bd = [] - offset_pts = contour.n_points - print(offset_pts) + points = contour.points + faces = contour.faces.reshape(-1,4) + # 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. + @jit(nopython=True) + def __get_boundary_faces(faces,points,bounds,tolerance): + assert faces.shape[1]==4 + 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): + for jj in range(3): + if (points[faces[ii,1],jj]bounds[2*jj+1]-tolerance[jj] and + points[faces[ii,2],jj]>bounds[2*jj+1]-tolerance[jj] and + points[faces[ii,3],jj]>bounds[2*jj+1]-tolerance[jj]): + output[ii] = 2*jj+1 + return output + @jit(nopython=True) + def __connect_faces_periodic(face_uni_lo,face_uni_hi,points,dist,idx,search_dist): + N = len(face_uni_hi) + output = np.zeros((N,4),dtype=face_uni_lo.dtype) + jj = 0 + for ii in range(N): + if dist[ii]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] + self._points = points + self._nfeatures = nregion + # + idx = np.flatnonzero(bd_face_tag<0) + idxbd = np.flatnonzero(bd_face_tag>=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._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] return @property @@ -967,6 +917,8 @@ class Features3d: '''Returns an array with volumes of all features.''' return np.add.reduceat(self._cell_volumes,self._offset[:-1]) + + 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 @@ -1074,28 +1026,42 @@ class Features3d: # checking any value afterwards. return np.searchsorted(self._offset,idx_cell,side='right')-1 - def faces_from_feature(self,features,concat=True): + def faces_from_feature(self,features,incl_boundary=False): '''Returns indices of cells which belong to given features.''' from collections.abc import Iterable + faces = [] if features is None: - idx = slice(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]) - elif concat: - idx = [] + faces.append(self._faces[idx,:]) + if incl_boundary: + idx = slice(self._offsetbd[features],self._offsetbd[features+1]) + faces.append(self._faces[idx,:]) + else: for feature in features: - rng_ = np.arange(self._offset[feature],self._offset[feature+1]) - idx.append(rng_) - if len(idx)==1: - idx = np.array(idx) - else: - idx = np.concatenate(idx,axis=0) - else: - idx = [] + idx = np.arange(self._offset[feature],self._offset[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): + '''Returns regionid cell array to be added to PolyData.''' + features = self.list_of_features(features) + regionid = [] + for feature in features: + nfaces = self._offset[feature+1]-self._offset[feature] + regionid.append(np.full(nfaces,feature,dtype=np.int64)) + if incl_boundary: for feature in features: - rng_ = slice(self._offset[feature],self._offset[feature+1]) - idx.append(rng_) - return idx + nfaces = self._offsetbd[feature+1]-self._offsetbd[feature] + regionid.append(np.full(nfaces,feature,dtype=np.int64)) + return np.concatenate(regionid,axis=0) def nfaces_of_feature(self,features): '''Returns number of cells which belong to given features. Can be used @@ -1217,10 +1183,15 @@ class Features3d: icolor = (icolor+1)%ncolor return ListedColormap(output) - def to_vtk(self,features): + def to_vtk(self,features,incl_regionid=False,incl_boundary=False): import pyvista as pv - idx = self.faces_from_feature(features) - return pv.PolyData(self._points,self._faces[idx,:]) + faces = self.faces_from_feature(features,incl_boundary=incl_boundary) + print(faces.shape) + 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 def to_vtk2(self,features): import pyvista as pv