Repository: PMEAL/porespy

Issue:
Add a `pc_to_seq` function
Not sure why we don't have this already.  The only complication is that it does not work with ibip results since the Pcs are not filled sequentially.  Not sure how to inidcate this.

Discussion:

Patch:
diff --git a/porespy/filters/_size_seq_satn.py b/porespy/filters/_size_seq_satn.py
index 6b83bc26d..ace920e73 100644
--- a/porespy/filters/_size_seq_satn.py
+++ b/porespy/filters/_size_seq_satn.py
@@ -8,6 +8,7 @@ __all__ = [
     'size_to_satn',
     'seq_to_satn',
     'pc_to_satn',
+    'pc_to_seq',
     'satn_to_seq',
 ]
 
@@ -220,6 +221,72 @@ def seq_to_satn(seq, im=None, mode='drainage'):
     return satn
 
 
+def pc_to_seq(pc, im, mode='drainage'):
+    r"""
+    Converts an image of capillary entry pressures to invasion sequence values
+
+    Parameters
+    ----------
+    pc : ndarray
+        A Numpy array with the value in each voxel indicating the capillary
+        pressure at which it was invaded. In order to accommodate the
+        possibility of both positive and negative capillary pressure values,
+        uninvaded voxels should be indicated by ``+inf`` and residual phase
+        by ``-inf``. Solid vs void phase is defined by ``im`` which is
+        mandatory.
+    im : ndarray
+        A Numpy array with ``True`` values indicating the void space
+    mode : str
+        Controls how the pressures are converted to sequence. The options are:
+
+        ============= ==============================================================
+        `mode`        Description
+        ============= ==============================================================
+        'drainage'    The pressures are assumed to have been filled from smallest to
+                      largest. Voxels with -np.inf are treated as though they are
+                      invaded by non-wetting fluid at the start of the process, and
+                      voxels with +np.inf are treated as though they are never
+                      invaded.
+        'imbibition'  The pressures are assumed to have been filled from largest to
+                      smallest. Voxels with -np.inf are treated as though they are
+                      already occupied by non-wetting fluid at the start of the
+                      process, and voxels with +np.inf are treated as though they
+                      are filled with wetting phase.
+        ============= ==============================================================
+
+    Returns
+    -------
+    seq : ndarray
+        A Numpy array the same shape as `pc`, with each voxel value indicating
+        the sequence at which it was invaded, according to the specified `mode`.
+        Uninvaded voxels are set to -1.
+
+    Notes
+    -----
+    Voxels with `+inf` are treated as though they were never invaded so are given a
+    sequence value of -1. Voxels with  `-inf` are treated as though they were
+    invaded by non-wetting phase at the start of the simulation so are given a
+    sequence number of 1 for both mode `drainage` and `imbibition`.
+
+    Examples
+    --------
+    `Click here
+    <https://porespy.org/examples/filters/reference/pc_to_seq.html>`_
+    to view online example.
+    """
+    inf = pc == np.inf  # save for later
+    if mode == 'drainage':
+        bins = np.unique(pc)
+    elif mode == 'imbibition':
+        pc[pc == -np.inf] = np.inf
+        bins = np.unique(pc)[-1::-1]
+    a = np.digitize(pc, bins=bins)
+    a[~im] = 0
+    a[np.where(inf)] = -1
+    a = make_contiguous(a, mode='symmetric')
+    return a
+
+
 def pc_to_satn(pc, im, mode='drainage'):
     r"""
     Converts an image of capillary entry pressures to saturation values
@@ -242,9 +309,9 @@ def pc_to_satn(pc, im, mode='drainage'):
         `mode`        Description
         ============= ==============================================================
         'drainage'    The pressures are assumed to have been filled from smallest to
-                      largest, ignoring +/- infs
+                      largest.
         'imbibition'  The pressures are assumed to have been filled from largest to
-                      smallest, ignoring +/- infs
+                      smallest
         ============= ==============================================================
 
     Returns
@@ -252,7 +319,9 @@ def pc_to_satn(pc, im, mode='drainage'):
     satn : ndarray
         A Numpy array the same shape as `pc`, with each voxel value indicating
         the global saturation at which it was invaded, according to the specified
-        `mode`.
+        `mode`. Voxels with  `-inf` are treated as though they were invaded
+        at the start of the simulation so are given a sequence number of 1 for both
+        mode `drainage` and `imbibition`.
 
     Notes
     -----
diff --git a/porespy/metrics/_funcs.py b/porespy/metrics/_funcs.py
index e44651570..328cf83f2 100644
--- a/porespy/metrics/_funcs.py
+++ b/porespy/metrics/_funcs.py
@@ -30,6 +30,7 @@ __all__ = [
     "pc_curve",
     "pc_curve_from_ibip",
     "pc_curve_from_mio",
+    "pc_map_to_pc_curve",
 ]
 
 
@@ -1055,10 +1056,14 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
             for n in seqs:
                 pbar.update()
                 mask = seq == n
-                # The following assumes only one size found, which was confirmed
-                r = sizes[mask][0]*voxel_size
-                pc = -2*sigma*np.cos(np.deg2rad(theta))/r
-                x.append(pc)
+                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)
@@ -1104,6 +1109,62 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
     return pc_curve
 
 
+def pc_map_to_pc_curve(pc, im, seq=None):
+    r"""
+    Converts a pc map into a capillary pressure curve
+
+    Parameters
+    ----------
+    pc : ndarray
+        A numpy array with each voxel containing the capillary pressure at which
+        it was invaded. `-inf` indicates voxels which are already filled with
+        non-wetting fluid, and `+inf` indicates voxels that are not invaded by
+        non-wetting fluid (e.g., trapped wetting phase). Solids should be
+        noted by `+inf` but this is also enforced inside the function using `im`.
+    im : ndarray
+        A numpy array with `True` values indicating the void space and `False`
+        elsewhere. This is necessary to define the total void volume of the domain
+        for computing the saturation.
+    seq : ndarray, optional
+        A numpy array with each voxel containing the sequence at which it was
+        invaded. This is required when analyzing results from invasion percolation
+        since the pressures in `pc` do not correspond to the sequence in which
+        they were filled.
+
+    Returns
+    -------
+    results : dataclass-like
+        A dataclass like object with the following attributes:
+
+        ================== =========================================================
+        Attribute          Description
+        ================== =========================================================
+        pc                 The capillary pressure
+        snwp               The fraction of void space filled by non-wetting
+                           phase at each pressure in ``pc``
+        ================== =========================================================
+
+    Notes
+    -----
+    To use this function with the results of `porosimetry` or `ibip` the sizes map
+    must be converted to a capillary pressure map first.  `drainage` and `invasion`
+    both return capillary pressure maps which can be passed directly as `pc`.
+    """
+    pc[~im] = np.inf  # Ensure solid voxels are set to inf invasion pressure
+    if seq is None:
+        pcs, counts = np.unique(pc, return_counts=True)
+    else:
+        vals, index, counts = np.unique(seq, return_index=True, return_counts=True)
+        pcs = pc.flatten()[index]
+    snwp = np.cumsum(counts[pcs < np.inf])/im.sum()
+    pcs = pcs[pcs < np.inf]
+
+    results = Results()
+    results.pc = pcs
+    results.snwp = snwp
+    return results
+
+
 def satn_profile(satn, s=None, im=None, axis=0, span=10, mode='tile'):
     r"""
     Computes a saturation profile from an image of fluid invasion


