implemented surface area calculation in binary data. Does not work too well.
This commit is contained in:
parent
d7658ea659
commit
56b413cfed
54
field.py
54
field.py
|
|
@ -790,20 +790,46 @@ class BinaryFieldNd:
|
|||
the volume of the corresponding region is returned.
|
||||
Note: it is more efficient to retrieve all volumes at once than querying
|
||||
single labels.'''
|
||||
from scipy import ndimage
|
||||
if self.labels is None:
|
||||
self.label()
|
||||
if label is None:
|
||||
return np.bincount(self.labels.ravel())
|
||||
# TBD: compare performance with
|
||||
# sizes = ndimage.sum(mask, label_im, range(nb_labels + 1))
|
||||
# http://scipy-lectures.org/advanced/image_processing/auto_examples/plot_find_object.html
|
||||
# return ndimage.sum_labels(self.data,labels=self.labels,index=range(0,self.nlabels+1)) # performs much worse
|
||||
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 area_feature(self,ignore_curved=False,ignore_pathological=False):
|
||||
from scipy import ndimage
|
||||
from time import time
|
||||
# Pad array in periodic directions
|
||||
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)
|
||||
data_ = np.pad(self.data,pw,mode='wrap')
|
||||
#
|
||||
area = np.zeros(self.nlabels+1)
|
||||
stypes = list(range(1,10))
|
||||
if ignore_pathological: del stypes[6:9]
|
||||
if ignore_curved: del stypes[3:6]
|
||||
print(stypes,ignore_curved,ignore_pathological)
|
||||
for stype in stypes:
|
||||
print(stype)
|
||||
svs_,w_ = BinaryFieldNd.surface_voxel_structure(stype)
|
||||
hit_ = np.zeros_like(data_)
|
||||
t = time()
|
||||
ii=0
|
||||
for s_ in svs_: # Loop over permutations
|
||||
ii+=1
|
||||
np.logical_or(hit_,ndimage.binary_hit_or_miss(data_,structure1=s_),out=hit_) # Check if this works inplace
|
||||
print(ii,time()-t,np.sum(hit_))
|
||||
area += w_*np.bincount(self.labels[hit_[sl_pad]],minlength=self.nlabels+1)
|
||||
return area
|
||||
|
||||
def feature_labels_by_volume(self,descending=True):
|
||||
'''Returns labels of connected regions sorted by volume.'''
|
||||
labels = np.argsort(self.volume_feature()[1:])+1
|
||||
|
|
@ -923,40 +949,42 @@ class BinaryFieldNd:
|
|||
raise ValueError('Invalid input. Accepting int,list,tuple,ndarray.')
|
||||
|
||||
@staticmethod
|
||||
def surface_voxel_structure(id):
|
||||
def surface_voxel_structure(stype):
|
||||
'''Get structure element and all its permutations for surface voxels
|
||||
as defined by
|
||||
Mullikin, J. C., & Verbeek, P. W. (1993). Surface area estimation of digitized planes. Bioimaging, 1(1), 6-16.
|
||||
Also returns weights, where S7-S9 are defined in
|
||||
Windreich, G., Kiryati, N., & Lohmann, G. (2003). Voxel-based surface area estimation: from theory to practice. Pattern Recognition, 36(11), 2531-2541.
|
||||
'''
|
||||
if id==1:
|
||||
if stype==1:
|
||||
S = np.array([[[ True, True, True], [ True, True, True], [ True, True, True]],[[ True, True, True], [ True, True, True], [ True, True, True]],[[False,False,False], [False,False,False], [False,False,False]]])
|
||||
W = 0.894
|
||||
elif id==2:
|
||||
elif stype==2:
|
||||
S = np.array([[[ True, True, True], [ True, True, True], [ True, True, True]],[[ True, True,False], [ True, True,False], [ True, True,False]],[[ True,False,False], [ True,False,False], [ True,False,False]]])
|
||||
W = 1.3409
|
||||
elif id==3:
|
||||
elif stype==3:
|
||||
S = np.array([[[ True, True, True], [ True, True,False], [ True,False,False]],[[ True, True,False], [ True, True,False], [False,False,False]],[[ True,False,False], [False,False,False], [False,False,False]]])
|
||||
W = 1.5879
|
||||
elif id==4:
|
||||
elif stype==4:
|
||||
S = np.array([[[ True, True, True], [ True, True, True], [ True, True, True]],[[False, True,False], [False, True,False], [False, True,False]],[[False,False,False], [False,False,False], [False,False,False]]])
|
||||
W = 2.0
|
||||
elif id==5:
|
||||
elif stype==5:
|
||||
S = np.array([[[ True,False,False], [False,False,False], [False,False,False]],[[ True, True,False], [ True, True,False], [False,False,False]],[[ True,False,False], [False,False,False], [False,False,False]]])
|
||||
W = 8./3.
|
||||
elif id==6:
|
||||
elif stype==6:
|
||||
S = np.array([[[False,False,False], [ True,False,False], [False,False,False]],[[False,False,False], [ True, True,False], [False,False,False]],[[False,False,False], [ True,False,False], [False,False,False]]])
|
||||
W = 10./3.
|
||||
elif id==7:
|
||||
elif stype==7:
|
||||
S = np.array([[[False,False,False], [False,False,False], [False,False,False]],[[ True, True, True], [ True, True, True], [ True, True, True]],[[False,False,False], [False,False,False], [False,False,False]]])
|
||||
W = 1.79
|
||||
elif id==8:
|
||||
elif stype==8:
|
||||
S = np.array([[[False,False,False], [False,False,False], [False,False,False]],[[False,False,False], [ True, True, True], [False,False,False]],[[False,False,False], [False,False,False], [False,False,False]]])
|
||||
W = 2.68
|
||||
elif id==9:
|
||||
elif stype==9:
|
||||
S = np.array([[[False,False,False], [False,False,False], [False,False,False]],[[False,False,False], [False, True,False], [False,False,False]],[[False,False,False], [False,False,False], [False,False,False]]])
|
||||
W = 4.08
|
||||
else:
|
||||
raise ValueError('Invalid structure type.')
|
||||
return (np.unique(BinaryFieldNd.rotations48(S).reshape(-1,3,3,3),axis=0),W)
|
||||
|
||||
@staticmethod
|
||||
|
|
|
|||
Loading…
Reference in New Issue