Touch ups to ibip, ibop, and qbop to get all tests passing

src/porespy/filters/__init__.py

@@ -5,7 +5,7 @@
 
 This module contains a variety of functions for altering images based on
 the structural characteristics, such as pore sizes.  A definition of a
-*filter* is a function that returns an image the shape as the original
+*filter* is a function that returns an image the same shape as the original
 image, but with altered values.
 
 .. currentmodule:: porespy

src/porespy/filters/_size_seq_satn.py

@@ -10,6 +10,7 @@
     'pc_to_satn',
     'pc_to_seq',
     'satn_to_seq',
+    'size_to_pc',
 ]
 
 
@@ -403,3 +404,88 @@ def satn_to_seq(satn, im=None, mode='drainage'):
         seq[~im] = 0
     seq[uninvaded] = -1
     return seq
+
+
+def size_to_pc(im, size, f=None, **kwargs):
+    r"""
+    Converts size map into capillary pressure map
+
+    Parameters
+    ----------
+    im : ndarray
+        A Numpy array with ``True`` values indicating the void space.
+    size : ndarray
+        The image containing invasion size values in each voxel. Solid
+        should be indicated as 0's and uninvaded voxels as -1.
+    f : function handle, optional
+        A function to compute the capillary pressure which receives `size` as
+        the first argument, followed by any additional `**kwargs`. If not
+        provided then the Washburn equation is used, which requires `theta` and
+        `sigma` to be specified as `kwargs`.
+    **kwargs : Key word arguments
+        All additional keyword arguments are passed on to `f`.
+
+    Returns
+    -------
+    pc : ndarray
+        An image with each voxel containing the capillary pressure at which it was
+        invaded. Any uninvaded voxels in `size` are set to `np.inf` which is meant
+        to indicate that these voxels are never invaded.
+
+    Notes
+    -----
+    The function `f` should be of the form:
+
+    .. code-block::
+
+        def func(r, a, b):
+            pc = ...  # Some equation for capillary pressure using r, a and b
+            return pc
+
+    """
+    if f is None:
+        def f(r, sigma, theta, voxel_size):
+            pc = -2*sigma*np.cos(np.deg2rad(theta))/(r*voxel_size)
+            return pc
+    pc = f(size, **kwargs)
+    pc[~im] = 0.0
+    pc[im == -1] = np.inf
+    return pc
+
+
+# def satn_to_time(im, satn, flow_rate):
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

src/porespy/metrics/_funcs.py

@@ -958,8 +958,7 @@ def phase_fraction(im, normed=True):
     return results
 
 
-def pc_curve(im, sizes=None, pc=None, seq=None,
-             sigma=1.0, theta=180, voxel_size=1.0):
+def pc_curve(im, pc, seq=None):
     r"""
     Produces a Pc-Snwp curve given a map of meniscus radii or capillary
     pressures at which each voxel was invaded
@@ -969,34 +968,14 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
     im : ndarray
         The voxel image of the porous media with ``True`` values indicating
         the void space
-    sizes : ndarray, optional
-        An image containing the sphere radii at which each voxel was invaded
-        during an invasion experiment. The curve is generated by scanning from
-        largest to smallest values and computing the corresponding saturation.
-    pc : ndarray, optional
+    pc : ndarray
         An image containing the capillary pressures at which each voxel was
-        invaded during an invasion experiment. The curve is generated by
-        scanning from smallest to larget values and computing the
-        corresponding saturation.
+        invaded during an invasion experiment. This image can be produced
+        using `size_to_pc` if not available.
     seq : ndarray, optional
         An image containing invasion sequence values, such as that returned
         from the ``ibip`` function. The curve is generated by scanning from
         lowest to highest values and computing the corresponding saturation.
-    sigma : float, optional
-        The surface tension of the fluid-fluid system of interest.
-        This argument is ignored if ``pc`` are specified, otherwise it
-        is used in the Washburn equation to convert ``sizes`` to capillary
-        pc.
-    theta : float
-        The contact angle measured through the invading phase in degrees.
-        This argument is ignored if ``pc`` are specified, otherwise it
-        is used in the Washburn equation to convert ``sizes`` to capillary
-        pressures.
-    voxel_size : float
-        The voxel resolution of the image.
-        This argument is ignored if ``pc`` are specified, otherwise it
-        is used in the Washburn equation to convert ``sizes`` to capillary
-        pressures.
 
     Returns
     -------
@@ -1013,18 +992,6 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
                             phase at each pressure in ``pc``
         ==================  ===================================================
 
-    Notes
-    -----
-    If ``sizes`` is provided, then the Washburn equation is used to convert
-    the radii to capillary pressures, using the given ``sigma`` and ``theta``
-    values, along with the ``voxel_size`` if the values are in voxel radii.
-    For more control over how capillary pressure model, it can be computed by
-    hand, for example:
-
-        $$ pc = \frac{-2*0.072*np.cos(np.deg2rad(180))}{sizes \cdot voxel_size} $$
-
-    then passed in as the ``pc`` argument.
-
     Examples
     --------
     `Click here
@@ -1033,66 +1000,23 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
 
     """
     tqdm = get_tqdm()