Test patch:
diff --git a/test/unit/test_filters_size_seq_satn.py b/test/unit/test_filters_size_seq_satn.py
index 2aca65070..b3ea28425 100644
--- a/test/unit/test_filters_size_seq_satn.py
+++ b/test/unit/test_filters_size_seq_satn.py
@@ -249,6 +249,17 @@ class SeqTest():
         assert satn[0, 0] == 0.0
         assert satn[0, 1] == 0.9
 
+    def test_pc_to_seq(self):
+        pc = 10.0*np.tile(np.atleast_2d(np.arange(0, 21)), [21, 1])
+        pc[:, 0] = 0
+        pc[:, -5] = np.inf
+        im = pc > 0
+        seq = ps.filters.pc_to_seq(pc=pc, im=im, mode='drainage')
+        assert seq[0, 0] == 0
+        assert seq[0, 1] == 1
+        assert seq[0, -1] == 19
+        assert seq[0, -5] == -1
+
 
 if __name__ == '__main__':
     t = SeqTest()
diff --git a/test/unit/test_metrics.py b/test/unit/test_metrics.py
index 8879d3846..595f59e10 100644
--- a/test/unit/test_metrics.py
+++ b/test/unit/test_metrics.py
@@ -309,7 +309,72 @@ class MetricsTest():
         satn[:25, :] = 0
         satn[-25:, :] = -1
         with pytest.raises(Exception):
