First implementation of BinaryFieldNd; passed 2d tests
This commit is contained in:
parent
6f027877e2
commit
296dacd82a
206
field.py
206
field.py
|
|
@ -603,44 +603,200 @@ def gaussian_filter_umean_channel(array,spacing,sigma,truncate=4.0):
|
|||
array = ndimage.gaussian_filter1d(array,sigma_img,axis=1,truncate=truncate,mode='mirror')
|
||||
return array
|
||||
|
||||
class VoxelThreshold:
|
||||
def __init__(self,data,threshold,invert=False):
|
||||
assert isinstance(data,np.ndarray),\
|
||||
"'data' must be a numpy array."
|
||||
self._dim = data.shape
|
||||
self._ndim = data.ndim
|
||||
if invert:
|
||||
self.data = data<threshold
|
||||
else:
|
||||
self.data = data>=threshold
|
||||
|
||||
|
||||
class BinaryFieldNd:
|
||||
def __init__(self,input):
|
||||
assert isinstance(input,np.ndarray) and input.dtype==np.dtype('bool'),\
|
||||
"'input' must be a numpy array of dtype('bool')."
|
||||
self.data = input
|
||||
self._dim = input.shape
|
||||
self._ndim = input.ndim
|
||||
self.labels = None
|
||||
self.nlabels = 0
|
||||
self.wrap = tuple(self._ndim*[None])
|
||||
self.set_structure(False)
|
||||
self.set_periodicity(self._ndim*[False])
|
||||
|
||||
@classmethod
|
||||
def from_field(cls,fld3d,threshold,invert=False):
|
||||
return cls(fld3d.data,threshold,invert=invert)
|
||||
def from_threshold(cls,fld,threshold,invert=False):
|
||||
if isinstance(fld,Field3d):
|
||||
fld = fld.data
|
||||
if invert:
|
||||
return cls(fld<threshold)
|
||||
else:
|
||||
return cls(fld>=threshold)
|
||||
|
||||
def fill_holes(self,periodicity=(False,False,False)):
|
||||
'''Fills topological holes in threshold regions.'''
|
||||
def set_periodicity(self,periodicity):
|
||||
assert all([isinstance(x,(bool,int)) for x in periodicity]),\
|
||||
"'periodicity' requires bool values."
|
||||
from scipy import ndimage
|
||||
binarr = ndimage.binary_fill_holes(self.data)
|
||||
for axis in range(self._ndim):
|
||||
if periodicity[axis]:
|
||||
n = binarr.shape[axis]
|
||||
binarr = np.roll(binarr,n//2,axis=axis)
|
||||
binarr = ndimage.binary_fill_holes(binarr)
|
||||
binarr = np.roll(binarr,-n//2,axis=axis)
|
||||
self.data = binarr
|
||||
assert len(periodicity)==self._ndim,\
|
||||
"Number of entries in 'periodicity' must match dimension of binary field."
|
||||
self.periodicity = tuple(bool(x) for x in periodicity)
|
||||
return
|
||||
|
||||
def set_structure(self,connect_diagonals):
|
||||
from scipy import ndimage
|
||||
if connect_diagonals:
|
||||
self.structure = ndimage.generate_binary_structure(self._ndim,self._ndim)
|
||||
else:
|
||||
self.structure = ndimage.generate_binary_structure(self._ndim,1)
|
||||
|
||||
def enable_diagonal_connections(self): self.set_structure(True)
|
||||
|
||||
def disable_diagonal_connections(self): self.set_structure(False)
|
||||
|
||||
def label(self):
|
||||
'''Labels connected regions in binary fields.'''
|
||||
from scipy import ndimage
|
||||
if any(self.periodicity):
|
||||
self.labels,self.nlabels,self.wrap = self._labels_periodic()
|
||||
else:
|
||||
self.labels,self.nlabels = ndimage.label(self.data,structure=self.structure)
|
||||
|
||||
def _labels_periodic(self,map_to_zero=False):
|
||||
'''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
|
||||
periodic boundary will be removed.'''
|
||||
from scipy import ndimage
|
||||
# Pad input data
|
||||
if map_to_zero:
|
||||
pw = tuple((1,1) if x else (0,0) for x in self.periodicity)
|
||||
sl_pad = tuple(slice(1,-1) if x else slice(None) for x in self.periodicity)
|
||||
else:
|
||||
pw = tuple((0,1) if x else (0,0) for x in self.periodicity)
|
||||
sl_pad = tuple(slice(0,-1) if x else slice(None) for x in self.periodicity)
|
||||
data_ = np.pad(self.data,pw,mode='wrap')
|
||||
# Compute labels on padded array
|
||||
labels_,nlabels_ = ndimage.label(data_,structure=self.structure)
|
||||
# Get a mapping of labels which differ at periodic overlap
|
||||
map_ = np.array(range(nlabels_+1),dtype=labels_.dtype)
|
||||
wrap_ = self._ndim*[None]
|
||||
for axis in range(self._ndim):
|
||||
if not self.periodicity[axis]: continue
|
||||
if map_to_zero:
|
||||
sl_lo = tuple(slice(0,2) if ii==axis else slice(None) for ii in range(self._ndim))
|
||||
sl_hi = tuple(slice(-2,None) if ii==axis else slice(None) for ii in range(self._ndim))
|
||||
lab_lo = labels_[sl_lo]
|
||||
lab_hi = labels_[sl_hi]
|
||||
li = (lab_lo!=lab_hi)
|
||||
for source_ in np.unique(lab_lo[li]):
|
||||
map_[source_] = 0
|
||||
for source_ in np.unique(lab_hi[li]):
|
||||
map_[source_] = 0
|
||||
else:
|
||||
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))
|
||||
lab_lo = labels_[sl_lo]
|
||||
lab_hi = labels_[sl_hi]
|
||||
# Initialize array to keep track of wrapping
|
||||
wrap_[axis] = np.zeros(nlabels_+1,dtype=bool)
|
||||
# Determine new label and map
|
||||
lab_new = np.minimum(lab_lo,lab_hi)
|
||||
for lab_ in [lab_lo,lab_hi]:
|
||||
li = (lab_!=lab_new)
|
||||
lab_li = lab_[li]
|
||||
lab_new_li = lab_new[li]
|
||||
for idx_ in np.unique(lab_li,return_index=True)[1]:
|
||||
source_ = lab_li[idx_] # the label to be changed
|
||||
target_ = lab_new_li[idx_] # the label which will be newly assigned
|
||||
while target_ != map_[target_]: # map it recursively
|
||||
target_ = map_[target_]
|
||||
map_[source_] = target_
|
||||
wrap_[axis][target_] = True
|
||||
# Remove gaps from target mapping
|
||||
idx_,map_ = np.unique(map_,return_index=True,return_inverse=True)[1:3]
|
||||
# Relabel and remove padding
|
||||
labels_ = map_[labels_[sl_pad]]
|
||||
nlabels_ = np.max(map_)
|
||||
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):
|
||||
'''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
|
||||
to a boundary. In the periodic sense, a hole is a region of zeros which is not
|
||||
connected to itself accross the periodic boundaries.'''
|
||||
from scipy import ndimage
|
||||
# Reimplementation of "binary_fill_holes" from ndimage
|
||||
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)
|
||||
|
||||
def probe(self,idx):
|
||||
'''Returns whether or not point at index is inside threshold region or not.'''
|
||||
'''Returns whether or not a point at idx is True or False.'''
|
||||
return self.data[tuple(idx)]
|
||||
|
||||
def volume(self):
|
||||
'''Returns volume of region above threshold.'''
|
||||
'''Returns the sum of True values.'''
|
||||
return np.sum(self.data)
|
||||
|
||||
def volume_feature(self,label=None):
|
||||
'''Returns volume of features, i.e. connected regions which have been
|
||||
labeled using the label() method. If 'label' is None all volumes
|
||||
are returned including the volume of the background region. The array
|
||||
is sorted by labels, i.e. vol[0] is the volume of the background region,
|
||||
vol[1] the volume of label 1, etc. If 'label' is an integer value, only
|
||||
the volume of the corresponding region is returned.
|
||||
Note: it is more efficient to retrieve all volumes at once than querying
|
||||
single labels.'''
|
||||
if self.labels is None:
|
||||
self.label()
|
||||
if label is None:
|
||||
return np.bincount(self.labels.ravel())
|
||||
else:
|
||||
return np.sum(self.labels==label)
|
||||
|
||||
def volume_domain(self):
|
||||
'''Returns volume of entire domain. Should be equal to sum(volume_feature()).'''
|
||||
return np.prod(self._dim)
|
||||
|
||||
def feature_labels_by_volume(self,descending=True):
|
||||
'''Returns labels of connected regions sorted by volume.'''
|
||||
labels = np.argsort(self.volume_feature()[1:])+1
|
||||
if descending: labels = labels[::-1]
|
||||
return labels
|
||||
|
||||
def discard_feature(self,selection):
|
||||
if self.labels is None:
|
||||
self.label()
|
||||
selection = self._select_feature(selection)
|
||||
# Map tagged regions to zero in order to discard them
|
||||
map_ = np.array(range(0,self.count+1),dtype=self.label.dtype)
|
||||
map_[selection] = 0
|
||||
# Remove gaps from target mapping
|
||||
map_ = np.unique(map_,return_inverse=True)[1]
|
||||
# Discard regions
|
||||
self.labels = map_[self.labels]
|
||||
self.nlabels = np.max(map_)
|
||||
self.data = self.labels>0
|
||||
|
||||
def _select_feature(self,selection):
|
||||
dtype = self.label.dtype
|
||||
if np.issubdtype(type(selection),np.integer):
|
||||
return np.array(selection,dtype=dtype)
|
||||
elif isinstance(selection,(list,tuple,np.ndarray)):
|
||||
selection = np.array(selection)
|
||||
if selection.dtype==np.dtype('bool'):
|
||||
assert selection.ndim==1 and selection.shape[0]==self.count+1,\
|
||||
"Boolean indexing must provide count+1 values."
|
||||
else:
|
||||
selection = selection.astype(dtype)
|
||||
assert np.max(selection)<=self.count and np.min(selection)>=0,\
|
||||
"Entry in selection is out-of-bounds."
|
||||
return selection
|
||||
else:
|
||||
raise ValueError('Invalid input. Accepting int,list,tuple,ndarray.')
|
||||
|
||||
class ConnectedRegions:
|
||||
def __init__(self,binarr,periodicity,connect_diagonals=False,fill_holes=False,bytes_label=32):
|
||||
assert isinstance(binarr,np.ndarray) and binarr.dtype==np.dtype('bool'),\
|
||||
|
|
|
|||
Loading…
Reference in New Issue