mwtkrayer
/
suspendtools
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: mwtkrayer:2acef17323644bd068ff22782f690384e05e1ce5
Branches Tags
mwtkrayer:master
..
compare: mwtkrayer:d95453ad45408c85cdc14dba27f29187efc493bb
Branches Tags
mwtkrayer:master

No commits in common. "2acef17323644bd068ff22782f690384e05e1ce5" and "d95453ad45408c85cdc14dba27f29187efc493bb" have entirely different histories.
2acef17323 ... d95453ad45

1 changed files with 200 additions and 169 deletions
Show Stats Expand all files Collapse all files

369
field.py
View File

@ -714,6 +714,12 @@ class Features3d:
if report: print('[Features3d.triangulate] computing isocontour using {}...'.format(contour_method)) 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) contour = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False,compute_gradients=True)
assert contour.is_all_triangles(), "Contouring produced non-triangle cells." 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 # Compute the connectivity of the triangulated surface: first we run an ordinary
# connectivity filter neglecting periodic wrapping. # connectivity filter neglecting periodic wrapping.
if report: print('[Features3d.triangulate] computing connectivity...') if report: print('[Features3d.triangulate] computing connectivity...')
@ -761,36 +767,70 @@ class Features3d:
while target_ != map_[target_]: # map it recursively while target_ != map_[target_]: # map it recursively
target_ = map_[target_] target_ = map_[target_]
map_[source_] = 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] point_labels = map_[point_labels]
self._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 # 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. # point of each cell determines the value for this cell.
if report: print('[Features3d.triangulate] sorting faces by labels...') if report: print('[Features3d.triangulate] identifying cell based labels...')
cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]] cell_labels = point_labels[contour.faces.reshape(contour.n_faces,4)[:,1]]
offset = np.zeros((self._nfeatures+1,),dtype=np.int64) cell_gradient = contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1]]
offset[1:] = np.cumsum(np.bincount(cell_labels)) print(cell_gradient.shape)
ind = np.argsort(cell_labels) # We are done with VTK for now. Since we are only dealing with triangles, let us
self._offset = offset # store them in a structure which is quicker to access and already sorted based
self._faces = contour.faces.reshape(contour.n_faces,4)[ind,:] # on labels.
self._points = contour.points # _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:]]
# Compute the volume and area per cell. For the volume computation, an arbitrary component # 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 # of the normal has to be chosen which defaults to the z-component and is set by
# 'cellvol_normal_component'. # 'cellvol_normal_component'.
if report: print('[Features3d.triangulate] calculating area and volume per cell...') if report: print('[Features3d.triangulate] calculating area and volume per cell...')
A = self._points[self._faces[:,1],:] A = self._vertices[:,0,:]
B = self._points[self._faces[:,2],:] B = self._vertices[:,1,:]
C = self._points[self._faces[:,3],:] C = self._vertices[:,2,:]
cn = np.cross(B-A,C-A) 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 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_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._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 return
@property @property
@ -819,12 +859,11 @@ class Features3d:
self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report) self.discard_features(np.flatnonzero(np.abs(self.volumes())<threshold),report=report)
return return
def discard_features(self,features,clean_points=False,report=False): def discard_features(self,features,report=False):
'''Removes features from triangulated data.''' '''Removes features from triangulated data.'''
features = self.list_of_features(features) features = self.list_of_features(features)
# Get index ranges which are to be deleted and also create an array # Get index ranges which are to be deleted and also create an array
# which determines the size of the block to be deleted # which determines the size of the block to be deleted
print('Discarding',features)
idx = [] idx = []
gapsize = np.zeros((self._nfeatures,),dtype=np.int) gapsize = np.zeros((self._nfeatures,),dtype=np.int)
for feature in features: for feature in features:
@ -833,9 +872,9 @@ class Features3d:
gapsize[feature] = len(rng_) gapsize[feature] = len(rng_)
idx = np.concatenate(idx,axis=0) idx = np.concatenate(idx,axis=0)
# Save former number of faces for report # Save former number of faces for report
nfaces = self._faces.shape[0] ncells = self._vertices.shape[0]
# Delete indexed elements from arrays # 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_areas = np.delete(self._cell_areas,idx,axis=0)
self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0) self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0)
# Correct offset # Correct offset
@ -845,20 +884,7 @@ class Features3d:
if report: if report:
print('[Features3d.discard_features]',end=' ') print('[Features3d.discard_features]',end=' ')
print('discarded {:} features with {:} faces.'.format( print('discarded {:} features with {:} faces.'.format(
len(features),nfaces-self._faces.shape[0])) len(features),ncells-self._vertices.shape[0]))
if clean_points:
self.clean_points(report=report)
return
def clean_points(self,report=False):
'''Removes points which are not referenced by any face.'''
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,:]
return return
def area(self,feature): def area(self,feature):
@ -889,26 +915,26 @@ class Features3d:
'''Builds a KD-tree for feature search.''' '''Builds a KD-tree for feature search.'''
from scipy import spatial from scipy import spatial
if kdaxis==0: if kdaxis==0:
min1 = np.amin(self._points[self._faces[:,1:],1],axis=1) min1 = np.amin(self._vertices[:,:,1],axis=1)
max1 = np.amax(self._points[self._faces[:,1:],1],axis=1) max1 = np.amax(self._vertices[:,:,1],axis=1)
min2 = np.amin(self._points[self._faces[:,1:],2],axis=1) min2 = np.amin(self._vertices[:,:,2],axis=1)
max2 = np.amax(self._points[self._faces[:,1:],2],axis=1) max2 = np.amax(self._vertices[:,:,2],axis=1)
elif kdaxis==1: elif kdaxis==1:
min1 = np.amin(self._points[self._faces[:,1:],0],axis=1) min1 = np.amin(self._vertices[:,:,0],axis=1)
max1 = np.amax(self._points[self._faces[:,1:],0],axis=1) max1 = np.amax(self._vertices[:,:,0],axis=1)
min2 = np.amin(self._points[self._faces[:,1:],2],axis=1) min2 = np.amin(self._vertices[:,:,2],axis=1)
max2 = np.amax(self._points[self._faces[:,1:],2],axis=1) max2 = np.amax(self._vertices[:,:,2],axis=1)
elif kdaxis==2: elif kdaxis==2:
min1 = np.amin(self._points[self._faces[:,1:],0],axis=1) min1 = np.amin(self._vertices[:,:,0],axis=1)
max1 = np.amax(self._points[self._faces[:,1:],0],axis=1) max1 = np.amax(self._vertices[:,:,0],axis=1)
min2 = np.amin(self._points[self._faces[:,1:],1],axis=1) min2 = np.amin(self._vertices[:,:,1],axis=1)
max2 = np.amax(self._points[self._faces[:,1:],1],axis=1) max2 = np.amax(self._vertices[:,:,1],axis=1)
else: else:
raise ValueError("Invalid kdaxis.") raise ValueError("Invalid ray axis.")
center = np.stack((0.5*(max1+min1),0.5*(max2+min2)),axis=1) center = np.stack((0.5*(max1+min1),0.5*(max2+min2)),axis=1)
radius = 0.5*np.amax(np.sqrt((max1-min1)**2+(max2-min2)**2)) radius = 0.5*np.amax(np.sqrt((max1-min1)**2+(max2-min2)**2))
self._kdtree = spatial.KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes, self._kdtree = spatial.KDTree(center,leafsize=leafsize,compact_nodes=compact_nodes,
copy_data=False,balanced_tree=balanced_tree) copy_data=False,balanced_tree=balanced_tree)
self._kdaxis = kdaxis self._kdaxis = kdaxis
self._kdradius = radius self._kdradius = radius
if report: if report:
@ -939,50 +965,54 @@ class Features3d:
if not self._kdtree: if not self._kdtree:
self.build_kdtree() self.build_kdtree()
if self._kdaxis==0: query_axis = [1,2] if self._kdaxis==0: query_axis = [1,2]
elif self._kdaxis==1: query_axis = [0,2] elif self._kdaxis==1: query_axis = [0,2]
elif self._kdaxis==2: query_axis = [0,1] elif self._kdaxis==2: query_axis = [0,1]
cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius) cand = self._kdtree.query_ball_point(coords[:,query_axis],self._kdradius)
if coords.shape[0]>12279:
print('Point coordinates:',coords[12279,:])
Ncoord = coords.shape[0] Ncoord = coords.shape[0]
raydir = np.zeros((3,),dtype=self._points.dtype) raydir = np.zeros((3,),dtype=self._vertices.dtype)
raydir[self._kdaxis] = 1.0 raydir[self._kdaxis] = 1.0
in_feat = np.empty((Ncoord,),dtype=np.int) in_feat = np.empty((Ncoord,),dtype=np.int)
is_hit = np.empty((Ncoord,),dtype=np.int) is_hit = np.empty((Ncoord,),dtype=np.int)
hit_dir_ = np.empty((Ncoord,),dtype=np.int) hit_dir_ = np.empty((Ncoord,),dtype=np.int)
face_ = np.empty((Ncoord,),dtype=np.int) face_ = np.empty((Ncoord,),dtype=np.int)
print('point 12156 = ',coords[12156])
for ii in range(Ncoord): for ii in range(Ncoord):
hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir, hit_idx,t,hit_dir = Features3d.ray_triangle_intersection(coords[ii],raydir,
self._points[self._faces[cand[ii],1],:], self._vertices[cand[ii],0,:],
self._points[self._faces[cand[ii],2],:], self._vertices[cand[ii],1,:],
self._points[self._faces[cand[ii],3],:]) self._vertices[cand[ii],2,:])
# if ii==12279: print('DEBUG',hit_idx,t,N) # if ii==12279: print('DEBUG',hit_idx,t,N)
if ii==12156: print('DEBUG',hit_idx,t,hit_dir)
if hit_idx is None: if hit_idx is None:
in_feat[ii] = -1 in_feat[ii] = -1
is_hit[ii] = -1
hit_dir_[ii]= 0
face_[ii] = -1
else: else:
idx = np.argmin(np.abs(t)) idx = np.argmin(np.abs(t))
if ii==12156: # print(cand[ii][hit_idx[idx]])
for kk in hit_idx: if ii==12279:
print(cand[ii][kk]) # 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(self._points[self._faces[cand[ii][hit_idx[idx]],1:],:]) print(t[idx],hit_dir[idx],raydir)
# print(cand[ii][hit_idx[idx]]) print(self._vertices[cand[ii][hit_idx[idx]],:,:])
# ifeat = self.feature_from_face(cand[ii][hit_idx[idx]]) hit_dir_[ii]= hit_dir[idx]
# print(self._offset[ifeat],self._offset[ifeat+1]) face_[ii] = cand[ii][hit_idx[idx]]
# stop if hit_dir[idx]<0:
# if hit_dir[idx]>0: in_feat[ii] = self.feature_from_cellidx(cand[ii][hit_idx[idx]])
# in_feat[ii] = self.feature_from_face(cand[ii][hit_idx[idx]]) is_hit[ii] = 1
# else: else:
# in_feat[ii] = -1 in_feat[ii] = -1
is_hit[ii] = 0
if report: if report:
print('[Features3d.inside_feature]',end=' ') print('[Features3d.inside_feature]',end=' ')
print('{} of {} points are located inside of features.'.format(sum(in_feat>=0),Ncoord)) print('{} of {} points are located inside of features.'.format(sum(in_feat>=0),Ncoord))
return in_feat return in_feat #,is_hit,hit_dir_,face_
def feature_from_face(self,idx_cell): def feature_from_cellidx(self,idx_cell):
'''Gets feature ID for a given cell.''' '''Gets feature ID for a given cell.'''
# This is an optimized solution for monotonically increasing arrays: # This is an optimized solution for monotonically increasing arrays:
# as long as offset is smaller than the cell index, the boolean operation # as long as offset is smaller than the cell index, the boolean operation
@ -990,13 +1020,13 @@ class Features3d:
# as soon as idx_cell is equal or larger than offset, the argmax function # 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 # short circuits and returns the index of the first occurence without
# checking any value afterwards. # checking any value afterwards.
return np.searchsorted(self._offset,idx_cell,side='right')-1 return np.argmax(self._offset>=idx_cell)
def faces_from_feature(self,features,concat=True): def cellidx_from_feature(self,features,concat=True):
'''Returns indices of cells which belong to given features.''' '''Returns indices of cells which belong to given features.'''
from collections.abc import Iterable from collections.abc import Iterable
if features is None: if features is None:
idx = slice(None) idx = None
elif not isinstance(features,Iterable): elif not isinstance(features,Iterable):
idx = slice(self._offset[features],self._offset[features+1]) idx = slice(self._offset[features],self._offset[features+1])
elif concat: elif concat:
@ -1004,10 +1034,7 @@ class Features3d:
for feature in features: for feature in features:
rng_ = np.arange(self._offset[feature],self._offset[feature+1]) rng_ = np.arange(self._offset[feature],self._offset[feature+1])
idx.append(rng_) idx.append(rng_)
if len(idx)==1: idx = np.concatenate(idx,axis=0)
idx = np.array(idx)
else:
idx = np.concatenate(idx,axis=0)
else: else:
idx = [] idx = []
for feature in features: for feature in features:
@ -1015,7 +1042,7 @@ class Features3d:
idx.append(rng_) idx.append(rng_)
return idx return idx
def nfaces_of_feature(self,features): def ncells_from_feature(self,features):
'''Returns number of cells which belong to given features. Can be used '''Returns number of cells which belong to given features. Can be used
to construct new offsets.''' to construct new offsets.'''
features = self.list_of_features(features) features = self.list_of_features(features)
@ -1109,8 +1136,27 @@ class Features3d:
if len(hit_idx)==0: return (None,None,None) if len(hit_idx)==0: return (None,None,None)
# Intersection point is R+t*dR # Intersection point is R+t*dR
t = (AO[hit_idx,:]*N[hit_idx,:]).sum(axis=-1)*invdet[hit_idx] 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]) return hit_idx,t,np.sign(t)*np.sign(det[hit_idx])
return hit_idx,t,N[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
def cmap_features(self,nfeatures,name='tab10'): def cmap_features(self,nfeatures,name='tab10'):
from matplotlib.colors import ListedColormap from matplotlib.colors import ListedColormap
@ -1135,31 +1181,24 @@ class Features3d:
icolor = (icolor+1)%ncolor icolor = (icolor+1)%ncolor
return ListedColormap(output) return ListedColormap(output)
def to_vtk(self,features): 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 import pyvista as pv
idx = self.faces_from_feature(features) faces,points = self.faces_points_representation(features,reduce_points=reduce_points)
return pv.PolyData(self._points,self._faces[idx,:]) output = pv.PolyData(points,faces)
if store_labels:
def to_vtk2(self,features): output.cell_arrays['label'] = np.empty(faces.shape[0],dtype=np.int)
import pyvista as pv offset = 0
idx = self.faces_from_feature(features) for ii,feature in enumerate(self.list_of_features(features)):
faces_ = self._faces[idx,:] ncells = self.ncells_from_feature(feature)
A = self._points[faces_[:,1],:] if remap_labels:
B = self._points[faces_[:,2],:] output.cell_arrays['label'][offset:offset+ncells] = ii
C = self._points[faces_[:,3],:] else:
cn = np.cross(B-A,C-A) output.cell_arrays['label'][offset:offset+ncells] = feature
cn = cn/np.sqrt(np.square(cn).sum(axis=-1,keepdims=True)) offset += ncells
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 return output
def vtk_faces(self,face_idx):
import pyvista as pv
return pv.PolyData(self._points,self._faces[face_idx,:])
class BinaryFieldNd: class BinaryFieldNd:
def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False): def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False):
assert isinstance(input,np.ndarray) and input.dtype==bool,\ assert isinstance(input,np.ndarray) and input.dtype==bool,\
@ -1168,6 +1207,7 @@ class BinaryFieldNd:
"'periodicity' requires bool values." "'periodicity' requires bool values."
assert len(periodicity)==input.ndim,\ assert len(periodicity)==input.ndim,\
"Number of entries in 'periodicity' must match dimension of binary field." "Number of entries in 'periodicity' must match dimension of binary field."
from scipy import ndimage
if has_ghost and deep: if has_ghost and deep:
self._data = input.copy() self._data = input.copy()
elif has_ghost: elif has_ghost:
@ -1182,7 +1222,10 @@ class BinaryFieldNd:
self.nlabels = None self.nlabels = None
self.wrap = tuple(self._ndim*[None]) self.wrap = tuple(self._ndim*[None])
self.periodicity = tuple(bool(x) for x in periodicity) self.periodicity = tuple(bool(x) for x in periodicity)
self.connect_diagonals = connect_diagonals if connect_diagonals:
self.structure = ndimage.generate_binary_structure(self._ndim,self._ndim)
else:
self.structure = ndimage.generate_binary_structure(self._ndim,1)
@property @property
def data(self): def data(self):
@ -1194,37 +1237,34 @@ class BinaryFieldNd:
return None return None
return self._labels[self._sldata] return self._labels[self._sldata]
def label(self,use_cc3d=False): def label(self):
'''Labels connected regions in binary fields.''' '''Labels connected regions in binary fields.'''
if use_cc3d: from scipy import ndimage
import cc3d if any(self.periodicity):
if self._ndim==2: connectivity = 8 if self.connect_diagonals else 4 self._labels,self.nlabels,self.wrap = self._labels_periodic()
elif self._ndim==3: connectivity = 18 if self.connect_diagonals else 6
else: raise RuntimeError("'use_cc3d' can only be used with 2D or 3D data.")
else: else:
from scipy import ndimage self._labels,self.nlabels = ndimage.label(self._data,structure=self.structure)
if self.connect_diagonals: structure = ndimage.generate_binary_structure(self._ndim,self._ndim)
else: structure = ndimage.generate_binary_structure(self._ndim,1) def _labels_periodic(self,map_to_zero=False):
# '''Label features in an array while taking into account periodic wrapping.
if use_cc3d: If map_to_zero=True, every feature which overlaps or is attached to the
self._labels,self.nlabels = cc3d.connected_components(self._data,connectivity=connectivity,return_N=True) periodic boundary will be removed.'''
else: from scipy import ndimage
self._labels,self.nlabels = ndimage.label(self._data,structure=structure) # Compute labels on padded array
if not any(self.periodicity): labels_,nlabels_ = ndimage.label(self._data,structure=self.structure)
return
# Get a mapping of labels which differ at periodic overlap # Get a mapping of labels which differ at periodic overlap
map_ = np.array(range(nlabels_+1),dtype=labels_.dtype) map_ = np.array(range(nlabels_+1),dtype=labels_.dtype)
wrap_ = self._ndim*[None] wrap_ = self._ndim*[None]
for axis in range(self._ndim): for axis in range(self._ndim):
if not self.periodicity[axis]: continue if not self.periodicity[axis]: continue
sl_lo = tuple(slice(0,1) if ii==axis else slice(None) for ii in range(self._ndim)) sl_lo = tuple(slice(0,1) if ii==axis else slice(None) for ii in range(self._ndim))
sl_hi = tuple(slice(-1,None) if ii==axis else slice(None) for ii in range(self._ndim)) sl_hi = tuple(slice(-1,None) if ii==axis else slice(None) for ii in range(self._ndim))
sl_pre = tuple(slice(-2,-1) if ii==axis else slice(None) for ii in range(self._ndim)) sl_pre = tuple(slice(-2,-1) if ii==axis else slice(None) for ii in range(self._ndim))
lab_lo = self._labels[sl_lo] lab_lo = labels_[sl_lo]
lab_hi = self._labels[sl_hi] lab_hi = labels_[sl_hi]
lab_pre = np.unique(self._labels[sl_pre]) # all labels in last (unwrapped) slice lab_pre = np.unique(labels_[sl_pre]) # all labels in last (unwrapped) slice
# Initialize array to keep track of wrapping # Initialize array to keep track of wrapping
wrap_[axis] = np.zeros(self.nlabels+1,dtype=bool) wrap_[axis] = np.zeros(nlabels_+1,dtype=bool)
# Determine new label and map # Determine new label and map
lab_new = np.minimum(lab_lo,lab_hi) lab_new = np.minimum(lab_lo,lab_hi)
for lab_ in [lab_lo,lab_hi]: for lab_ in [lab_lo,lab_hi]:
@ -1234,56 +1274,47 @@ class BinaryFieldNd:
for idx_ in np.unique(lab_li,return_index=True)[1]: for idx_ in np.unique(lab_li,return_index=True)[1]:
source_ = lab_li[idx_] # the label to be changed source_ = lab_li[idx_] # the label to be changed
target_ = lab_new_li[idx_] # the label which will be newly assigned target_ = lab_new_li[idx_] # the label which will be newly assigned
while target_ != map_[target_]: # map it recursively if map_to_zero and source_ in lab_pre:
target_ = map_[target_] map_[source_] = 0
map_[source_] = target_ map_[target_] = 0
if source_ in lab_pre: # check if source is not a ghost else:
wrap_[axis][target_] = True while target_ != map_[target_]: # map it recursively
target_ = map_[target_]
map_[source_] = target_
if source_ in lab_pre: # check if source is not a ghost
wrap_[axis][target_] = True
# Remove gaps from target mapping # Remove gaps from target mapping
idx_,map_ = np.unique(map_,return_index=True,return_inverse=True)[1:3] idx_,map_ = np.unique(map_,return_index=True,return_inverse=True)[1:3]
# Relabel # Relabel and remove padding
self._labels = map_[self._labels] labels_ = map_[labels_]
self.nlabels = np.max(map_) nlabels_ = np.max(map_)
self.wrap = tuple(None if x is None else x[idx_] for x in wrap_) assert nlabels_==len(idx_)-1, "DEBUG assertion"
return self.wrap = tuple(None if x is None else x[idx_] for x in self.wrap)
# for axis in range(self._ndim):
# if wrap_[axis] is not None:
# wrap_[axis] = wrap_[axis][idx_]
return labels_,nlabels_,tuple(wrap_)
def fill_holes(self,keep_wall_attached=True,use_cc3d=False,return_mask=False): def fill_holes(self):
'''Fill the holes in binary objects while taking into account periodicity. '''Fill the holes in binary objects while taking into account periodicity.
In the non-periodic sense, a hole is a region of zeros which does not connect In the non-periodic sense, a hole is a region of zeros which does not connect
to a boundary. When keep_wall_attached==False, only regions are kept which to a boundary. In the periodic sense, a hole is a region of zeros which is not
fully connect from top to bottom wall. connected to itself accross the periodic boundaries.'''
In the periodic sense, a hole is a region of zeros which is not from scipy import ndimage
connected to itself accross the opposite periodic boundary.''' # Reimplementation of "binary_fill_holes" from ndimage
if return_mask: mask = self._data.copy() mask = np.logical_not(self._data) # only modify locations which are "False" at the moment
tmp = np.zeros(mask.shape,bool) # create empty array to "grow from boundaries"
ndimage.binary_dilation(tmp,structure=None,iterations=-1,
mask=mask,output=self._data,border_value=1,
origin=0) # everything connected to the boundary is now True in self._data
# Remove holes which overlap the boundaries
if any(self.periodicity):
self._data = self._labels_periodic(map_to_zero=True)[0]>0
# Invert to get the final result
np.logical_not(self._data,self._data) np.logical_not(self._data,self._data)
if use_cc3d:
import cc3d
if self._ndim==2: connectivity = 8 if self.connect_diagonals else 4
elif self._ndim==3: connectivity = 18 if self.connect_diagonals else 6
else: raise RuntimeError("'use_cc3d' can only be used with 2D or 3D data.")
labels_,nlabels_ = cc3d.connected_components(self._data,connectivity=connectivity,return_N=True)
else:
from scipy import ndimage
if self.connect_diagonals: structure = ndimage.generate_binary_structure(self._ndim,self._ndim)
else: structure = ndimage.generate_binary_structure(self._ndim,1)
labels_,nlabels_ = ndimage.label(self._data,structure=structure)
labels_keep = set()
for axis in range(3):
sl_lo = tuple(slice(None) if ii!=axis else 0 for ii in range(3))
sl_hi = tuple(slice(None) if ii!=axis else -1 for ii in range(3))
if self.periodicity[axis]:
labels_keep |= set(np.unique(labels_[sl_lo])) & set(np.unique(labels_[sl_hi]))
elif keep_wall_attached:
labels_keep |= set(np.unique(labels_[sl_lo])) | set(np.unique(labels_[sl_hi]))
labels_keep.discard(0)
self._data = np.isin(labels_,tuple(labels_keep),invert=True)
# If labels have been computed already, recompute them to stay consistent # If labels have been computed already, recompute them to stay consistent
if self._labels is not None: self.label() if self._labels is not None:
# Compute mask if it is to be returned, otherwise we are done self.label()
if return_mask:
np.logical_xor(mask,self._data,out=mask)
return mask
return
def probe(self,idx,probe_label=False): def probe(self,idx,probe_label=False):
'''Returns whether or not a point at idx is True or False.''' '''Returns whether or not a point at idx is True or False.'''