From d7658ea6597ce75bda1595445ed1cff612d64e25 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Sat, 7 Aug 2021 01:12:25 +0200 Subject: [PATCH] added structure elements and weights for voxel sruface area estimation --- field.py | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 64 insertions(+), 4 deletions(-) diff --git a/field.py b/field.py index ea46c0d..1007b46 100644 --- a/field.py +++ b/field.py @@ -889,14 +889,17 @@ class BinaryFieldNd: data_,scal_ = self.isolate_feature(lab_,array=array) data_,scal_ = data_[0],scal_[0] # Fill interior in case we filled holes which is not in scal_ - data_ = ndimage.binary_erosion(data_,structure=self.structure,iterations=2) - scal_[data_] = 1.0 - scal_ = np.pad(scal_,pw,mode='reflect',reflect_type='odd') + li = ndimage.binary_erosion(data_,structure=self.structure,iterations=2) + scal_[li] = 1.0 + li = np.logical_not(ndimage.binary_dilation(data_,structure=self.structure,iterations=2)) + scal_[li] = -1.0 + # scal_ = np.pad(scal_,pw,mode='reflect',reflect_type='odd') + scal_ = np.pad(scal_,pw,mode='constant',constant_values=-1.0) pd = Field3d(scal_,origin,spacing).vtk_contour(0.0) else: data_ = self.isolate_feature(lab_)[0] data_ = np.pad(data_.astype(float),pw,mode='constant',constant_values=-1.0) - pd = Field3d(data_,origin,spacing).vtk_contour(0.5).smooth(1000) + pd = Field3d(data_,origin,spacing).vtk_contour(0.5) output.append(pd) return output @@ -919,6 +922,63 @@ class BinaryFieldNd: else: raise ValueError('Invalid input. Accepting int,list,tuple,ndarray.') + @staticmethod + def surface_voxel_structure(id): + '''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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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 + return (np.unique(BinaryFieldNd.rotations48(S).reshape(-1,3,3,3),axis=0),W) + + @staticmethod + def rotations48(a): + '''Generates all possible permutations for numpy array. + Code taken from + https://stackoverflow.com/questions/33190042/how-to-calculate-all-24-rotations-of-3d-array + ''' + import itertools + # Get all combinations of axes that are permutable + n = a.ndim + axcomb = np.array(list(itertools.permutations(range(n), n))) + # Initialize output array + out = np.zeros((6,2,2,2,) + a.shape,dtype=a.dtype) + # Run loop through all axes for flipping and permuting each axis + for i,ax in enumerate(axcomb): + for j,fx in enumerate([1,-1]): + for k,fy in enumerate([1,-1]): + for l,fz in enumerate([1,-1]): + out[i,j,k,l] = np.transpose(a[::fx,::fy,::fz],ax) + return out + class ChunkIterator: '''Iterates through all chunks. 'snapshot' must be an instance of a class which returns a Field3d from the method call