storing triangulation data as long list of vertices now. Seems to have improved performance due to better memory access pattern.
This commit is contained in:
parent
5a88eebf4a
commit
0723c0b769
148
field.py
148
field.py
|
|
@ -696,6 +696,7 @@ class Features3d:
|
|||
import pyvista as pv
|
||||
import vtk
|
||||
from scipy import ndimage, spatial
|
||||
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."
|
||||
# Wrap data for VTK using pyvista
|
||||
|
|
@ -707,22 +708,21 @@ 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)
|
||||
assert contour_.is_all_triangles(), "Contouring produced non-triangle cells."
|
||||
contour = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False)
|
||||
assert contour.is_all_triangles(), "Contouring produced non-triangle cells."
|
||||
# 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
|
||||
alg.SetInputData(contour_)
|
||||
alg.SetInputData(contour)
|
||||
alg.SetExtractionModeToAllRegions()
|
||||
alg.SetColorRegions(True)
|
||||
alg.Update()
|
||||
contour_ = pv.filters._get_output(alg)
|
||||
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...')
|
||||
points = contour_.points
|
||||
bd_points_xyz = np.zeros((0,3),dtype=points.dtype)
|
||||
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
|
||||
|
|
@ -732,19 +732,19 @@ class Features3d:
|
|||
period[axis] = bound_pos-datavtk.origin[axis]
|
||||
bound_dist = 1e-5*datavtk.spacing[axis]
|
||||
# Lower boundary
|
||||
li = np.abs(points[:,axis]-datavtk.origin[axis])<bound_dist
|
||||
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,points[li,:],axis=0)
|
||||
bd_points_xyz = np.append(bd_points_xyz,contour.points[li,:],axis=0)
|
||||
# Upper boundary
|
||||
li = np.abs(points[:,axis]-bound_pos)<bound_dist
|
||||
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,points[li,:]-period,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']
|
||||
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())
|
||||
|
|
@ -760,37 +760,35 @@ class Features3d:
|
|||
#
|
||||
map_ = np.unique(map_,return_inverse=True)[1]
|
||||
point_labels = map_[point_labels]
|
||||
nfeatures = np.max(map_)+1 # starts with zero
|
||||
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...')
|
||||
ncells = contour_.n_faces
|
||||
faces = contour_.faces.reshape(ncells,4)[:,:]
|
||||
cell_labels = point_labels[faces[:,1]]
|
||||
cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]]
|
||||
# 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 = np.ascontiguousarray(
|
||||
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'.
|
||||
if report: print('[Features3d.triangulate] calculating area and volume per cell...')
|
||||
A = points[faces[:,1],:]
|
||||
B = points[faces[:,2],:]
|
||||
C = points[faces[:,3],:]
|
||||
A = self._vertices[:,0,:]
|
||||
B = self._vertices[:,1,:]
|
||||
C = self._vertices[:,2,:]
|
||||
cn = np.cross(B-A,C-A)
|
||||
cc = (A+B+C)/3
|
||||
area = 0.5*np.sqrt(np.square(cn).sum(axis=1))
|
||||
vol = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component]
|
||||
# Now the label is known per cell. We only need to find all cells with the same label
|
||||
# and group them. Therefore we sort the array holding the faces by feature label
|
||||
# and save the offset in another array to easily extract them whenever necessary.
|
||||
if report: print('[Features3d.triangulate] finalizing internal data...')
|
||||
offset = np.zeros((nfeatures+1,),dtype=np.int64)
|
||||
offset[1:] = np.cumsum(np.bincount(cell_labels))
|
||||
ind = np.argsort(cell_labels)
|
||||
self._offset = offset
|
||||
self._faces = faces[ind,:]
|
||||
self._points = points
|
||||
self._cell_areas = area[ind]
|
||||
self._cell_volumes = vol[ind]
|
||||
self._nfeatures = nfeatures
|
||||
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
|
||||
|
|
@ -799,12 +797,6 @@ class Features3d:
|
|||
@property
|
||||
def nfeatures(self): return self._nfeatures
|
||||
|
||||
@property
|
||||
def faces(self): return np.split(self._faces,self._offset)
|
||||
|
||||
@property
|
||||
def points(self): return self._points
|
||||
|
||||
@property
|
||||
def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1])
|
||||
|
||||
|
|
@ -822,7 +814,7 @@ class Features3d:
|
|||
if is_relative:
|
||||
vol_domain = np.prod(self.spacing*self.dimensions)
|
||||
threshold = threshold*vol_domain
|
||||
self.discard_features(np.flatnonzero(self.volumes()<threshold),report=report)
|
||||
self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report)
|
||||
return
|
||||
|
||||
def discard_features(self,labels,clean_points=False,report=False):
|
||||
|
|
@ -840,9 +832,9 @@ class Features3d:
|
|||
gapsize[label] = len(rng_)
|
||||
idx = np.concatenate(idx,axis=0)
|
||||
# Save former number of faces for report
|
||||
nfaces = self._faces.shape[0]
|
||||
ncells = self._vertices.shape[0]
|
||||
# Delete indexed elements from arrays
|
||||
self._faces = np.delete(self._faces,idx,axis=0)
|
||||
self._vertices = np.delete(self._vertices,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)
|
||||
# Correct offset
|
||||
|
|
@ -851,9 +843,7 @@ class Features3d:
|
|||
if report:
|
||||
print('[Features3d.discard_features]',end=' ')
|
||||
print('discarded {:} features with {:} faces.'.format(
|
||||
len(labels),nfaces-self._faces.shape[0]))
|
||||
if clean_points:
|
||||
self.clean_points(report=report)
|
||||
len(labels),ncells-self._vertices.shape[0]))
|
||||
return
|
||||
|
||||
def area(self,feature):
|
||||
|
|
@ -862,7 +852,7 @@ class Features3d:
|
|||
if feature is None:
|
||||
return np.sum(self._cell_areas)
|
||||
else:
|
||||
return np.sum(self.cell_areas[feature-1])
|
||||
return np.sum(self._cell_areas[self._offset[feature:feature+2]])
|
||||
|
||||
def areas(self):
|
||||
'''Returns an array with surface areas of all features.'''
|
||||
|
|
@ -874,7 +864,7 @@ class Features3d:
|
|||
if feature is None:
|
||||
return np.sum(self._cell_volumes)
|
||||
else:
|
||||
return np.sum(self.cell_volumes[feature-1])
|
||||
return np.sum(self._cell_volumes[self._offset[feature:feature+2]])
|
||||
|
||||
def volumes(self):
|
||||
'''Returns an array with volumes of all features.'''
|
||||
|
|
@ -907,7 +897,6 @@ class Features3d:
|
|||
|
||||
# pts = np.random.random((10000,3))
|
||||
|
||||
print(self._points[self._faces[:,1:],0].shape)
|
||||
t = time()
|
||||
|
||||
# xmin = np.amin(self._points[self._faces[:,1:],0],axis=1)
|
||||
|
|
@ -915,8 +904,8 @@ class Features3d:
|
|||
# zmin = np.amin(self._points[self._faces[:,1:],2],axis=1)
|
||||
# zmax = np.amax(self._points[self._faces[:,1:],2],axis=1)
|
||||
# print(time()-t)
|
||||
xmin,xmax = minmax(self._points[self._faces[:,1:],0])
|
||||
zmin,zmax = minmax(self._points[self._faces[:,1:],2])
|
||||
xmin,xmax = minmax(self._vertices[:,:,0])
|
||||
zmin,zmax = minmax(self._vertices[:,:,2])
|
||||
print(time()-t)
|
||||
# print(xmin.shape,xmin2.shape)
|
||||
# print(np.all(np.isclose(xmin,xmin2)))
|
||||
|
|
@ -957,9 +946,9 @@ class Features3d:
|
|||
print(is_inside.shape)
|
||||
for ii in range(Npts):
|
||||
hit_idx,t,N = Features3d.ray_triangle_intersection(pts[ii],dR,
|
||||
self._points[self._faces[bla[ii],1]],
|
||||
self._points[self._faces[bla[ii],2]],
|
||||
self._points[self._faces[bla[ii],3]])
|
||||
self._vertices[bla[ii],0,:],
|
||||
self._vertices[bla[ii],1,:],
|
||||
self._vertices[bla[ii],2,:])
|
||||
if hit_idx is None:
|
||||
is_inside[ii] = False
|
||||
else:
|
||||
|
|
@ -1019,37 +1008,36 @@ class Features3d:
|
|||
t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx]
|
||||
return hit_idx,t,N[hit_idx]
|
||||
|
||||
def clean_points(self,report=False):
|
||||
nfaces_ = self._faces.shape[0]
|
||||
ind,inv = np.unique(self._faces[:,1:],return_inverse=True)
|
||||
if report:
|
||||
print('[Features3d.clean_points]',end=' ')
|
||||
print('removed {:} orphan points.'.format(self._points.shape[0]-len(ind)))
|
||||
self._faces[:,1:4] = inv.reshape(nfaces_,3)
|
||||
self._points = self._points[ind,:]
|
||||
def faces_points_representation(self,labels,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)
|
||||
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)
|
||||
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)
|
||||
return faces,points
|
||||
|
||||
def to_vtk(self,labels,reduce_points=False):
|
||||
import pyvista as pv
|
||||
from collections.abc import Iterable
|
||||
if labels is None:
|
||||
return pv.PolyData(self._points,self._faces)
|
||||
if not isinstance(labels,Iterable):
|
||||
labels = [labels]
|
||||
points = self._points
|
||||
faces = []
|
||||
for label in labels:
|
||||
sl_ = slice(self._offset[label],self._offset[label+1])
|
||||
faces.append(self._faces[sl_,:])
|
||||
faces = np.concatenate(faces,axis=0)
|
||||
# if reduce_points:
|
||||
# nfaces_ = self._faces[idx].shape[0]
|
||||
# ind,inv = np.unique(self._faces[:,1:4],return_inverse=True)
|
||||
# faces_ = np.zeros((nfaces_,4),dtype=self._faces.dtype)
|
||||
# faces_[:,0] = 3
|
||||
# faces_[:,1:4] = inv.reshape(nfaces_,3)
|
||||
# points_ = self._points[ind,:]
|
||||
# return pv.PolyData(points_,faces_)
|
||||
# else:
|
||||
faces,points = self.faces_points_representation(labels,reduce_points=reduce_points)
|
||||
return pv.PolyData(points,faces)
|
||||
|
||||
class BinaryFieldNd:
|
||||
|
|
|
|||
Loading…
Reference in New Issue