fixed false usage of _dim -> dim() method should be called. UNTESTED

field.py

@@ -128,7 +128,7 @@ class Field3d:
         return cls(data,origin,spacing)
 
     @classmethod
-    def allocate(cls,dim,origin,spacing,fill=None,dtype=np.float64,pseudo=False):
+    def allocate(cls,dim,origin,spacing,fill=None,dtype=np.float64):
         '''Allocates an empty field, or a field filled with 'fill'.'''
         assert isinstance(dim,(tuple,list,np.ndarray)) and len(dim)==3,\
             "'dim' must be a tuple/list of length 3."
@@ -136,17 +136,25 @@ class Field3d:
             "'origin' must be a tuple/list of length 3."
         assert isinstance(spacing,(tuple,list,np.ndarray)) and len(spacing)==3,\
             "'spacing' must be a tuple/list of length 3."
-        if pseudo:
-            data = np.empty((0,0,0))
-            r = cls(data,origin,spacing)
-            r._dim = dim
-            return r
+        if fill is None:
+            data = np.empty(dim,dtype=dtype)
         else:
-            if fill is None:
-                data = np.empty(dim,dtype=dtype)
-            else:
-                data = np.full(dim,fill,dtype=dtype)
-            return cls(data,origin,spacing)
+            data = np.full(dim,fill,dtype=dtype)
+        return cls(data,origin,spacing)
+    
+    @classmethod
+    def pseudo_field(cls,dim,origin,spacing):
+        '''Creates a Field3d instance without allocating any memory.'''
+        assert isinstance(dim,(tuple,list,np.ndarray)) and len(dim)==3,\
+            "'dim' must be a tuple/list of length 3."
+        assert isinstance(origin,(tuple,list,np.ndarray)) and len(origin)==3,\
+            "'origin' must be a tuple/list of length 3."
+        assert isinstance(spacing,(tuple,list,np.ndarray)) and len(spacing)==3,\
+            "'spacing' must be a tuple/list of length 3."
+        data = np.empty((0,0,0))
+        r = cls(data,origin,spacing)
+        r._dim = dim
+        return r
 
     def save(self,file,name='Field3d',truncate=False):
         import h5py
@@ -163,6 +171,9 @@ class Field3d:
     def copy(self):
         return Field3d(self.data.copy(),self.origin,self.spacing)
 
+    def pseudo_copy(self):
+        return self.pseudo_field(self.dim(),self.origin,self.spacing)
+
     def insert_subfield(self,subfield):
         assert all([abs(subfield.spacing[ii]-self.spacing[ii])<self.eps_collapse 
                         for ii in range(3)]), "spacing differs. Got {}, have {}".format(subfield.spacing,self.spacing)
@@ -638,15 +649,39 @@ class ConnectedRegions:
         self.label = map_tgt[self.label]
         self.count = np.max(map_tgt)
 
+    @classmethod
+    def from_field(cls,fld3d,val,periodicity,connect_diagonals=False,bytes_label=32,invert_threshold=False):
+        if invert_threshold:
+            return cls(fld3d.data<val,periodicity,
+                       connect_diagonals=connect_diagonals,bytes_label=bytes_label)
+        else:
+            return cls(fld3d.data>=val,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
@@ -658,9 +693,30 @@ class ConnectedRegions:
         self.label = map_tgt[self.label]
         self.count = np.max(map_tgt)
 
+    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 isinstance(selection,int):
+        if np.issubdtype(type(selection),np.integer):
             return np.array(selection,dtype=dtype)
         elif isinstance(selection,(list,tuple,np.ndarray)):
             selection = np.array(selection)
@@ -673,7 +729,7 @@ class ConnectedRegions:
                         "Entry in selection is out-of-bounds."
             return selection
         else:
-            raise ValueError('Invalid input.')
+            raise ValueError('Invalid input. Accepting int,list,tuple,ndarray.')
 
 
 class ChunkIterator: