its getting better
This commit is contained in:
parent
16e67ad666
commit
e95b7dc74a
317
field.py
317
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]+tolerance[jj] and
|
||||
points[faces[ii,2],jj]<bounds[2*jj]+tolerance[jj] and
|
||||
points[faces[ii,3],jj]<bounds[2*jj]+tolerance[jj]):
|
||||
output[ii] = 2*jj
|
||||
elif (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]<search_dist:
|
||||
output[jj,:] = (3,face_uni_hi[ii],face_uni_lo[idx[ii]],face_uni_hi[ii])
|
||||
jj += 1
|
||||
return output[:jj,:]
|
||||
#####
|
||||
tol_overlap = 1e-5*self.spacing
|
||||
bd_face_tag = __get_boundary_faces(faces,points,tuple(contour.bounds),tol_overlap)
|
||||
#####
|
||||
faces_connect = [faces]
|
||||
for axis in range(3):
|
||||
# Build a nearest-neighbor search KD-trees for boundary points so that we can connect
|
||||
# them to the boundary faces when needed
|
||||
# t__ = time()
|
||||
pos_bd_lo = self.origin[axis]
|
||||
pos_bd_hi = self.origin[axis]+self.spacing[axis]*(self.dimensions[axis]-1)
|
||||
search_dist = 1e-5*self.spacing[axis]
|
||||
idx_lo = np.flatnonzero(np.abs(contour.points[:,axis]-pos_bd_lo)<search_dist)
|
||||
idx_hi = np.flatnonzero(np.abs(contour.points[:,axis]-pos_bd_lo)<search_dist)
|
||||
if not self.periodicity[axis]: continue
|
||||
sl_pln = [0,1,2]
|
||||
del sl_pln[axis]
|
||||
# print(len(idx_lo),len(idx_hi))
|
||||
# kd_lo = KDTree(contour.points[np.ix_(idx_lo,sl_pln)],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
|
||||
# kd_hi = KDTree(contour.points[np.ix_(idx_hi,sl_pln)],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
|
||||
# print('KD tree build:',time()-t__)
|
||||
# Compute the contour on the boundary: the normal should point outwards.
|
||||
sl_lo = 3*[slice(None)]
|
||||
sl_lo[axis] = 0
|
||||
sl_hi = 3*[slice(None)]
|
||||
sl_hi[axis] = -1
|
||||
origin = self.origin.copy()
|
||||
origin[axis] -= self.spacing[axis]
|
||||
dimensions = self.dimensions.copy()
|
||||
dimensions[axis] = 2
|
||||
#
|
||||
tmp = np.empty(dimensions,dtype=self._input.dtype,order='F')
|
||||
tmp[tuple(sl_hi)] = self._input[tuple(sl_lo)]
|
||||
tmp[tuple(sl_lo)] = -1e-30
|
||||
print(origin)
|
||||
print(self.spacing)
|
||||
print(dimensions)
|
||||
#
|
||||
planevtk = pv.UniformGrid()
|
||||
planevtk.dimensions = dimensions
|
||||
planevtk.spacing = self.spacing
|
||||
planevtk.origin = origin
|
||||
planevtk.point_arrays['data'] = tmp.ravel('F')
|
||||
# Contour for lower boundary
|
||||
contour_bd = planevtk.contour([self._threshold],method=contour_method)
|
||||
faces_bd = contour_bd.faces.reshape(-1,4).copy()
|
||||
points_bd = contour_bd.points.copy()
|
||||
print(points_bd)
|
||||
# Find points which connect to main contour and update faces
|
||||
kd = KDTree(points_bd[:,sl_pln],leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
|
||||
kddist = 1e-1*self.spacing[axis]
|
||||
print(kddist)
|
||||
# ptidx = kd.query(contour.points[np.ix_(idx_lo,sl_pln)],k=1)#,distance_upper_bound=kddist)
|
||||
ptidx = kd.query(contour.points[np.ix_(idx_lo,sl_pln)],k=1)#,distance_upper_bound=kddist)
|
||||
# print(ptidx[0])
|
||||
# print(np.min(contour.points[idx_lo,0]),np.max(contour.points[idx_lo,0]))
|
||||
print(np.min(points_bd[:,0]),np.max(points_bd[:,0]))
|
||||
print('connections:',np.sum(ptidx[0]<kddist))
|
||||
print('boundary points:',points_bd[:,sl_pln].shape)
|
||||
print('selected points:',contour.points[np.ix_(idx_lo,sl_pln)].shape)
|
||||
stop
|
||||
#
|
||||
if self.periodicity[axis]:
|
||||
sl_swap = [1,2,3]
|
||||
del sl_swap[axis]
|
||||
self._faces_bd += [contour_bd.faces.reshape(-1,4).copy()]
|
||||
self._faces_bd[-1][:,sl_swap] = self._faces_bd[-1][:,sl_swap[::-1]]
|
||||
self._points_bd += [contour_bd.points.copy()]
|
||||
self._points_bd[-1][axis] += self.spacing[axis]*(self.dimensions[axis]-1)
|
||||
else:
|
||||
origin[axis] = self.origin[axis]+self.spacing[axis]*(self.dimensions[axis]-1)
|
||||
tmp[tuple(sl_lo)] = self._input[tuple(sl_hi)]
|
||||
tmp[tuple(sl_hi)] = -1e-30
|
||||
planevtk.origin = origin
|
||||
planevtk.point_arrays['data'] = tmp.ravel('F')
|
||||
contour_bd = planevtk.contour([self._threshold],method=contour_method)
|
||||
self._faces_bd += [contour_bd.faces.reshape(-1,4).copy()]
|
||||
self._points_bd += [contour_bd.points.copy()]
|
||||
print('boundary contour:',time()-t__)
|
||||
#
|
||||
# Compute the connectivity of the triangulated surface: first we run an ordinary
|
||||
# connectivity filter neglecting periodic wrapping.
|
||||
if report: print('[Features3d.triangulate] computing connectivity...')
|
||||
alg = vtk.vtkPolyDataConnectivityFilter() # ~twice as fast as default vtkConnectivityFilter
|
||||
face_uni_lo = np.unique(faces[bd_face_tag==2*axis,1:].ravel())
|
||||
face_uni_hi = np.unique(faces[bd_face_tag==2*axis+1,1:].ravel())
|
||||
kdt = KDTree(points[np.ix_(face_uni_lo,sl_pln)],leafsize=10,compact_nodes=False,
|
||||
copy_data=False,balanced_tree=False)
|
||||
dist,idx = kdt.query(points[np.ix_(face_uni_hi,sl_pln)],k=1,
|
||||
distance_upper_bound=tol_overlap[axis])
|
||||
faces_connect += [__connect_faces_periodic(face_uni_lo,face_uni_hi,points,
|
||||
dist,idx,tol_overlap[axis])]
|
||||
#####
|
||||
ncells = faces.shape[0]
|
||||
contour.faces = np.concatenate(faces_connect,axis=0)
|
||||
#####
|
||||
alg = vtk.vtkPolyDataConnectivityFilter()
|
||||
alg.SetInputData(contour)
|
||||
alg.SetExtractionModeToAllRegions()
|
||||
alg.SetColorRegions(True)
|
||||
alg.Update()
|
||||
contour = pv.filters._get_output(alg)
|
||||
# Now determine the boundary points, i.e. points on the boundary which have a corresponding
|
||||
# vertex on the other side of the domain
|
||||
if report: print('[Features3d.triangulate] correcting connectivity for periodic boundaries...')
|
||||
bd_points_xyz = np.zeros((0,3),dtype=contour.points.dtype)
|
||||
bd_points_idx = np.zeros((0,),dtype=np.int)
|
||||
for axis in range(3):
|
||||
if not self.periodicity[axis]: continue
|
||||
# Compute position of boundary, period of domain and set a threshold for "on boundary" condition
|
||||
bound_pos = datavtk.origin[axis]+datavtk.spacing[axis]*(datavtk.dimensions[axis]-1)
|
||||
period = np.zeros((3,))
|
||||
period[axis] = bound_pos-datavtk.origin[axis]
|
||||
bound_dist = 1e-5*datavtk.spacing[axis]
|
||||
# Lower boundary
|
||||
li = np.abs(contour.points[:,axis]-datavtk.origin[axis])<bound_dist
|
||||
bd_points_idx = np.append(bd_points_idx,np.nonzero(li)[0],axis=0)
|
||||
bd_points_xyz = np.append(bd_points_xyz,contour.points[li,:],axis=0)
|
||||
# Upper boundary
|
||||
li = np.abs(contour.points[:,axis]-bound_pos)<bound_dist
|
||||
bd_points_idx = np.append(bd_points_idx,np.nonzero(li)[0],axis=0)
|
||||
bd_points_xyz = np.append(bd_points_xyz,contour.points[li,:]-period,axis=0)
|
||||
# Construct a KD Tree for efficient neighborhood search
|
||||
kd = spatial.KDTree(bd_points_xyz,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
|
||||
# Construct a map to get new labels for regions. We search for pairs of points with
|
||||
# a very small distance inbetween. Then we check if their labels differ and choose
|
||||
# the lower label as the new joint one.
|
||||
point_labels = contour.point_arrays['RegionId']
|
||||
nfeatures = np.max(point_labels)+1
|
||||
map_ = np.array(range(nfeatures),dtype=point_labels.dtype)
|
||||
overlap_dist = 1e-4*np.sqrt(np.square(datavtk.spacing).sum())
|
||||
for (ii,jj) in kd.query_pairs(r=overlap_dist):
|
||||
label_ii = point_labels[bd_points_idx[ii]]
|
||||
label_jj = point_labels[bd_points_idx[jj]]
|
||||
if label_ii!=label_jj:
|
||||
source_ = np.maximum(label_ii,label_jj)
|
||||
target_ = np.minimum(label_ii,label_jj)
|
||||
while target_ != map_[target_]: # map it recursively
|
||||
target_ = map_[target_]
|
||||
map_[source_] = target_
|
||||
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
|
||||
# RegionIds 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.
|
||||
regionid = contour.point_arrays['RegionId'][contour.faces.reshape(-1,4)[:ncells,1]].ravel()
|
||||
nregion = np.amax(regionid)+1
|
||||
# Store the face/vertex data in class arrays. Sort the faces by RegionId for quick access
|
||||
# later on.
|
||||
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
|
||||
# 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]
|
||||
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
|
||||
|
|
|
|||
Loading…
Reference in New Issue