-    if seq is not None:
-        seqs = np.unique(seq)
-        seqs = seqs[seqs > 0]
-        x = []
-        y = []
-        with tqdm(seqs, **settings.tqdm) as pbar:
-            for n in seqs:
-                pbar.update()
-                mask = seq == n
-                if (pc is not None) and (sizes is not None):
-                    raise Exception("Only one of pc or sizes can be specified")
-                elif pc is not None:
-                    pressure = pc[mask][0]
-                elif sizes is not None:
-                    r = sizes[mask][0]*voxel_size
-                    pressure = -2*sigma*np.cos(np.deg2rad(theta))/r
-                x.append(pressure)
-                snwp = ((seq <= n)*(seq > 0) *
-                        (im == 1)).sum(dtype=np.int64)/im.sum(dtype=np.int64)
-                y.append(snwp)
-        pc_curve = Results()
-        pc_curve.pc = x
-        pc_curve.snwp = y
-    elif sizes is not None:
-        if im is None:
-            im = ~(sizes == 0)
-        sz = np.unique(sizes)[:0:-1]
-        sz = np.hstack((sz[0]*2, sz))
-        x = []
-        y = []
-        with tqdm(sz, **settings.tqdm) as pbar:
-            for n in sz:
-                pbar.update()
-                r = n*voxel_size
-                pc = -2*sigma*np.cos(np.deg2rad(theta))/r
-                x.append(pc)
-                snwp = ((sizes >= n)*(im == 1)).sum(dtype=np.int64) \
-                    / im.sum(dtype=np.int64)
-                y.append(snwp)
-        pc_curve = Results()
-        pc_curve.pc = x
-        pc_curve.snwp = y
-    elif pc is not None:
-        Ps = np.unique(pc[im])
-        # Utilize the fact that -inf and +inf will be at locations 0 & -1 in Ps
-        if Ps[-1] == np.inf:
-            Ps[-1] = Ps[-2]*2
-        if Ps[0] == -np.inf:
-            Ps[0] = Ps[1] - np.abs(Ps[1]/2)
-        else:
-            # Add a point at begining to ensure curve starts a 0, if no residual
-            Ps = np.hstack((Ps[0] - np.abs(Ps[0]/2), Ps))
-        y = []
-        Vp = im.sum(dtype=np.int64)
-        temp = pc[im]
-        for p in tqdm(Ps, **settings.tqdm):
-            y.append((temp <= p).sum(dtype=np.int64)/Vp)
-        pc_curve = Results()
-        pc_curve.pc = Ps
-        pc_curve.snwp = np.array(y)
+    Ps = np.unique(pc[im])
+    # Utilize the fact that -inf and +inf will be at locations 0 & -1 in Ps
+    if Ps[-1] == np.inf:
+        Ps[-1] = Ps[-2]*2
+    if Ps[0] == -np.inf:
+        Ps[0] = Ps[1] - np.abs(Ps[1]/2)
+    else:
+        # Add a point at begining to ensure curve starts a 0, if no residual
+        Ps = np.hstack((Ps[0] - np.abs(Ps[0]/2), Ps))
+    y = []
+    Vp = im.sum(dtype=np.int64)
+    temp = pc[im]
+    for p in tqdm(Ps, **settings.tqdm):
+        y.append((temp <= p).sum(dtype=np.int64)/Vp)
+    pc_curve = Results()
+    pc_curve.pc = Ps
+    pc_curve.snwp = np.array(y)
     return pc_curve
 
 

src/porespy/simulations/_ibip.py

@@ -105,8 +105,8 @@ def ibip(
     if dt is None:  # Find dt if not given
         dt = edt(im)
     # Initialize inv image with -1 in the solid, and 0's in the void
-    inv = -1*(~im)
-    sizes = -1*(~im)
+    seq = -1*(~im)
+    sizes = -1.0*(~im)
     scratch = np.copy(bd)
     for step in tqdm(range(1, maxiter), **settings.tqdm):
         # Find insertion points
@@ -119,8 +119,8 @@ def ibip(
         dt_thresh = dt >= r_max
         # Extract the actual coordinates of the insertion sites
         pt = _where(edge*dt_thresh)
-        inv = _insert_disk_at_points(
-            im=inv,
+        seq = _insert_disk_at_points(
+            im=seq,
             coords=pt,
             r=int(r_max),
             v=step,
@@ -131,7 +131,7 @@ def ibip(
                 im=sizes,
                 coords=pt,
                 r=int(r_max),
-                v=int(r_max),
+                v=r_max,
                 smooth=True,
             )
         dt, bd = _update_dt_and_bd(dt, bd, pt)
@@ -143,18 +143,18 @@ def ibip(
             smooth=False,
         )
     # Convert inv image so that uninvaded voxels are set to -1 and solid to 0
-    temp = inv == 0  # Uninvaded voxels are set to -1 after _ibip
-    inv[~im] = 0
-    inv[temp] = -1
-    inv = make_contiguous(im=inv, mode='symmetric')
+    temp = seq == 0  # Uninvaded voxels are set to -1 after _ibip
+    seq[~im] = 0
+    seq[temp] = -1
+    seq = make_contiguous(im=seq, mode='symmetric')
     # Deal with invasion sizes similarly
     temp = sizes == 0
     sizes[~im] = 0
     sizes[temp] = -1
     results = Results()
     results.im_size = np.copy(sizes)
-    results.im_seq = np.copy(inv)
-    results.im_satn = seq_to_satn(inv)
+    results.im_seq = np.copy(seq)
+    results.im_satn = seq_to_satn(seq)
     return results