Add rescale_network and param.* metadata

src/porespy/networks/_funcs.py

@@ -15,6 +15,7 @@
     "label_phases",
     "label_boundaries",
     "map_to_regions",
+    "rescale_network",
 ]
 
 
@@ -367,3 +368,91 @@ def label_boundaries(
                 network['pore.boundary'] += hits
                 network['pore.' + labels[i][j]] = hits
     return network
+
+
+_LENGTH_NAMES = {
+    'coords', 'local_peak', 'global_peak', 'geometric_centroid', 'perimeter',
+}
+_LENGTH_SUFFIXES = ('_diameter', '_length', '_radius', '_perimeter')
+_AREA_NAMES = {'area'}
+_AREA_SUFFIXES = ('_area',)
+_VOLUME_NAMES = {'volume'}
+_VOLUME_SUFFIXES = ('_volume',)
+_SIZE_FACTOR_SUFFIXES = ('size_factor', 'size_factors')
+
+
+def _scale_exponent(key, ndim):
+    if '.' not in key:
+        return 0
+    name = key.rsplit('.', 1)[1]
+    if name in _LENGTH_NAMES or name.endswith(_LENGTH_SUFFIXES):
+        return 1
+    if name in _AREA_NAMES or name.endswith(_AREA_SUFFIXES):
+        return 2
+    if name in _VOLUME_NAMES or name.endswith(_VOLUME_SUFFIXES):
+        return ndim
+    if name.endswith(_SIZE_FACTOR_SUFFIXES):
+        return 1
+    return 0
+
+
+def rescale_network(network, voxel_size):
+    r"""
+    Rescale a network's geometric properties to a new voxel size.
+
+    Networks extracted by porespy carry their original ``voxel_size`` in
+    ``param.voxel_size``. This function uses that value to compute a scale
+    factor and rescales every length/area/volume field accordingly, so the
+    output is identical to what a fresh extraction at the new voxel size
+    would have produced (modulo float precision).
+
+    Parameters
+    ----------
+    network : dict
+        A network produced by ``snow2``, ``regions_to_network``, or
+        ``magnet``. Must contain ``param.voxel_size`` and ``param.ndim``.
+    voxel_size : scalar
+        The new voxel size, in the same units as the original.
+
+    Returns
+    -------
+    rescaled : dict
+        A new network dict with all dimensional fields rescaled and
+        ``param.voxel_size`` updated to the new value. The input is not
+        mutated.
+
+    Notes
+    -----
+    Anisotropic voxel sizes (i.e., a tuple stored in ``param.voxel_size``)
+    are not supported; if all axes share the same value, that scalar is
+    used, otherwise an exception is raised.
+
+    """
+    if 'param.voxel_size' not in network or 'param.ndim' not in network:
+        msg = (
+            "Network is missing 'param.voxel_size' and/or 'param.ndim'. "
+            "Networks extracted before this metadata was added cannot be "
+            "rescaled automatically; set both keys manually before calling "
+            "this function."
+        )
+        raise KeyError(msg)
+
+    vs_old = np.asarray(network['param.voxel_size'])
+    if vs_old.ndim > 0:
+        if not np.allclose(vs_old, vs_old.flat[0]):
+            raise NotImplementedError(
+                "Anisotropic voxel sizes are not supported by rescale_network"
+            )
+        vs_old = float(vs_old.flat[0])
+    else:
+        vs_old = float(vs_old)
+    factor = float(voxel_size) / vs_old
+    ndim = int(np.asarray(network['param.ndim']).flat[0])
+
+    rescaled = dict(network)
+    for key, value in network.items():
+        p = _scale_exponent(key, ndim)
+        if p:
+            rescaled[key] = np.asarray(value) * factor**p
+    rescaled['param.voxel_size'] = np.asarray(voxel_size)
+    return rescaled

src/porespy/networks/_getnet_orig.py

@@ -277,5 +277,7 @@ def regions_to_network(
         A = np.array(t_area)*(voxel_size**2)
         net['throat.cross_sectional_area'] = A
         net['throat.equivalent_diameter'] = (4*A/np.pi)**(1/2)
+    net['param.voxel_size'] = np.asarray(voxel_size)
+    net['param.ndim'] = np.asarray(ND)
 
     return net

src/porespy/networks/_getnet_para.py

@@ -259,6 +259,9 @@ def regions_to_network_parallel(
         for single_key in (f"{multi_root}_{i}" for i in range(n)):
             del net[single_key]
 
+    net['param.voxel_size'] = np.asarray(voxel_size)
+    net['param.ndim'] = np.asarray(im.ndim)
+
     return net
 
 

src/porespy/networks/_magnet.py

@@ -715,6 +715,8 @@ def junctions_to_network(sk, juncs, throats, dt, throat_area, voxel_size=1):
     net["pore.inscribed_diameter"] = p_diameter * voxel_size
     net["pore.equivalent_diameter"] = p_diameter_equivalent * voxel_size
     net["pore.index"] = np.arange(0, Np)
+    net["param.voxel_size"] = np.asarray(voxel_size)
+    net["param.ndim"] = np.asarray(sk.ndim)
     return net
 
 

test/unit/test_magnet.py

@@ -109,6 +109,20 @@ def test_junctions(self):
         except Exception:
             pass
 
+    def test_metadata_in_extracted_network(self):
+        magnet = ps.networks.magnet(self.blobs3D, voxel_size=2e-6)
+        assert magnet.network["param.voxel_size"] == 2e-6
+        assert magnet.network["param.ndim"] == 3
+
+    def test_rescale_network_matches_fresh_extraction(self):
+        net1 = ps.networks.magnet(self.blobs3D, voxel_size=1).network
+        net1_scaled = ps.networks.rescale_network(net1, voxel_size=4.0)
+        net2 = ps.networks.magnet(self.blobs3D, voxel_size=4.0).network
+        for key in net2:
+            if key.startswith("param."):
+                continue
+            assert np.allclose(net1_scaled[key], net2[key]), f"mismatch on {key}"
+
     def test_throat_area(self):
         im = self.blobs3D
         magnet = ps.networks.magnet(im, throat_area=True)

test/unit/test_network_extraction.py

@@ -119,6 +119,64 @@ def test_planar_2d_image(self):
         assert np.allclose(net1['pore.coords'][:, 1], net2['pore.coords'][:, 2])
         assert np.allclose(net1['pore.coords'][:, 0], net3['pore.coords'][:, 1])
 
+    def test_metadata_in_extracted_network(self):
+        im = self.snow3d.regions*self.im3d
+        net = ps.networks.regions_to_network(im, voxel_size=1.5e-6)
+        assert net['param.voxel_size'] == 1.5e-6
+        assert net['param.ndim'] == 3
+
+        net2d = ps.networks.regions_to_network(self.snow.regions*self.im)
+        assert net2d['param.ndim'] == 2
+
+        snow = ps.networks.snow2(self.im3d, voxel_size=2e-6)
+        assert snow.network['param.voxel_size'] == 2e-6
+        assert snow.network['param.ndim'] == 3
+
+    def test_rescale_network_matches_fresh_extraction(self):
+        im = self.snow3d.regions*self.im3d
+        vs = 3.0
+        net1 = ps.networks.regions_to_network(im, voxel_size=1)
+        net1_scaled = ps.networks.rescale_network(net1, voxel_size=vs)
+        net2 = ps.networks.regions_to_network(im, voxel_size=vs)
+        for key in net2:
+            if key.startswith('param.'):
+                continue
+            assert np.allclose(net1_scaled[key], net2[key]), f"mismatch on {key}"
+
+    def test_rescale_network_2d(self):
+        im = self.snow.regions*self.im
+        net1 = ps.networks.regions_to_network(im, voxel_size=1)
+        net2 = ps.networks.regions_to_network(im, voxel_size=2.5)
+        net1_scaled = ps.networks.rescale_network(net1, voxel_size=2.5)
+        for key in net2:
+            if key.startswith('param.'):
+                continue
+            assert np.allclose(net1_scaled[key], net2[key]), f"mismatch on {key}"
+
+    def test_rescale_network_does_not_mutate(self):
+        im = self.snow3d.regions*self.im3d
+        net = ps.networks.regions_to_network(im, voxel_size=1)
+        coords_before = net['pore.coords'].copy()
+        _ = ps.networks.rescale_network(net, voxel_size=5)
+        assert np.array_equal(net['pore.coords'], coords_before)
+        assert net['param.voxel_size'] == 1
+
+    def test_rescale_network_is_idempotent(self):
+        im = self.snow3d.regions*self.im3d
+        net = ps.networks.regions_to_network(im, voxel_size=1)
+        once = ps.networks.rescale_network(net, voxel_size=2)
+        twice = ps.networks.rescale_network(once, voxel_size=7)
+        direct = ps.networks.rescale_network(net, voxel_size=7)
+        for key in direct:
+            assert np.allclose(twice[key], direct[key]), f"mismatch on {key}"
+
+    def test_rescale_network_missing_metadata_raises(self):
+        im = self.snow3d.regions*self.im3d
+        net = ps.networks.regions_to_network(im, voxel_size=1)
+        del net['param.voxel_size']
+        with pytest.raises(KeyError):
+            ps.networks.rescale_network(net, voxel_size=2)
+
 
 if __name__ == '__main__':
     t = NetworkExtractionTest()