field.py

@@ -603,260 +603,44 @@ 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):
-class BinaryFieldNd:
+        assert isinstance(data,np.ndarray),\
-    def __init__(self,input):
+               "'data' must be a numpy array."
-        assert isinstance(input,np.ndarray) and input.dtype==np.dtype('bool'),\
+        self._dim  = data.shape
-               "'input' must be a numpy array of dtype('bool')."
+        self._ndim = data.ndim
-        self.data        = input
+        if invert:
-        self._dim        = input.shape
+            self.data = data<threshold
-        self._ndim       = input.ndim
+        else:
-        self.labels      = None
+            self.data = data>=threshold
-        self.nlabels     = 0
-        self.wrap        = tuple(self._ndim*[None])
-        self._feat_slice = None
-        self.set_structure(False)
-        self.set_periodicity(self._ndim*[False])
 
     @classmethod
-    def from_threshold(cls,fld,threshold,invert=False):
+    def from_field(cls,fld3d,threshold,invert=False):
-        if isinstance(fld,Field3d):
+        return cls(fld3d.data,threshold,invert=invert)
-            fld = fld.data
-        if invert:
-            return cls(fld<threshold)
-        else:
-            return cls(fld>=threshold)
 
-    def set_periodicity(self,periodicity):
+    def fill_holes(self,periodicity=(False,False,False)):
+        '''Fills topological holes in threshold regions.'''
         assert all([isinstance(x,(bool,int)) for x in periodicity]),\
                 "'periodicity' requires bool values."
-        assert len(periodicity)==self._ndim,\
+        from scipy import ndimage
-                "Number of entries in 'periodicity' must match dimension of binary field."
+        binarr = ndimage.binary_fill_holes(self.data)
-        self.periodicity = tuple(bool(x) for x in periodicity)
+        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
         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)
-        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.
-           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))
-                sl_pre = tuple(slice(-2,-1)   if ii==axis else slice(None) for ii in range(self._ndim))
-                lab_lo = labels_[sl_lo]
-                lab_hi = labels_[sl_hi]
-                lab_pre = np.unique(labels_[sl_pre]) # all labels in last (unwrapped) slice
-                # 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_
-                        if source_ in lab_pre: # check if source is not a ghost
-                            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_)
-        assert nlabels_==len(idx_)-1, "DEBUG assertion"
-        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)
-        # If labels have been computed already, recompute them to stay consistent
-        if self.labels is not None:
-            self.label()
-
     def probe(self,idx):
-        '''Returns whether or not a point at idx is True or False.'''
+        '''Returns whether or not point at index is inside threshold region or not.'''
         return self.data[tuple(idx)]
 
     def volume(self):
-        '''Returns the sum of True values.'''
+        '''Returns volume of region above threshold.'''
         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())
-            # 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
-        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(self.nlabels+1),dtype=self.labels.dtype)
-        map_[selection] = 0
-        # Remove gaps from target mapping
-        idx_,map_ = np.unique(map_,return_index=True,return_inverse=True)[1:3]
-        # Discard regions
-        self.labels  = map_[self.labels]
-        self.nlabels = np.max(map_)
-        for axis in range(self._ndim):
-            if self.wrap[axis] is not None:
-                self.wrap[axis] = self.wrap[axis][idx_]
-        self.data = self.labels>0
-
-    def isolate_feature(self,selection):
-        from scipy import ndimage
-        if self.labels is None:
-            self.label()
-        selection = self._select_feature(selection)
-        output = []
-        for lab_ in selection:
-            # Extract feature of interest
-            if lab_==0:
-                data_ = np.logical_not(self.data)
-            else:
-                data_ = (self.labels[self._feat_slice[lab_-1]]==lab_)
-            # If feature is wrapped periodically, duplicate it and extract
-            # largest one
-            iswrapped = False
-            rep_ = self._ndim*[1]
-            for axis in range(self._ndim):
-                if self.wrap[axis] is not None and self.wrap[axis][lab_]:
-                    rep_[axis] = 2
-                    iswrapped = True
-            if iswrapped:
-                data_  = np.tile(data_,rep_)
-                l_,nl_ = ndimage.label(data_,structure=self.structure)
-                vol_   = np.bincount(l_.ravel())
-                il_    = np.argmax(vol_[1:])+1
-                sl_    = ndimage.find_objects(l_==il_)[0]
-                data_  = data_[sl_]
-            # Add to output
-            output.append(data_)
-        return output
-
-    def triangulate_feature(self,selection,origin=(0,0,0),spacing=(1,1,1)):
-        assert self._ndim==3, "Triangulation requires 3D data."
-        if self.labels is None:
-            self.label()
-        selection = self._select_feature(selection)
-        output = []
-        pw = tuple((1,1) for ii in range(self._ndim))
-        for lab_ in selection:
-            data_ = np.pad(self.isolate_feature(lab_)[0],pw,mode='constant')
-            output.append(Field3d(data_,origin,spacing).vtk_contour(0.5))
-        return output
-
-    def _select_feature(self,selection):
-        dtype = self.labels.dtype
-        if selection is None:
-            return np.array(range(1,self.nlabels+1))
-        elif np.issubdtype(type(selection),np.integer):
-            return np.array(selection,dtype=dtype,ndmin=1)
-        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.nlabels+1,\
-                        "Boolean indexing must provide count+1 values."
-            else:
-                selection = selection.astype(dtype)
-                assert np.max(selection)<=self.nlabels 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'),\