implementing particle masking, unfinished

field.py

@@ -587,18 +587,56 @@ 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,boolarr,periodicity,connect_diagonals=False,bytes_label=32):
+    def __init__(self,binarr,periodicity,connect_diagonals=False,fill_holes=False,bytes_label=32):
-        assert isinstance(boolarr,np.ndarray) and boolarr.dtype==np.dtype('bool'),\
+        assert isinstance(binarr,np.ndarray) and binarr.dtype==np.dtype('bool'),\
-               "'boolarr' must be a numpy array of dtype('bool')."
+               "'binarr' 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  = boolarr.shape
+        self._dim  = binarr.shape
-        self._ndim = boolarr.ndim
+        self._ndim = binarr.ndim
         assert self._ndim in (2,3),\
-               "'boolarr' must be either two or three dimensional."
+               "'binarr' 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
@@ -625,7 +663,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(boolarr,structure=connectivity,output=self.label)
+        ndimage.label(binarr,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)
@@ -635,9 +673,9 @@ class ConnectedRegions:
             # Merge the first and last plane and compute connectivity
             sl = self._ndim*[slice(None)]
             sl[axis] = (-1,0)
-            boolarr_ = boolarr[tuple(sl)]
+            binarr_  = binarr[tuple(sl)]
-            label_   = np.empty(boolarr_.shape,dtype=dtype_label)
+            label_   = np.empty(binarr_.shape,dtype=dtype_label)
-            ndimage.label(boolarr_,structure=connectivity,output=label_)
+            ndimage.label(binarr_,structure=connectivity,output=label_)
             for val_ in np.unique(label_):
                 # Get all global labels which are associated to a region 
                 # connected over the boundary
@@ -659,13 +697,18 @@ class ConnectedRegions:
         self.count = np.max(map_tgt)
 
     @classmethod
-    def from_field(cls,fld3d,val,periodicity,connect_diagonals=False,bytes_label=32,invert_threshold=False):
+    def from_field(cls,fld3d,threshold,periodicity,connect_diagonals=False,bytes_label=32,invert=False):
-        if invert_threshold:
+        voxthr = VoxelThreshold.from_field(fld3d,threshold,invert=invert)
-            return cls(fld3d.data<val,periodicity,
+        return cls.from_voxelthresh(voxthr,periodicity,
-                       connect_diagonals=connect_diagonals,bytes_label=bytes_label)
+                            connect_diagonals=connect_diagonals,
-        else:
+                            bytes_label=bytes_label)
-            return cls(fld3d.data>=val,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
@@ -702,6 +745,10 @@ 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,6 +123,52 @@ 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()
@@ -130,7 +176,6 @@ 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."