diff --git a/src/porespy/networks/_funcs.py b/src/porespy/networks/_funcs.py index 28b3c1599..34132fffd 100644 --- a/src/porespy/networks/_funcs.py +++ b/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 diff --git a/src/porespy/networks/_getnet_orig.py b/src/porespy/networks/_getnet_orig.py index 1f6abe149..755d6f27c 100644 --- a/src/porespy/networks/_getnet_orig.py +++ b/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 diff --git a/src/porespy/networks/_getnet_para.py b/src/porespy/networks/_getnet_para.py index bd78ec938..c4436dce3 100644 --- a/src/porespy/networks/_getnet_para.py +++ b/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 diff --git a/src/porespy/networks/_magnet.py b/src/porespy/networks/_magnet.py index af746b1a4..60efe83aa 100644 --- a/src/porespy/networks/_magnet.py +++ b/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 diff --git a/test/unit/test_magnet.py b/test/unit/test_magnet.py index c34f27660..331ea4f91 100644 --- a/test/unit/test_magnet.py +++ b/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) diff --git a/test/unit/test_network_extraction.py b/test/unit/test_network_extraction.py index ab4f439d2..35388af81 100644 --- a/test/unit/test_network_extraction.py +++ b/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()