field.py

@@ -587,56 +587,18 @@ 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):
-        assert isinstance(data,np.ndarray),\
-               "'data' must be a numpy array."
-        self._dim  = data.shape
-        self._ndim = data.ndim
-        if invert:
-            self.data = data<threshold
-        else:
-            self.data = data>=threshold
-
-    @classmethod
-    def from_field(cls,fld3d,threshold,invert=False):
-        return cls(fld3d.data,threshold,invert=invert)
-
-    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."
-        from scipy import ndimage
-        binarr = ndimage.binary_fill_holes(self.data)
-        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 probe(self,idx):
-        '''Returns whether or not point at index is inside threshold region or not.'''
-        return self.data[tuple(idx)]
-
-    def volume(self):
-        '''Returns volume of region above threshold.'''
-        return np.sum(self.data)
-
 class ConnectedRegions:
-    def __init__(self,binarr,periodicity,connect_diagonals=False,fill_holes=False,bytes_label=32):
+    def __init__(self,boolarr,periodicity,connect_diagonals=False,bytes_label=32):
-        assert isinstance(binarr,np.ndarray) and binarr.dtype==np.dtype('bool'),\
+        assert isinstance(boolarr,np.ndarray) and boolarr.dtype==np.dtype('bool'),\
-               "'binarr' must be a numpy array of dtype('bool')."
+               "'boolarr' must be a numpy array of dtype('bool')."
         assert all([isinstance(x,(bool,int)) for x in periodicity]),\
                 "'periodicity' requires bool values."
         assert bytes_label in (8,16,32,64),\
                "'bytes_label' must be one of {8,16,32,64}."
-        self._dim  = binarr.shape
+        self._dim  = boolarr.shape
-        self._ndim = binarr.ndim
+        self._ndim = boolarr.ndim
         assert self._ndim in (2,3),\
-               "'binarr' must be either two or three dimensional."
+               "'boolarr' must be either two or three dimensional."
         assert len(periodicity)==self._ndim,\
                "Length of 'periodicity' must match number of dimensions of data."
         from scipy import ndimage
@@ -663,7 +625,7 @@ class ConnectedRegions:
         # this does not take into account periodic wrapping
         dtype_label = np.dtype('uint'+str(bytes_label))
         self.label = np.empty(self._dim,dtype=dtype_label)
-        ndimage.label(binarr,structure=connectivity,output=self.label)
+        ndimage.label(boolarr,structure=connectivity,output=self.label)
         self.count = np.max(self.label)
         # Merge labels if there are periodic overlaps
         map_tgt = np.array(range(0,self.count+1),dtype=dtype_label)
@@ -673,9 +635,9 @@ class ConnectedRegions:
             # Merge the first and last plane and compute connectivity
             sl = self._ndim*[slice(None)]
             sl[axis] = (-1,0)
-            binarr_  = binarr[tuple(sl)]
+            boolarr_ = boolarr[tuple(sl)]
-            label_   = np.empty(binarr_.shape,dtype=dtype_label)
+            label_   = np.empty(boolarr_.shape,dtype=dtype_label)
-            ndimage.label(binarr_,structure=connectivity,output=label_)
+            ndimage.label(boolarr_,structure=connectivity,output=label_)
             for val_ in np.unique(label_):
                 # Get all global labels which are associated to a region 
                 # connected over the boundary
@@ -697,18 +659,13 @@ class ConnectedRegions:
         self.count = np.max(map_tgt)
 
     @classmethod
-    def from_field(cls,fld3d,threshold,periodicity,connect_diagonals=False,bytes_label=32,invert=False):
+    def from_field(cls,fld3d,val,periodicity,connect_diagonals=False,bytes_label=32,invert_threshold=False):
-        voxthr = VoxelThreshold.from_field(fld3d,threshold,invert=invert)
+        if invert_threshold:
-        return cls.from_voxelthresh(voxthr,periodicity,
+            return cls(fld3d.data<val,periodicity,
-                            connect_diagonals=connect_diagonals,
+                       connect_diagonals=connect_diagonals,bytes_label=bytes_label)
-                            bytes_label=bytes_label)
+        else:
-
+            return cls(fld3d.data>=val,periodicity,
-    @classmethod
+                       connect_diagonals=connect_diagonals,bytes_label=bytes_label)
-    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
@@ -745,10 +702,6 @@ class ConnectedRegions:
         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."

particle.py

@@ -123,52 +123,6 @@ class Particles:
                 key = ('x','y','z')[axis]
                 self.attr[key] %= self.period[axis]
         return
-    def position_duplicates(self,ipart,padding=0.0):
-        pos  = np.array((px[ipart],py[ipart],pz[ipart]))
-        rp   = pr[ipart]+padding
-        posd = [pos.copy()]
-        for axis in range(3):
-            if self.period[axis] is not None:
-                nposd = len(posd)
-                if pos[axis]-rp<0.0:
-                    for ii in range(nposd):
-                        tmp = posd[ii].copy()
-                        tmp[axis] = np.mod(tmp[axis]-rp,self.period[axis])
-                        posd.append(tmp)
-                if pos[axis]+rp>self.period[axis]:
-                    for ii in range(nposd):
-                        tmp = posd[ii].copy()
-                        tmp[axis] = np.mod(tmp[axis]+rp,self.period[axis])
-                        posd.append(tmp)
-        return posd
-    def mask_field(self,fld,cval=np.nan,padding=0.0):
-      '''Fills grid points which lie inside of solid phase with values.'''
-      for ipart in range(0,part.num):
-        # Slice a box from the field around the particle
-        xp,yp,zp,rp = (px[ipart],py[ipart],pz[ipart],pr[ipart])
-        rp += padding
-        idxlo = np.array(fld.nearest_gridpoint(np.array(xp,yp,zp)-rp,lower=True))
-        idxhi = idxlo+2*rp/fld.spacing
-        
-
-        # Get bounding box of particle
-        idx_x = np.nonzero((xg>=xp-rp) & (xg<=xp+rp))[0]
-        idx_y = np.nonzero((yg>=yp-rp) & (yg<=yp+rp))[0]
-        idx_z = np.nonzero((zg>=zp-rp) & (zg<=zp+rp))[0]
-        # Triple for loop
-        for ii in range(idx_x[0],idx_x[-1]+1):
-          Dx = xg[ii]-xp
-          for jj in range(idx_y[0],idx_y[-1]+1):
-            Dy = yg[jj]-yp
-            for kk in range(idx_z[0],idx_z[-1]+1):
-              Dz = zg[kk]-zp
-              isInside = Dx*Dx+Dy*Dy+Dz*Dz <= rp*rp
-              if isInside:
-                if reconstruct:
-                  self.field[key][ii,jj,kk] = coeff_lin + coeff_rotx*Dx + coeff_roty*Dy + coeff_rotz*Dz
-                else:
-                  self.field[key][ii,jj,kk] = cval
-        return
     def to_vtk(self,deep=False):
         import pyvista as pv
         position = np.vstack([self.attr[key] for key in ('x','y','z')]).transpose()
@@ -176,6 +130,7 @@ class Particles:
         for key in self.attr:
             mesh[key] = self.attr[key]
         return mesh
+
     def glyph(self,theta_resolution=30,phi_resolution=30,deep=False):
         import pyvista as pv
         assert self.has_attribute('r'), "Attribute 'r' required."