field.py

@@ -251,9 +251,9 @@ class Field3d:
         '''Returns a binary array indicating which grid points are above,
            or in case of invert=True below, a given threshold.'''
         if invert:
-            return self.data<val
-        else:
             return self.data>=val
+        else:
+            return self.data<val
 
     def coordinate(self,idx,axis=None):
         if axis is None:
@@ -656,10 +656,19 @@ class BinaryFieldNd:
         self.labels      = None
         self.nlabels     = 0
         self.wrap        = tuple(self._ndim*[None])
-        self._featsl     = None
+        self._feat_slice = None
         self.set_structure(False)
         self.set_periodicity(self._ndim*[False])
 
+    @classmethod
+    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 set_periodicity(self,periodicity):
         assert all([isinstance(x,(bool,int)) for x in periodicity]),\
                 "'periodicity' requires bool values."
@@ -686,7 +695,7 @@ class BinaryFieldNd:
             self.labels,self.nlabels,self.wrap = self._labels_periodic()
         else:
             self.labels,self.nlabels = ndimage.label(self.data,structure=self.structure)
-        self._featsl = ndimage.find_objects(self.labels)
+        self._feat_slice = ndimage.find_objects(self.labels)
 
     def _labels_periodic(self,map_to_zero=False):
         '''Label features in an array while taking into account periodic wrapping.
@@ -811,7 +820,6 @@ class BinaryFieldNd:
         return labels
 
     def discard_feature(self,selection):
-        '''Removes a feature from data.'''
         if self.labels is None:
             self.label()
         selection = self._select_feature(selection)
@@ -844,8 +852,8 @@ class BinaryFieldNd:
                 data_ = np.logical_not(self.data)
                 if has_array: data2_ = array
             else:
-                data_ = (self.labels[self._featsl[lab_-1]]==lab_)
+                data_ = (self.labels[self._feat_slice[lab_-1]]==lab_)
-                if has_array: data2_ = array[self._featsl[lab_-1]]
+                if has_array: data2_ = array[self._feat_slice[lab_-1]]
             # If feature is wrapped periodically, duplicate it and extract
             # largest one
             iswrapped = False
@@ -860,7 +868,6 @@ class BinaryFieldNd:
                 vol_   = np.bincount(l_.ravel())
                 il_    = np.argmax(vol_[1:])+1
                 sl_    = ndimage.find_objects(l_==il_)[0]
-                print(sl_)
                 data_  = data_[sl_]
                 if has_array:
                     data2_ = np.tile(data2_,rep_)[sl_]                    
@@ -889,17 +896,14 @@ 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_
-                li = ndimage.binary_erosion(data_,structure=self.structure,iterations=2)
+                data_ = ndimage.binary_erosion(data_,structure=self.structure,iterations=2)
-                scal_[li] = 1.0
+                scal_[data_] = 1.0
-                li = np.logical_not(ndimage.binary_dilation(data_,structure=self.structure,iterations=2))
+                scal_ = np.pad(scal_,pw,mode='reflect',reflect_type='odd')
-                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)
+                pd = Field3d(data_,origin,spacing).vtk_contour(0.5).smooth(1000)
             output.append(pd)
         return output
 
@@ -922,62 +926,168 @@ class BinaryFieldNd:
         else:
             raise ValueError('Invalid input. Accepting int,list,tuple,ndarray.')
 
-    @staticmethod
+class ConnectedRegions:
-    def surface_voxel_structure(id):
+    def __init__(self,binarr,periodicity,connect_diagonals=False,fill_holes=False,bytes_label=32):
-        '''Get structure element and all its permutations for surface voxels
+        assert isinstance(binarr,np.ndarray) and binarr.dtype==np.dtype('bool'),\
-           as defined by 
+               "'binarr' must be a numpy array of dtype('bool')."
-           Mullikin, J. C., & Verbeek, P. W. (1993). Surface area estimation of digitized planes. Bioimaging, 1(1), 6-16.
+        assert all([isinstance(x,(bool,int)) for x in periodicity]),\
-           Also returns weights, where S7-S9 are defined in 
+                "'periodicity' requires bool values."
-           Windreich, G., Kiryati, N., & Lohmann, G. (2003). Voxel-based surface area estimation: from theory to practice. Pattern Recognition, 36(11), 2531-2541.
+        assert bytes_label in (8,16,32,64),\
-        '''
+               "'bytes_label' must be one of {8,16,32,64}."
-        if id==1:
+        self._dim  = binarr.shape
-            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]]])
+        self._ndim = binarr.ndim
-            W = 0.894
+        assert self._ndim in (2,3),\
-        elif id==2:
+               "'binarr' must be either two or three dimensional."
-            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]]])
+        assert len(periodicity)==self._ndim,\
-            W = 1.3409
+               "Length of 'periodicity' must match number of dimensions of data."
-        elif id==3:
+        from scipy import ndimage
-            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]]])
+        # Construct connectivity stencil
-            W = 1.5879
+        if self._ndim==2:
-        elif id==4:
+            connectivity = np.ones((3,3),dtype='bool')
-            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]]])
+            if not connect_diagonals:
-            W = 2.0
+                connectivity[0,0] = False
-        elif id==5:
+                connectivity[0,2] = False
-            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]]])
+                connectivity[2,0] = False
-            W = 8./3.
+                connectivity[0,2] = False
-        elif id==6:
+        else:
-            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]]])
+            connectivity = np.ones((3,3,3),dtype='bool')
-            W = 10./3.
+            if not connect_diagonals:
-        elif id==7:
+                connectivity[0,0,0] = False
-            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]]])
+                connectivity[2,0,0] = False
-            W = 1.79
+                connectivity[0,2,0] = False
-        elif id==8:
+                connectivity[0,0,2] = False
-            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]]])
+                connectivity[2,2,0] = False
-            W = 2.68
+                connectivity[0,2,2] = False
-        elif id==9:
+                connectivity[2,0,2] = False
-            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]]])
+                connectivity[2,2,2] = False
-            W = 4.08
+        # Compute labels:
-        return (np.unique(BinaryFieldNd.rotations48(S).reshape(-1,3,3,3),axis=0),W)
+        # this does not take into account periodic wrapping
-    
+        dtype_label = np.dtype('uint'+str(bytes_label))
-    @staticmethod
+        self.label = np.empty(self._dim,dtype=dtype_label)
-    def rotations48(a):
+        ndimage.label(binarr,structure=connectivity,output=self.label)
-        '''Generates all possible permutations for numpy array.
+        self.count = np.max(self.label)
-           Code taken from
+        # Merge labels if there are periodic overlaps
-           https://stackoverflow.com/questions/33190042/how-to-calculate-all-24-rotations-of-3d-array
+        map_tgt = np.array(range(0,self.count+1),dtype=dtype_label)
-        '''
+        for axis in range(self._ndim):
-        import itertools
+            if not periodicity[axis]:
-        # Get all combinations of axes that are permutable
+                continue
-        n = a.ndim
+            # Merge the first and last plane and compute connectivity
-        axcomb = np.array(list(itertools.permutations(range(n), n)))
+            sl = self._ndim*[slice(None)]
-        # Initialize output array
+            sl[axis] = (-1,0)
-        out = np.zeros((6,2,2,2,) + a.shape,dtype=a.dtype)
+            binarr_  = binarr[tuple(sl)]
-        # Run loop through all axes for flipping and permuting each axis
+            label_   = np.empty(binarr_.shape,dtype=dtype_label)
-        for i,ax in enumerate(axcomb):
+            ndimage.label(binarr_,structure=connectivity,output=label_)
-            for j,fx in enumerate([1,-1]):
+            for val_ in np.unique(label_):
-                for k,fy in enumerate([1,-1]):
+                # Get all global labels which are associated to a region 
-                    for l,fz in enumerate([1,-1]):
+                # connected over the boundary
-                        out[i,j,k,l] = np.transpose(a[::fx,::fy,::fz],ax) 
+                global_labels = list(np.unique(self.label[tuple(sl)][label_==val_]))
-        return out
+                # If there is only one label, nothing needs to be done
+                if len(global_labels)==1:
+                    continue
+                # Determine target label:
+                # this needs to be done recursively because the original
+                # target may already be reassigned
+                tgt = global_labels[0]
+                while tgt!=map_tgt[tgt]:
+                    tgt=map_tgt[tgt]
+                map_tgt[global_labels[1:]] = tgt
+        # Remove gaps from target mapping
+        map_tgt = np.unique(map_tgt,return_inverse=True)[1]
+        # Remap labels
+        self.label = map_tgt[self.label]
+        self.count = np.max(map_tgt)
+
+    @classmethod
+    def from_field(cls,fld3d,threshold,periodicity,connect_diagonals=False,bytes_label=32,invert=False):
+        voxthr = VoxelThreshold.from_field(fld3d,threshold,invert=invert)
+        return cls.from_voxelthresh(voxthr,periodicity,
+                            connect_diagonals=connect_diagonals,
+                            bytes_label=bytes_label)
+
+    @classmethod
+    def from_voxelthresh(cls,voxthr,periodicity,connect_diagonals=False,bytes_label=32):
+        return cls(voxthr.data,periodicity,
+                       connect_diagonals=connect_diagonals,
+                       bytes_label=bytes_label)
+
+
+    def volume(self,label=None):
+        '''Returns volume of labeled regions. 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 label is None:
+            return np.bincount(self.label.ravel())
+        else:
+            return np.sum(self.label==label)
+
+    def volume_domain(self):
+        '''Returns volume of entire domain. Should be equal to sum(volume()).'''
+        return np.prod(self._dim)
+
+    def labels_by_volume(self,descending=True):
+        '''Returns labels of connected regions sorted by volume.'''
+        labels = np.argsort(self.volume()[1:])+1
+        if descending:
+            labels = labels[::-1]
+        return labels
+
+    def discard_regions(self,selection):
+        selection = self._parse_selection(selection)
+        # Map tagged regions to zero in order to discard them
+        map_tgt = np.array(range(0,self.count+1),dtype=self.label.dtype)
+        map_tgt[selection] = 0
+        # Remove gaps from target mapping
+        map_tgt = np.unique(map_tgt,return_inverse=True)[1]
+        # Discard regions
+        self.label = map_tgt[self.label]
+        self.count = np.max(map_tgt)
+
+    def probe(self,idx):
+        '''Returns label for given index.'''
+        return self.label[tuple(idx)]
+
+    def vtk_contour(self,fld3,val,selection):
+        '''Computes contours of a Field3d only within selected structures.'''
+        assert isinstance(fld3,Field3d), "'fld3' must be a Field3d instance."
+        assert tuple(self._dim)==tuple(fld3.dim()), \
+            "'fld3' must be of dimension {}.".format(self._dim)
+        selection = self._parse_selection(selection)
+        from scipy import ndimage
+        # Create binary map of selection
+        map_tgt = np.zeros(self.count+1,dtype='bool')
+        map_tgt[selection] = True
+        binary_map = map_tgt[self.label]
+        # Add an extra cell to get the contour interpolation right
+        print(np.sum(binary_map))
+        binary_map = ndimage.binary_dilation(binary_map)
+        print(np.sum(binary_map))
+        # Extract the subfield
+        fld_con = fld3.copy()
+        fld_con.data[~binary_map] = np.nan
+        # Compute the contour
+        return fld_con.vtk_contour(val)
+
+    def _parse_selection(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 ChunkIterator:
     '''Iterates through all chunks. 'snapshot' must be an instance