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

Compare commits

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

2 Commits
d95453ad45 ... 2acef17323

Author SHA1 Message Date
Michael Krayer 2acef17323 faster binary hole filling 2021-08-16 21:46:40 +02:00
Michael Krayer c5908e00f5 update: inside/outside detection does not work properly 2021-08-16 12:58:23 +02:00
1 changed files with 169 additions and 200 deletions
Show Stats Expand all files Collapse all files

369
field.py
View File

@ -714,12 +714,6 @@ 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...')
@ -767,70 +761,36 @@ 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] identifying cell based labels...') if report: print('[Features3d.triangulate] sorting faces by 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]]
cell_gradient = contour.point_arrays['Gradients'][contour.faces.reshape(contour.n_faces,4)[:,1]] offset = np.zeros((self._nfeatures+1,),dtype=np.int64)
print(cell_gradient.shape) offset[1:] = np.cumsum(np.bincount(cell_labels))
# We are done with VTK for now. Since we are only dealing with triangles, let us ind = np.argsort(cell_labels)
# store them in a structure which is quicker to access and already sorted based self._offset = offset
# on labels. self._faces = contour.faces.reshape(contour.n_faces,4)[ind,:]
# _vertices has dim(ncells,nvert,ndim), where nvert=3 (three vertices per triangle) self._points = contour.points
# 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._vertices[:,0,:] A = self._points[self._faces[:,1],:]
B = self._vertices[:,1,:] B = self._points[self._faces[:,2],:]
C = self._vertices[:,2,:] C = self._points[self._faces[:,3],:]
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
@ -859,11 +819,12 @@ 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,report=False): def discard_features(self,features,clean_points=False,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:
@ -872,9 +833,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
ncells = self._vertices.shape[0] nfaces = self._faces.shape[0]
# Delete indexed elements from arrays # Delete indexed elements from arrays
self._vertices = np.delete(self._vertices,idx,axis=0) self._faces = np.delete(self._faces,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
@ -884,7 +845,20 @@ 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),ncells-self._vertices.shape[0])) len(features),nfaces-self._faces.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):
@ -915,26 +889,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._vertices[:,:,1],axis=1) min1 = np.amin(self._points[self._faces[:,1:],1],axis=1)
max1 = np.amax(self._vertices[:,:,1],axis=1) max1 = np.amax(self._points[self._faces[:,1:],1],axis=1)
min2 = np.amin(self._vertices[:,:,2],axis=1) min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
max2 = np.amax(self._vertices[:,:,2],axis=1) max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
elif kdaxis==1: elif kdaxis==1:
min1 = np.amin(self._vertices[:,:,0],axis=1) min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
max1 = np.amax(self._vertices[:,:,0],axis=1) max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
min2 = np.amin(self._vertices[:,:,2],axis=1) min2 = np.amin(self._points[self._faces[:,1:],2],axis=1)
max2 = np.amax(self._vertices[:,:,2],axis=1) max2 = np.amax(self._points[self._faces[:,1:],2],axis=1)
elif kdaxis==2: elif kdaxis==2:
min1 = np.amin(self._vertices[:,:,0],axis=1) min1 = np.amin(self._points[self._faces[:,1:],0],axis=1)
max1 = np.amax(self._vertices[:,:,0],axis=1) max1 = np.amax(self._points[self._faces[:,1:],0],axis=1)
min2 = np.amin(self._vertices[:,:,1],axis=1) min2 = np.amin(self._points[self._faces[:,1:],1],axis=1)
max2 = np.amax(self._vertices[:,:,1],axis=1) max2 = np.amax(self._points[self._faces[:,1:],1],axis=1)
else: else:
raise ValueError("Invalid ray axis.") raise ValueError("Invalid kdaxis.")
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:
@ -965,54 +939,50 @@ 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._vertices.dtype) raydir = np.zeros((3,),dtype=self._points.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._vertices[cand[ii],0,:], self._points[self._faces[cand[ii],1],:],
self._vertices[cand[ii],1,:], self._points[self._faces[cand[ii],2],:],
self._vertices[cand[ii],2,:]) self._points[self._faces[cand[ii],3],:])
# 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))
# print(cand[ii][hit_idx[idx]]) if ii==12156:
if ii==12279: for kk in hit_idx:
# 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(cand[ii][kk])
print(t[idx],hit_dir[idx],raydir) # print(self._points[self._faces[cand[ii][hit_idx[idx]],1:],:])
print(self._vertices[cand[ii][hit_idx[idx]],:,:]) # print(cand[ii][hit_idx[idx]])
hit_dir_[ii]= hit_dir[idx] # ifeat = self.feature_from_face(cand[ii][hit_idx[idx]])
face_[ii] = cand[ii][hit_idx[idx]] # print(self._offset[ifeat],self._offset[ifeat+1])
if hit_dir[idx]<0: # stop
in_feat[ii] = self.feature_from_cellidx(cand[ii][hit_idx[idx]]) # if hit_dir[idx]>0:
is_hit[ii] = 1 # in_feat[ii] = self.feature_from_face(cand[ii][hit_idx[idx]])
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 #,is_hit,hit_dir_,face_ return in_feat
def feature_from_cellidx(self,idx_cell): def feature_from_face(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
@ -1020,13 +990,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.argmax(self._offset>=idx_cell) return np.searchsorted(self._offset,idx_cell,side='right')-1
def cellidx_from_feature(self,features,concat=True): def faces_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 = None idx = slice(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:
@ -1034,7 +1004,10 @@ 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_)
idx = np.concatenate(idx,axis=0) if len(idx)==1:
idx = np.array(idx)
else:
idx = np.concatenate(idx,axis=0)
else: else:
idx = [] idx = []
for feature in features: for feature in features:
@ -1042,7 +1015,7 @@ class Features3d:
idx.append(rng_) idx.append(rng_)
return idx return idx
def ncells_from_feature(self,features): def nfaces_of_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)
@ -1136,27 +1109,8 @@ 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
@ -1181,24 +1135,31 @@ class Features3d:
icolor = (icolor+1)%ncolor icolor = (icolor+1)%ncolor
return ListedColormap(output) return ListedColormap(output)
def to_vtk(self,features,reduce_points=False,store_labels=False,remap_labels=False): def to_vtk(self,features):
'''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
faces,points = self.faces_points_representation(features,reduce_points=reduce_points) idx = self.faces_from_feature(features)
output = pv.PolyData(points,faces) return pv.PolyData(self._points,self._faces[idx,:])
if store_labels:
output.cell_arrays['label'] = np.empty(faces.shape[0],dtype=np.int) def to_vtk2(self,features):
offset = 0 import pyvista as pv
for ii,feature in enumerate(self.list_of_features(features)): idx = self.faces_from_feature(features)
ncells = self.ncells_from_feature(feature) faces_ = self._faces[idx,:]
if remap_labels: A = self._points[faces_[:,1],:]
output.cell_arrays['label'][offset:offset+ncells] = ii B = self._points[faces_[:,2],:]
else: C = self._points[faces_[:,3],:]
output.cell_arrays['label'][offset:offset+ncells] = feature cn = np.cross(B-A,C-A)
offset += ncells cn = cn/np.sqrt(np.square(cn).sum(axis=-1,keepdims=True))
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,\
@ -1207,7 +1168,6 @@ 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:
@ -1222,10 +1182,7 @@ 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)
if connect_diagonals: self.connect_diagonals = 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):
@ -1237,34 +1194,37 @@ class BinaryFieldNd:
return None return None
return self._labels[self._sldata] return self._labels[self._sldata]
def label(self): def label(self,use_cc3d=False):
'''Labels connected regions in binary fields.''' '''Labels connected regions in binary fields.'''
from scipy import ndimage if use_cc3d:
if any(self.periodicity): import cc3d
self._labels,self.nlabels,self.wrap = self._labels_periodic() 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.")
else: else:
self._labels,self.nlabels = ndimage.label(self._data,structure=self.structure) from scipy import ndimage
if self.connect_diagonals: structure = ndimage.generate_binary_structure(self._ndim,self._ndim)
def _labels_periodic(self,map_to_zero=False): else: structure = ndimage.generate_binary_structure(self._ndim,1)
'''Label features in an array while taking into account periodic wrapping. #
If map_to_zero=True, every feature which overlaps or is attached to the if use_cc3d:
periodic boundary will be removed.''' self._labels,self.nlabels = cc3d.connected_components(self._data,connectivity=connectivity,return_N=True)
from scipy import ndimage else:
# Compute labels on padded array self._labels,self.nlabels = ndimage.label(self._data,structure=structure)
labels_,nlabels_ = ndimage.label(self._data,structure=self.structure) if not any(self.periodicity):
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 = labels_[sl_lo] lab_lo = self._labels[sl_lo]
lab_hi = labels_[sl_hi] lab_hi = self._labels[sl_hi]
lab_pre = np.unique(labels_[sl_pre]) # all labels in last (unwrapped) slice lab_pre = np.unique(self._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(nlabels_+1,dtype=bool) wrap_[axis] = np.zeros(self.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]:
@ -1274,47 +1234,56 @@ 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
if map_to_zero and source_ in lab_pre: while target_ != map_[target_]: # map it recursively
map_[source_] = 0 target_ = map_[target_]
map_[target_] = 0 map_[source_] = target_
else: if source_ in lab_pre: # check if source is not a ghost
while target_ != map_[target_]: # map it recursively wrap_[axis][target_] = True
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 and remove padding # Relabel
labels_ = map_[labels_] self._labels = map_[self._labels]
nlabels_ = np.max(map_) self.nlabels = np.max(map_)
assert nlabels_==len(idx_)-1, "DEBUG assertion" self.wrap = tuple(None if x is None else x[idx_] for x in wrap_)
self.wrap = tuple(None if x is None else x[idx_] for x in self.wrap) return
# 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): def fill_holes(self,keep_wall_attached=True,use_cc3d=False,return_mask=False):
'''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. In the periodic sense, a hole is a region of zeros which is not to a boundary. When keep_wall_attached==False, only regions are kept which
connected to itself accross the periodic boundaries.''' fully connect from top to bottom wall.
from scipy import ndimage In the periodic sense, a hole is a region of zeros which is not
# Reimplementation of "binary_fill_holes" from ndimage connected to itself accross the opposite periodic boundary.'''
mask = np.logical_not(self._data) # only modify locations which are "False" at the moment if return_mask: mask = self._data.copy()
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: if self._labels is not None: self.label()
self.label() # Compute mask if it is to be returned, otherwise we are done
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.'''