-            prof1 = ps.metrics.satn_profile(satn=satn, s=0.5)
+            _ = ps.metrics.satn_profile(satn=satn, s=0.5)
+
+    def test_pc_map_to_pc_curve_drainage_with_trapping_and_residual(self):
+        vx = 50e-6
+        im = ps.generators.blobs(shape=[200, 200], porosity=0.5, blobiness=2, seed=0)
+        mio = ps.filters.porosimetry(im)
+        trapped = im*(~ps.filters.fill_blind_pores(im))
+        residual = im*(~trapped)*(mio < mio.mean())
+        pc = -2*0.072*np.cos(np.radians(110))/(mio*vx)
+        pc[trapped] = np.inf
+        pc[residual] = -np.inf
+        d = ps.metrics.pc_map_to_pc_curve(pc, im)
+        assert d.snwp[0] == residual.sum()/im.sum()
+        assert d.snwp[-1] == (im.sum() - trapped.sum())/im.sum()
+
+    def test_pc_map_to_pc_curve_invasion_with_trapping(self):
+        vx = 50e-6
+        im = ps.generators.blobs(shape=[200, 200], porosity=0.5, blobiness=2, seed=0)
+        ibip = ps.simulations.ibip(im=im)
+        pc = -2*0.072*np.cos(np.radians(110))/(ibip.inv_sizes*vx)
+        trapped = ibip.inv_sequence == -1
+        # residual = pc*im > 500
+        pc[trapped] = np.inf
+        seq = ibip.inv_sequence
+        d = ps.metrics.pc_map_to_pc_curve(pc=pc, im=im, seq=seq)
+        # assert d.snwp[0] == residual.sum()/im.sum()
+        assert d.snwp[-1] == (im.sum() - trapped.sum())/im.sum()
+
+    def test_pc_map_to_pc_curve_compare_invasion_to_drainage(self):
+        vx = 50e-6
+        im = ps.generators.blobs(shape=[200, 200], porosity=0.6, blobiness=1, seed=0)
+        im = ps.filters.fill_blind_pores(im, conn=8, surface=True)
+
+        # Do drainage without sequence
+        dt = edt(im)
+        mio = ps.filters.porosimetry(im, sizes=np.unique(dt)[1:].astype(int))
+        pc1 = -2*0.072*np.cos(np.radians(110))/(mio*vx)
+        d1 = ps.metrics.pc_map_to_pc_curve(pc=pc1, im=im)
+
+        # Ensure drainage works with sequence
+        seq = ps.filters.pc_to_seq(pc1, im)
+        d3 = ps.metrics.pc_map_to_pc_curve(pc=pc1, im=im, seq=seq)
+
+        # Using the original ibip, which requires that sequence be supplied
+        ibip = ps.simulations.ibip(im=im)
+        pc2 = -2*0.072*np.cos(np.radians(110))/(ibip.inv_sizes*vx)
+        pc2[ibip.inv_sequence < 0] = np.inf
+        seq = ibip.inv_sequence
+        d2 = ps.metrics.pc_map_to_pc_curve(pc=pc2, im=im, seq=seq)
+
+        # Ensure they all return the same Pc values
+        assert_allclose(np.unique(d1.pc), np.unique(d2.pc), rtol=1e-10)
+        assert_allclose(np.unique(d2.pc), np.unique(d3.pc), rtol=1e-10)
+        assert_allclose(np.unique(d1.pc), np.unique(d3.pc), rtol=1e-10)
+
+        # Ensure the high and low saturations are all the same
+        assert d1.snwp[0] == d2.snwp[0]
+        assert d1.snwp[-1] == d2.snwp[-1]
+        assert d2.snwp[0] == d3.snwp[0]
+        assert d2.snwp[-1] == d3.snwp[-1]
+
+        # These graphs should lie perfectly on top of each other
+        # import matplotlib.pyplot as plt
+        # plt.step(d1.pc, d1.snwp, 'r-o', where='post')
+        # plt.step(d3.pc, d3.snwp, 'b--', where='post')
+        # plt.step(d2.pc, d2.snwp, 'g.-', where='post')
 
 
 if __name__ == '__main__':