From f6734d98c45278427d373fd3be573d116655ba23 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Mon, 6 Mar 2023 13:25:33 -0500 Subject: [PATCH 01/16] adding rectangular_pillars function --- porespy/generators/__init__.py | 1 + porespy/generators/_micromodels.py | 125 +++++++++++++++++++++++++++++ 2 files changed, 126 insertions(+) create mode 100644 porespy/generators/_micromodels.py diff --git a/porespy/generators/__init__.py b/porespy/generators/__init__.py index 9957a27e8..98ecce335 100644 --- a/porespy/generators/__init__.py +++ b/porespy/generators/__init__.py @@ -51,3 +51,4 @@ from ._noise import fractal_noise from ._borders import * from ._fractals import * +from ._micromodels import * diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py new file mode 100644 index 000000000..84814f3e7 --- /dev/null +++ b/porespy/generators/_micromodels.py @@ -0,0 +1,125 @@ +# import porespy as ps +import numpy as np +import scipy.ndimage as spim +import scipy.spatial as sptl +from porespy.tools import ps_rect, ps_round, extend_slice +from porespy.generators import lattice_spheres, line_segment + + +__all__ = [ + 'rectangular_pillars', +] + + +def cross(r, t=0): + cr = np.zeros([2*r+1, 2*r+1], dtype=bool) + cr[r-t:r+t+1, :] = True + cr[:, r-t:r+t+1] = True + return cr + + +def ex(r, t=0): + x = np.eye(2*r + 1).astype(bool) + x += np.fliplr(x) + x = spim.binary_dilation(x, structure=ps_rect(w=2*t+1, ndim=2)) + return x + + +def rectangular_pillars(shape=[5, 5], spacing=30, Rmin=2, Rmax=20, lattice='sc'): + r""" + A 2D micromodel with rectangular pillars arranged on a regular lattice + + Parameters + ---------- + shape : list + The number of pillars in the x and y directions. The size of the of the + image will be dictated by the ``spacing`` argument. + spacing : int + The number of pixels between neighboring pores centers. + Rmin : int + The minimum size of the openings between pillars in pixels + Rmax : int + The maximum size of the openings between pillars in pixels + lattice : str + The type of lattice to use. Options are: + + ======== =================================================================== + lattice description + ======== =================================================================== + 'sc' A simple cubic lattice where the pillars are aligned vertically and + horizontally with the standard grid. In this case the meaning of + ``spacing``, ``Rmin`` and ``Rmax`` directly refers to the number of + pixels. + 'tri' A triangular matrix, which is esentially a cubic matrix rotated 45 + degrees. In this case the mean of ``spacing``, ``Rmin`` and ``Rmax`` + refer to the length of a pixel. + ======== =================================================================== + + Returns + ------- + ims : dataclass + Several images are generated in internally, so they are all returned as + attributes of a dataclass-like object. The attributes are as follows: + + ========== ================================================================= + attribute description + ========== ================================================================= + im A 2D image whose size is dictated by the number of pillars + (given by ``shape``) and the ``spacing`` between them. + centers An image the same size as ``im`` with ``True`` values marking + the center of each pore body. + edges An image the same size as ``im`` with ``True`` values marking + the edges connecting the pore centers. Note that the ``centers`` + have been removed from this image. + ========== ================================================================= + + Examples + -------- + `Click here + `_ + to view online example. + """ + if lattice == 'sc': + strel = cross + Rmax = Rmax + 1 + else: + strel = ex + shape = np.array(shape) - 1 + Rmin = int(Rmin*np.sin(np.deg2rad(45))) + Rmax = int((Rmax-2)*np.sin(np.deg2rad(45))) + centers = ~lattice_spheres( + shape=[shape[0]*spacing+1, shape[1]*spacing+1], + spacing=spacing, + r=1, + offset=0, + lattice=lattice) + Rmin = max(1, Rmin) + crds = np.where(centers) + tri = sptl.Delaunay(np.vstack(crds).T) + edges = np.zeros_like(centers, dtype=bool) + for s in tri.simplices: + s2 = s.tolist() + s2.append(s[0]) + for i in range(len(s)): + P1, P2 = tri.points[s2[i]], tri.points[s2[i+1]] + L = np.sqrt(np.sum(np.square(np.subtract(P1, P2)))) + if ((lattice == 'tri') and (L < spacing)) \ + or ((lattice == 'sc') and (L <= spacing)): + crds = line_segment(P1, P2) + edges[tuple(crds)] = True + temp = spim.binary_dilation(centers, structure=ps_rect(w=1, ndim=2)) + edges = edges*~temp + if lattice == 'sc': + labels, N = spim.label(edges, structure=ps_round(r=1, ndim=2, smooth=False)) + else: + labels, N = spim.label(edges, structure=ps_rect(w=3, ndim=2)) + slices = spim.find_objects(labels) + throats = np.zeros_like(edges, dtype=int) + for i, s in enumerate(slices): + r = np.random.randint(Rmin, Rmax) + s2 = extend_slice(s, throats.shape, pad=2*r+1) + mask = labels[s2] == (i + 1) + t = spim.binary_dilation(mask, structure=strel(r=r, t=1)) + throats[s2] += t + micromodel = throats > 0 + return micromodel, edges, centers From 8059d0c24ac86bc2915b1a0fbfa6f76f755bfedf Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Mon, 6 Mar 2023 20:49:14 -0500 Subject: [PATCH 02/16] adding magnet2 from stash --- porespy/networks/__init__.py | 1 + porespy/networks/_magnet.py | 300 +++++++++++++++++++++++++++++++++++ 2 files changed, 301 insertions(+) create mode 100644 porespy/networks/_magnet.py diff --git a/porespy/networks/__init__.py b/porespy/networks/__init__.py index 50b336fb5..2df1a1d1f 100644 --- a/porespy/networks/__init__.py +++ b/porespy/networks/__init__.py @@ -36,3 +36,4 @@ from ._size_factors import diffusive_size_factor_AI from ._size_factors import create_model from ._size_factors import find_conns +from ._magnet import * diff --git a/porespy/networks/_magnet.py b/porespy/networks/_magnet.py new file mode 100644 index 000000000..ab1e1f1c9 --- /dev/null +++ b/porespy/networks/_magnet.py @@ -0,0 +1,300 @@ +import numpy as np +from skimage.segmentation import find_boundaries +from skimage.morphology import skeletonize_3d +import porespy as ps +from edt import edt +import scipy.ndimage as spim +from porespy.tools import get_tqdm, Results +import pandas as pd +# import openpnm as op +# import matplotlib.pyplot as plt + + +__all__ = [ + 'magnet2', +] + + +tqdm = get_tqdm() + + +def analyze_skeleton_2(sk, dt): + # Blur the DT + dt2 = spim.gaussian_filter(dt, sigma=0.4) + # Dilate the skeleton (probably not needed) + # strel = ps.tools.ps_round(r=1, ndim=im.ndim, smooth=False) + # sk2 = spim.binary_dilation(sk, structure=strel) + # Run maximum filter on dt + strel = ps.tools.ps_round(r=3, ndim=sk.ndim, smooth=False) + dt3 = spim.maximum_filter(dt2, footprint=strel) + # Multiply skeleton by smoothed and filtered dt + sk3 = sk*dt3 + # Find peaks on sk3 + strel = ps.tools.ps_round(r=5, ndim=sk.ndim, smooth=False) + peaks = (spim.maximum_filter(sk3, footprint=strel) == dt3)*sk + return peaks + + +def magnet2(im): + im = ps.filters.fill_blind_pores(im, surface=True) + if im.ndim == 3: + im = ps.filters.trim_floating_solid(im, conn=2*im.ndim, surface=True) + sk = skeletonize_3d(im) > 0 + sk_orig = np.copy(sk) + dt = edt(im) + dt = spim.maximum_filter(dt, size=3) + spheres = np.zeros_like(im, dtype=int) + centers = np.zeros_like(im, dtype=int) + jcts = ps.filters.analyze_skeleton(sk) + peaks = analyze_skeleton_2(sk, dt) + + # %% Insert spheres and center points into image, and delete underlying skeleton + crds = np.vstack(np.where(jcts.endpts + jcts.juncs + peaks)).T + inds = np.argsort(dt[tuple(crds.T)])[-1::-1] + crds = crds[inds, :] + count = 0 + for i, row in enumerate(tqdm(crds)): + r = int(dt[tuple(row)]) + if spheres[tuple(row)] == 0: + count += 1 + ps.tools._insert_disk_at_points( + im=sk, + coords=np.atleast_2d(row).T, + r=r, + v=False, + smooth=False, + overwrite=True) + ps.tools._insert_disk_at_points( + im=centers, + coords=np.atleast_2d(row).T, + r=1, + v=1, + smooth=True, + overwrite=False) + ps.tools._insert_disk_at_points( + im=spheres, + coords=np.atleast_2d(row).T, + r=r, + v=count, + smooth=False, + overwrite=False) + + # %% Add skeleton to edges/intersections of overlapping spheres + temp = find_boundaries(spheres, mode='thick') + sk += temp*sk_orig + + # %% Analyze image to extract pore and throat info + pore_labels = np.copy(spheres) + centers = centers*pore_labels + strel = ps.tools.ps_rect(w=3, ndim=centers.ndim) + throat_labels, Nt = spim.label(input=sk > 0, structure=strel) + pore_slices = spim.find_objects(pore_labels) + throat_slices = spim.find_objects(throat_labels) + + # %% Get pore coordinates and diameters + coords = [] + pore_diameters = [] + for i, p in enumerate(pore_slices): + inds = np.vstack(np.where(centers[p] == (i + 1))).T[0, :] + pore_diameters.append(2*dt[p][tuple(inds)]) + inds = inds + np.array([s.start for s in p]) + coords.append(inds.tolist()) + pore_diameters = np.array(pore_diameters, dtype=float) + coords = np.vstack(coords).astype(float) + + # %% Get throat connections and diameters + conns = [] + throat_diameters = [] + for i, t in enumerate(throat_slices): + s = ps.tools.extend_slice(t, shape=im.shape, pad=1) + mask = throat_labels[s] == (i + 1) + mask_dil = spim.binary_dilation(mask, structure=strel)*sk_orig[s] + neighbors = np.unique(pore_labels[s]*mask_dil)[1:] + Dt = 2*dt[s][mask].min() + if len(neighbors) == 2: + conns.append(neighbors.tolist()) + throat_diameters.append(Dt) + elif len(neighbors) > 2: + inds = np.argsort(pore_diameters[neighbors-1])[-1::-1] + inds = neighbors[inds] + temp = [[inds[0], inds[j+1]] for j in range(len(inds)-1)] + conns.extend(temp) + # The following is a temporary shortcut and needs to be done properly + temp = [Dt for _ in range(len(inds)-1)] + throat_diameters.extend(temp) + else: + pass + throat_diameters = np.array(throat_diameters, dtype=float) + # Move to upper triangular and increment to 0 indexing + conns = np.sort(np.vstack(conns), axis=1) - 1 + # Remove duplicate throats + hits = pd.DataFrame(conns).duplicated().to_numpy() + conns = conns[~hits, :] + throat_diameters = throat_diameters[~hits] + + # %% Store in openpnm compatible dictionary + net = {} + if coords.shape[1] == 2: + coords = np.vstack((coords[:, 0], coords[:, 1], np.zeros_like(coords[:, 0]))).T + net['pore.coords'] = coords + net['throat.conns'] = conns + net['pore.diameter'] = pore_diameters + net['throat.diameter'] = throat_diameters + net['pore.all'] = np.ones([coords.shape[0], ], dtype=bool) + net['throat.all'] = np.ones([conns.shape[0], ], dtype=bool) + net['pore.xmin'] = coords[:, 0] < 0.1*(coords[:, 0].max() - coords[:, 0].min()) + net['pore.xmax'] = coords[:, 0] > 0.9*(coords[:, 0].max() - coords[:, 0].min()) + + results = Results() + results.network = net + results.centers = centers + results.spheres = spheres + results.skeleton = sk_orig + results.im = im + return results + + + + + +# %% +if __name__ == "__main__": + import openpnm as op + import matplotlib.pyplot as plt + np.random.seed(0) + im = ps.generators.blobs([200, 200, 200], blobiness=0.5, porosity=0.7) + im = ps.filters.fill_blind_pores(im, conn=2*im.ndim, surface=True) + im = ps.filters.trim_floating_solid(im, conn=2*im.ndim, surface=True) + net = magnet2(im) + net2 = ps.networks.snow2(im, boundary_width=0) + + # %% + pn_m = op.io.network_from_porespy(net.network) + pn_s = op.io.network_from_porespy(net2.network) + print(pn_m) + print(pn_s) + pn_s['pore.diameter'] = pn_s['pore.inscribed_diameter'] + pn_s['throat.diameter'] = pn_s['throat.inscribed_diameter'] + coords = pn_s.coords + pn_s['pore.xmin'] = coords[:, 0] < 0.1*(coords[:, 0].max() - coords[:, 0].min()) + pn_s['pore.xmax'] = coords[:, 0] > 0.9*(coords[:, 0].max() - coords[:, 0].min()) + h = op.utils.check_network_health(pn_s) + op.topotools.trim(network=pn_s, pores=h['disconnected_pores']) + h = op.utils.check_network_health(pn_m) + op.topotools.trim(network=pn_m, pores=h['disconnected_pores']) + pn_s.regenerate_models() + pn_m.regenerate_models() + pn_s.add_model_collection(op.models.collections.geometry.snow) + pn_s.regenerate_models() + pn_m.add_model_collection(op.models.collections.geometry.magnet) + pn_m.regenerate_models() + + # %% + if 0: + for i in range(100): + Dt = pn_m['throat.diameter'] == pn_m['throat.diameter'].max() + Lt = pn_m['throat.length'] == 1e-15 + T = np.where(Dt*Lt)[0][0] + P1, P2 = pn_m.conns[T] + op.topotools.merge_pores(network=pn_m, pores=[P1, P2]) + + # %% + fig, ax = plt.subplots(2, 2) + kw = {'edgecolor': 'k', 'bins': 20, 'alpha': 0.5, 'density': True, 'cumulative': True} + ax[0][0].hist(pn_s['pore.diameter'], color='b', label='snow', **kw) + ax[0][0].hist(pn_m['pore.diameter'], color='r', label='magnet', **kw) + ax[0][0].set_xlabel('Pore Diameter') + ax[0][0].legend() + ax[0][1].hist(pn_s['throat.diameter'], color='b', label='snow', **kw) + ax[0][1].hist(pn_m['throat.diameter'], color='r', label='magnet', **kw) + ax[0][1].set_xlabel('Throat Diameter') + ax[0][1].legend() + ax[1][0].hist(pn_s['throat.length'], color='b', label='snow', **kw) + ax[1][0].hist(pn_m['throat.length'], color='r', label='magnet', **kw) + ax[1][0].set_xlabel('Throat Length') + ax[1][0].legend() + ax[1][1].hist(pn_s['pore.coordination_number'], color='b', label='snow', **kw) + ax[1][1].hist(pn_m['pore.coordination_number'], color='r', label='magnet', **kw) + ax[1][1].set_xlabel('Coordination Number') + ax[1][1].legend() + + # %% + w_s = op.phase.Water(network=pn_s) + w_s['pore.diffusivity'] = 1.0 + w_s.add_model_collection(op.models.collections.physics.standard) + w_s.regenerate_models() + w_m = op.phase.Water(network=pn_m) + w_m['pore.diffusivity'] = 1.0 + w_m.add_model_collection(op.models.collections.physics.standard) + w_m.regenerate_models() + + # %% + fig, ax = plt.subplots(2, 2) + kw = {'edgecolor': 'k', 'bins': 20, 'alpha': 0.5, 'density': True, 'cumulative': True} + ax[0][0].hist(w_s['throat.entry_pressure'], color='b', label='snow', **kw) + ax[0][0].hist(w_m['throat.entry_pressure'], color='r', label='magnet', **kw) + ax[0][0].set_xlabel('Throat Entry Pressure') + ax[0][0].legend() + ax[0][1].hist(w_s['throat.hydraulic_conductance'], color='b', label='snow', **kw) + ax[0][1].hist(w_m['throat.hydraulic_conductance'], color='r', label='magnet', **kw) + ax[0][1].set_xlabel('Throat Hydraulic Conductance') + ax[0][1].legend() + ax[1][0].plot(pn_s['throat.diameter'], pn_s['pore.diameter'][pn_s.conns][:, 1], 'b.', label='snow') + ax[1][0].plot(pn_m['throat.diameter'], pn_m['pore.diameter'][pn_m.conns][:, 1], 'r.', label='magnet') + ax[1][0].plot([0, 20], [0, 20], 'k-') + ax[1][0].set_xlabel('Throat Diameter') + ax[1][0].set_ylabel('Pore Diameter') + ax[1][0].legend() + + # %% + sf_s = op.algorithms.StokesFlow(network=pn_s, phase=w_s) + sf_s.set_value_BC(pores=pn_s.pores('xmin'), values=1.0) + sf_s.set_value_BC(pores=pn_s.pores('xmax'), values=0.0) + sf_s.run() + print(sf_s.rate(pores=pn_s.pores('xmin'), mode='group')) + + sf_m = op.algorithms.StokesFlow(network=pn_m, phase=w_m) + sf_m.set_value_BC(pores=pn_m.pores('xmin'), values=1.0) + sf_m.set_value_BC(pores=pn_m.pores('xmax'), values=0.0) + sf_m.run() + print(sf_m.rate(pores=pn_m.pores('xmin'), mode='group')) + + # %% + pc_s = op.algorithms.Drainage(network=pn_s, phase=w_s) + pc_s.set_inlet_BC(pores=pn_s.pores('xmin')) + pc_s.run() + + pc_m = op.algorithms.Drainage(network=pn_m, phase=w_m) + pc_m.set_inlet_BC(pores=pn_m.pores('xmin')) + pc_m.run() + + ax[1][1].plot(pc_s.pc_curve().pc,pc_s.pc_curve().snwp, 'b-o', label='snow') + ax[1][1].plot(pc_m.pc_curve().pc,pc_m.pc_curve().snwp, 'r-o', label='magnet') + ax[1][1].legend() + ax[1][1].set_xlabel('Capillary Pressure') + ax[1][1].set_ylabel('Non-Wetting Phase Saturation') + ax[1][1].legend() + + # %% + fd_s = op.algorithms.FickianDiffusion(network=pn_s, phase=w_s) + fd_s.set_value_BC(pores=pn_s.pores('xmin'), values=1.0) + fd_s.set_value_BC(pores=pn_s.pores('xmax'), values=0.0) + fd_s.run() + Deff = fd_s.rate(pores=pn_s.pores('xmin'))*im.shape[0]/(im.shape[1]*im.shape[2]) + taux_s = (im.sum()/im.size)/Deff + print(taux_s) + + fd_m = op.algorithms.FickianDiffusion(network=pn_m, phase=w_m) + fd_m.set_value_BC(pores=pn_m.pores('xmin'), values=1.0) + fd_m.set_value_BC(pores=pn_m.pores('xmax'), values=0.0) + fd_m.run() + Deff = fd_m.rate(pores=pn_m.pores('xmin'))*im.shape[0]/(im.shape[1]*im.shape[2]) + taux_m = (im.sum()/im.size)/Deff + print(taux_m) + + + + + + + From 5f5a67bdde25801f57247c25cee9a89fc1de6b26 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Mon, 6 Mar 2023 21:00:45 -0500 Subject: [PATCH 03/16] add ability to pass in custom skel --- porespy/networks/_magnet.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/porespy/networks/_magnet.py b/porespy/networks/_magnet.py index ab1e1f1c9..b68935f22 100644 --- a/porespy/networks/_magnet.py +++ b/porespy/networks/_magnet.py @@ -35,11 +35,12 @@ def analyze_skeleton_2(sk, dt): return peaks -def magnet2(im): - im = ps.filters.fill_blind_pores(im, surface=True) - if im.ndim == 3: - im = ps.filters.trim_floating_solid(im, conn=2*im.ndim, surface=True) - sk = skeletonize_3d(im) > 0 +def magnet2(im, sk=None): + if sk is None: + im = ps.filters.fill_blind_pores(im, surface=True) + if im.ndim == 3: + im = ps.filters.trim_floating_solid(im, conn=2*im.ndim, surface=True) + sk = skeletonize_3d(im) > 0 sk_orig = np.copy(sk) dt = edt(im) dt = spim.maximum_filter(dt, size=3) From 9640826945bc68913b380e7b5d39c0a4ab56c568 Mon Sep 17 00:00:00 2001 From: jgostick Date: Tue, 7 Mar 2023 20:24:05 -0500 Subject: [PATCH 04/16] adding junction code to analyze_skeleton_2 func --- porespy/networks/_magnet.py | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/porespy/networks/_magnet.py b/porespy/networks/_magnet.py index b68935f22..222806e40 100644 --- a/porespy/networks/_magnet.py +++ b/porespy/networks/_magnet.py @@ -1,9 +1,11 @@ import numpy as np +import scipy as sp from skimage.segmentation import find_boundaries -from skimage.morphology import skeletonize_3d +from skimage.morphology import skeletonize_3d, square, cube import porespy as ps from edt import edt import scipy.ndimage as spim +from porespy.filters import reduce_peaks from porespy.tools import get_tqdm, Results import pandas as pd # import openpnm as op @@ -19,11 +21,28 @@ def analyze_skeleton_2(sk, dt): + # kernel for convolution + if sk.ndim == 2: + a = square(3) + else: + a = cube(3) + # compute convolution directly or via fft, whichever is fastest + conv = sp.signal.convolve(sk*1.0, a, mode='same', method='auto') + conv = np.rint(conv).astype(int) # in case of fft, accuracy is lost + # find junction points of skeleton + juncs = (conv >= 4) * sk + # find endpoints of skeleton + end_pts = (conv == 2) * sk + # reduce cluster of junctions to single pixel at centre + juncs_r = reduce_peaks(juncs) + # results object + pt = Results() + pt.juncs = juncs + pt.endpts = end_pts + pt.juncs_r = juncs_r + # Blur the DT dt2 = spim.gaussian_filter(dt, sigma=0.4) - # Dilate the skeleton (probably not needed) - # strel = ps.tools.ps_round(r=1, ndim=im.ndim, smooth=False) - # sk2 = spim.binary_dilation(sk, structure=strel) # Run maximum filter on dt strel = ps.tools.ps_round(r=3, ndim=sk.ndim, smooth=False) dt3 = spim.maximum_filter(dt2, footprint=strel) @@ -32,7 +51,8 @@ def analyze_skeleton_2(sk, dt): # Find peaks on sk3 strel = ps.tools.ps_round(r=5, ndim=sk.ndim, smooth=False) peaks = (spim.maximum_filter(sk3, footprint=strel) == dt3)*sk - return peaks + pt.peaks = peaks + return pt def magnet2(im, sk=None): @@ -46,8 +66,8 @@ def magnet2(im, sk=None): dt = spim.maximum_filter(dt, size=3) spheres = np.zeros_like(im, dtype=int) centers = np.zeros_like(im, dtype=int) - jcts = ps.filters.analyze_skeleton(sk) - peaks = analyze_skeleton_2(sk, dt) + jcts = analyze_skeleton_2(sk, dt) + peaks = jcts.peaks # %% Insert spheres and center points into image, and delete underlying skeleton crds = np.vstack(np.where(jcts.endpts + jcts.juncs + peaks)).T @@ -132,6 +152,7 @@ def magnet2(im, sk=None): hits = pd.DataFrame(conns).duplicated().to_numpy() conns = conns[~hits, :] throat_diameters = throat_diameters[~hits] + sk = sk_orig # %% Store in openpnm compatible dictionary net = {} From 632e15304ea3dddf124a81e4c2738beef983d634 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 9 Mar 2023 10:08:11 -0500 Subject: [PATCH 05/16] adding tqdm to for loops since one is annoying slow for big images --- porespy/generators/_micromodels.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py index 84814f3e7..919196c2a 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/generators/_micromodels.py @@ -2,8 +2,9 @@ import numpy as np import scipy.ndimage as spim import scipy.spatial as sptl -from porespy.tools import ps_rect, ps_round, extend_slice +from porespy.tools import ps_rect, ps_round, extend_slice, get_tqdm from porespy.generators import lattice_spheres, line_segment +from porespy import settings __all__ = [ @@ -11,6 +12,9 @@ ] +tqdm = get_tqdm() + + def cross(r, t=0): cr = np.zeros([2*r+1, 2*r+1], dtype=bool) cr[r-t:r+t+1, :] = True @@ -97,7 +101,8 @@ def rectangular_pillars(shape=[5, 5], spacing=30, Rmin=2, Rmax=20, lattice='sc') crds = np.where(centers) tri = sptl.Delaunay(np.vstack(crds).T) edges = np.zeros_like(centers, dtype=bool) - for s in tri.simplices: + msg = 'Adding edges of triangulation to image' + for s in tqdm(tri.simplices, msg, **settings.tqdm): s2 = s.tolist() s2.append(s[0]) for i in range(len(s)): @@ -115,7 +120,8 @@ def rectangular_pillars(shape=[5, 5], spacing=30, Rmin=2, Rmax=20, lattice='sc') labels, N = spim.label(edges, structure=ps_rect(w=3, ndim=2)) slices = spim.find_objects(labels) throats = np.zeros_like(edges, dtype=int) - for i, s in enumerate(slices): + msg = 'Dilating edges to random widths' + for i, s in enumerate(tqdm(slices, msg, **settings.tqdm)): r = np.random.randint(Rmin, Rmax) s2 = extend_slice(s, throats.shape, pad=2*r+1) mask = labels[s2] == (i + 1) From fbf366e138b9a9a91d21b2fb4442321822bf8b59 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 9 Mar 2023 10:20:34 -0500 Subject: [PATCH 06/16] making lattice argument more forgiving --- porespy/generators/_micromodels.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py index 919196c2a..441999cdc 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/generators/_micromodels.py @@ -83,14 +83,18 @@ def rectangular_pillars(shape=[5, 5], spacing=30, Rmin=2, Rmax=20, lattice='sc') `_ to view online example. """ - if lattice == 'sc': + if lattice.startswith('s'): strel = cross Rmax = Rmax + 1 - else: + lattice = 'sc' # In case user specified s, sq or square, etc. + elif lattice.startswith('t'): strel = ex shape = np.array(shape) - 1 Rmin = int(Rmin*np.sin(np.deg2rad(45))) Rmax = int((Rmax-2)*np.sin(np.deg2rad(45))) + lattice = 'tri' # In case user specified t, or triangle, etc. + else: + raise Exception(f"Unrecognized lattice type {lattice}") centers = ~lattice_spheres( shape=[shape[0]*spacing+1, shape[1]*spacing+1], spacing=spacing, From 2f86611352abac12e1acc922ee4c3bf871d95128 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 9 Mar 2023 10:37:13 -0500 Subject: [PATCH 07/16] adding extra args for return_edges and return_centers --- porespy/generators/_micromodels.py | 38 +++++++++++++++++++++++++----- 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py index 441999cdc..95264d360 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/generators/_micromodels.py @@ -2,7 +2,7 @@ import numpy as np import scipy.ndimage as spim import scipy.spatial as sptl -from porespy.tools import ps_rect, ps_round, extend_slice, get_tqdm +from porespy.tools import ps_rect, ps_round, extend_slice, get_tqdm, Results from porespy.generators import lattice_spheres, line_segment from porespy import settings @@ -29,7 +29,15 @@ def ex(r, t=0): return x -def rectangular_pillars(shape=[5, 5], spacing=30, Rmin=2, Rmax=20, lattice='sc'): +def rectangular_pillars( + shape=[5, 5], + spacing=30, + Rmin=5, + Rmax=15, + lattice='sc', + return_edges=False, + return_centers=False +): r""" A 2D micromodel with rectangular pillars arranged on a regular lattice @@ -59,11 +67,20 @@ def rectangular_pillars(shape=[5, 5], spacing=30, Rmin=2, Rmax=20, lattice='sc') refer to the length of a pixel. ======== =================================================================== + return_edges : boolean, optional, default is ``False`` + If ``True`` then an image of of the edges between each pore center is also + returned along with the micromodel + return_centers : boolean, optional, default is ``False`` + If ``True`` then an image with marks located at each pore center is also + return along with the micromodel + Returns ------- - ims : dataclass - Several images are generated in internally, so they are all returned as - attributes of a dataclass-like object. The attributes are as follows: + im or ims : ndarray or dataclass + If ``return_centers`` and ``return_edges`` are both ``False``, then only + an ndarray of the micromodel is returned. If either or both are ``True`` + then a ``dataclass-like`` object is return with multiple images attached + as attributes: ========== ================================================================= attribute description @@ -132,4 +149,13 @@ def rectangular_pillars(shape=[5, 5], spacing=30, Rmin=2, Rmax=20, lattice='sc') t = spim.binary_dilation(mask, structure=strel(r=r, t=1)) throats[s2] += t micromodel = throats > 0 - return micromodel, edges, centers + if (not return_edges) and (not return_centers): + return micromodel + else: + ims = Results() + ims.im = micromodel + if return_edges: + ims.edges = edges + if return_centers: + ims.centers = centers + return ims From 13c36efefa865b300303880b0fe1fc0916cce8eb Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 9 Mar 2023 10:37:22 -0500 Subject: [PATCH 08/16] adding examle notebook [wip] --- .../reference/rectangular_pillars.ipynb | 116 ++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 examples/generators/reference/rectangular_pillars.ipynb diff --git a/examples/generators/reference/rectangular_pillars.ipynb b/examples/generators/reference/rectangular_pillars.ipynb new file mode 100644 index 000000000..2f2884462 --- /dev/null +++ b/examples/generators/reference/rectangular_pillars.ipynb @@ -0,0 +1,116 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "96e8a38c", + "metadata": {}, + "outputs": [], + "source": [ + "import porespy as ps\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "id": "dfbea3b7", + "metadata": {}, + "source": [ + "## Default Arguments\n", + "\n", + "The function returns a sample image without supplying any arguments. This is a useful way to begin experimenting with the function. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c57eeb7f", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Adding edges of triangulation to image: 0%| | 0/50 [00:00" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=[5, 5])\n", + "ax.imshow(im, interpolation='none')\n", + "ax.axis(False);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b481f8b5", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.15" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 956717102d055f50cf0f9200f991b69ddf268e7b Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Sat, 11 Mar 2023 15:48:02 -0500 Subject: [PATCH 09/16] adding preliminary random pillar generator --- porespy/generators/_micromodels.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py index 95264d360..ca3a52120 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/generators/_micromodels.py @@ -159,3 +159,20 @@ def rectangular_pillars( if return_centers: ims.centers = centers return ims + + +def random_cylindrical_pillars(shape=[100, 100]): + from nanomesh import Mesher2D + from porespy.generators import borders + + im = np.ones([50, 50], dtype=float) + bd = borders(im.shape, mode='faces') + im[bd] = 0.0 + + mesher = Mesher2D(im) + mesher.generate_contour(max_edge_dist=50, level=0.999) + + mesh = mesher.triangulate(opts='q20a5ne') + # mesh.plot_pyvista(jupyter_backend='static', show_edges=True) + tri = mesh.triangle_dict + return tri From 14635baa417f28be17bb372522de2c08e48ccc2d Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Sat, 11 Mar 2023 22:16:34 -0500 Subject: [PATCH 10/16] put tri into image [wip] --- porespy/generators/_micromodels.py | 40 +++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py index ca3a52120..0d808bf75 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/generators/_micromodels.py @@ -161,18 +161,52 @@ def rectangular_pillars( return ims -def random_cylindrical_pillars(shape=[100, 100]): +def random_cylindrical_pillars(shape=[100, 100], f=0.45): from nanomesh import Mesher2D from porespy.generators import borders - im = np.ones([50, 50], dtype=float) + im = np.ones(shape, dtype=float) bd = borders(im.shape, mode='faces') im[bd] = 0.0 mesher = Mesher2D(im) mesher.generate_contour(max_edge_dist=50, level=0.999) - mesh = mesher.triangulate(opts='q20a5ne') + mesh = mesher.triangulate(opts='q1a50ne') # mesh.plot_pyvista(jupyter_backend='static', show_edges=True) tri = mesh.triangle_dict + + r_max = np.inf*np.ones([tri['vertices'].shape[0], ]) + for e in tri['edges']: + L = np.sqrt(np.sum(np.diff(tri['vertices'][e], axis=0)**2)) + if tri['vertex_markers'][e[0]] == 0: + r_max[e[0]] = min(r_max[e[0]], L/2) + if tri['vertex_markers'][e[1]] == 0: + r_max[e[1]] = min(r_max[e[1]], L/2) + + mask = np.ravel(tri['vertex_markers'] == 0) + r = f*(2*r_max[mask]) + + coords = np.vstack((tri['vertices'][mask].T, r)).T + + + + + + + + + + + + + + + + + + + + + return tri From 425816a682b43b5dfb2b48c5f2c4135964c75a85 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Mon, 13 Mar 2023 01:17:07 -0400 Subject: [PATCH 11/16] working on micromodel generators --- porespy/generators/_micromodels.py | 75 +++++++++++++++++++----------- 1 file changed, 48 insertions(+), 27 deletions(-) diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py index bc66a2274..7a0bc2714 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/generators/_micromodels.py @@ -3,6 +3,7 @@ import scipy.ndimage as spim import scipy.spatial as sptl from porespy.tools import ps_rect, ps_round, extend_slice, get_tqdm, Results +from porespy.tools import _insert_disks_at_points from porespy.generators import lattice_spheres, line_segment from porespy import settings @@ -161,13 +162,42 @@ def rectangular_pillars( return ims +def points_to_spheres(im): + from scipy.spatial import distance_matrix + if im.ndim == 3: + x, y, z = np.where(im > 0) + coords = np.vstack((x, y, z)).T + else: + x, y = np.where(im > 0) + coords = np.vstack((x, y)) + if im.dtype == bool: + dmap = distance_matrix(coords.T, coords.T) + mask = dmap < 1 + dmap[mask] = np.inf + r = np.around(dmap.min(axis=0)/2, decimals=0).astype(int) + else: + r = im[x, y].flatten() + im_spheres = np.zeros_like(im, dtype=bool) + im_spheres = _insert_disks_at_points( + im_spheres, + coords=coords, + radii=r, + v=True, + smooth=False, + ) + return im_spheres + + def random_cylindrical_pillars( - shape=[100, 100], + shape=[1500, 1500], f=0.45, + a=1500, ): from nanomesh import Mesher2D from porespy.generators import borders, spheres_from_coords + if len(shape) != 2: + raise Exception("Shape must be 2D") im = np.ones(shape, dtype=float) bd = borders(im.shape, mode='faces') im[bd] = 0.0 @@ -175,7 +205,7 @@ def random_cylindrical_pillars( mesher = Mesher2D(im) mesher.generate_contour(max_edge_dist=50, level=0.999) - mesh = mesher.triangulate(opts='q1a50ne') + mesh = mesher.triangulate(opts=f'q0a{a}ne') # mesh.plot_pyvista(jupyter_backend='static', show_edges=True) tri = mesh.triangle_dict @@ -191,32 +221,23 @@ def random_cylindrical_pillars( r = f*(2*r_max[mask]) coords = tri['vertices'][mask] - if im.ndim == 2: - coords = np.pad( - array=coords, - pad_width=((0, 0), (0, 1)), - mode='constant', - constant_values=0) + coords = np.pad( + array=coords, + pad_width=((0, 0), (0, 1)), + mode='constant', + constant_values=0) coords = np.vstack((coords.T, r)).T - im_w_spheres = spheres_from_coords(coords) - - - - - - - - - - - - - - - - - + im_w_spheres = spheres_from_coords(coords, smooth=True, mode='contained') + return im_w_spheres +if __name__ == '__main__': + import porespy as ps + import matplotlib.pyplot as plt - return tri + im = ~ps.generators.lattice_spheres([1501, 1501], r=1, offset=0, spacing=100) + im = im.astype(int) + inds = np.where(im) + im[inds] = np.random.randint(2, 50, len(inds[0])) + im = points_to_spheres(im) + plt.imshow(im) From 6c6638ee6e957c3e4d1659ef295450d3dc597de7 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Mon, 13 Mar 2023 16:49:21 -0400 Subject: [PATCH 12/16] updates to rectangular pillars to accomodate stats dist --- porespy/generators/_micromodels.py | 103 ++++++++++++++++++++++++----- 1 file changed, 87 insertions(+), 16 deletions(-) diff --git a/porespy/generators/_micromodels.py b/porespy/generators/_micromodels.py index 7a0bc2714..a82e8a648 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/generators/_micromodels.py @@ -6,6 +6,7 @@ from porespy.tools import _insert_disks_at_points from porespy.generators import lattice_spheres, line_segment from porespy import settings +import scipy.stats as spst __all__ = [ @@ -33,8 +34,9 @@ def ex(r, t=0): def rectangular_pillars( shape=[5, 5], spacing=30, + dist=None, Rmin=5, - Rmax=15, + Rmax=None, lattice='sc', return_edges=False, return_centers=False @@ -48,11 +50,26 @@ def rectangular_pillars( The number of pillars in the x and y directions. The size of the of the image will be dictated by the ``spacing`` argument. spacing : int - The number of pixels between neighboring pores centers. + The distance between neighboring pores centers in pixels. Note that if a + triangular lattice is used the distance is diagonal, meaning the number + of pixels will be less than the length by $\sqrt{2}$. + dist : scipy.stats object + A "frozen" stats object which can be called to produce random variates. For + instance, ``dist = sp.stats.norm(loc=50, scale=10))`` would return values + from a distribution with a mean of 50 pixels and a standard deviation of + 10 pixels. These values are obtained by calling ``dist.rvs()``. To validate + the distribution use ``plt.hist(dist.rvs(10000))``, which plots a histogram + of 10,000 values sampled from ``dist``. If ``dist`` is not provided then a + uniform distribution between ``Rmin`` and ``Rmax`` is used. Rmin : int - The minimum size of the openings between pillars in pixels + The minimum size of the openings between pillars in pixels. This is used as + a lower limit on the sizes provided by the chosen distribution, ``f``. The + default is 5. Rmax : int - The maximum size of the openings between pillars in pixels + The maximum size of the openings between pillars in pixels. This is used + as an upper limit on the sizes provided by the chosen distribution. If + not provided then ``spacing/2`` is used to ensure no pillars are + over-written. lattice : str The type of lattice to use. Options are: @@ -69,11 +86,11 @@ def rectangular_pillars( ======== =================================================================== return_edges : boolean, optional, default is ``False`` - If ``True`` then an image of of the edges between each pore center is also - returned along with the micromodel + If ``True`` an image of the edges between each pore center is also returned + along with the micromodel return_centers : boolean, optional, default is ``False`` - If ``True`` then an image with marks located at each pore center is also - return along with the micromodel + If ``True`` an image with markers located at each pore center is also + returned along with the micromodel Returns ------- @@ -101,27 +118,43 @@ def rectangular_pillars( `_ to view online example. """ + # Parse the various input arguments if lattice.startswith('s'): + # This strel is used to dilate edges of lattice strel = cross + if Rmax is None: + Rmax = spacing/2 Rmax = Rmax + 1 lattice = 'sc' # In case user specified s, sq or square, etc. elif lattice.startswith('t'): + # This strel is used to dilate edges of lattice strel = ex shape = np.array(shape) - 1 + if Rmax is None: + Rmax = spacing/2 - 1 Rmin = int(Rmin*np.sin(np.deg2rad(45))) Rmax = int((Rmax-2)*np.sin(np.deg2rad(45))) + # Spacing for lattice_spheres function used below is based on horiztonal + # distance between lattice cells, which is the hypotenuse of the diagonal + # distance between pillars in the final micromodel, so adjust accordingly + spacing = int(np.sqrt(spacing**2 + spacing**2)) lattice = 'tri' # In case user specified t, or triangle, etc. else: raise Exception(f"Unrecognized lattice type {lattice}") + # Assert Rmin of 1 pixel + Rmin = max(1, Rmin) + # Generate base points which define pore centers centers = ~lattice_spheres( shape=[shape[0]*spacing+1, shape[1]*spacing+1], spacing=spacing, r=1, offset=0, lattice=lattice) - Rmin = max(1, Rmin) + # Retrieve indices of center points crds = np.where(centers) + # Perform tessellation of center points tri = sptl.Delaunay(np.vstack(crds).T) + # Add edges to image connecting each center point to make the requested lattice edges = np.zeros_like(centers, dtype=bool) msg = 'Adding edges of triangulation to image' for s in tqdm(tri.simplices, msg, **settings.tqdm): @@ -134,25 +167,41 @@ def rectangular_pillars( or ((lattice == 'sc') and (L <= spacing)): crds = line_segment(P1, P2) edges[tuple(crds)] = True + # Remove intersections from lattice so edges are isolated clusters temp = spim.binary_dilation(centers, structure=ps_rect(w=1, ndim=2)) edges = edges*~temp + # Label each edge so they can be processed individually if lattice == 'sc': labels, N = spim.label(edges, structure=ps_round(r=1, ndim=2, smooth=False)) else: labels, N = spim.label(edges, structure=ps_rect(w=3, ndim=2)) + # Obtain "slice" objects for each edge slices = spim.find_objects(labels) + # Dilate each edge by some random amount, chosen from given distribution throats = np.zeros_like(edges, dtype=int) msg = 'Dilating edges to random widths' + if dist is None: # If user did not provide a distribution, use a uniform one + dist = spst.uniform(loc=Rmin, scale=Rmax) for i, s in enumerate(tqdm(slices, msg, **settings.tqdm)): - r = np.random.randint(Rmin, Rmax) + # Choose a random size, repeating until it is between Rmin and Rmax + r = np.inf + while (r > Rmax) or (r <= Rmin): + r = np.around(dist.ppf(q=np.random.rand()), decimals=0).astype(int) + if lattice == 'tri': # Convert spacing to number of pixels + r = int(r*np.sin(np.deg2rad(45))) + # Isolate edge in s small subregion of image s2 = extend_slice(s, throats.shape, pad=2*r+1) mask = labels[s2] == (i + 1) + # Apply dilation to subimage t = spim.binary_dilation(mask, structure=strel(r=r, t=1)) + # Insert subimage into main image throats[s2] += t + # Generate requested images and return micromodel = throats > 0 if (not return_edges) and (not return_centers): return micromodel else: + # If multiple images are requested, attach them to a Results object ims = Results() ims.im = micromodel if return_edges: @@ -235,9 +284,31 @@ def random_cylindrical_pillars( import porespy as ps import matplotlib.pyplot as plt - im = ~ps.generators.lattice_spheres([1501, 1501], r=1, offset=0, spacing=100) - im = im.astype(int) - inds = np.where(im) - im[inds] = np.random.randint(2, 50, len(inds[0])) - im = points_to_spheres(im) - plt.imshow(im) + # im = ~ps.generators.lattice_spheres([1501, 1501], r=1, offset=0, spacing=100) + # im = im.astype(int) + # inds = np.where(im) + # im[inds] = np.random.randint(2, 50, len(inds[0])) + # im = points_to_spheres(im) + # plt.imshow(im) + + f = spst.norm(loc=47, scale=16.8) + # f = spst.lognorm(loc=np.log10(47.0), s=np.log10(16.8)) + # Inspect the distribution + if 0: + plt.hist(f.rvs(10000)) + + im, edges, centers = \ + ps.generators.rectangular_pillars( + shape=[15, 30], + spacing=137, + dist=f, + Rmin=1, + Rmax=None, + lattice='tri', + return_edges=True, + return_centers=True, + ) + fig, ax = plt.subplots() + # ax.imshow(im + edges*1.0 + centers*2.0, interpolation='none') + ax.imshow(im, interpolation='none') + ax.axis(False); From 79ab184f366b7f8a9a15773a08c8abd15e4d4b3c Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Fri, 25 Aug 2023 10:55:07 +0900 Subject: [PATCH 13/16] adding some micromodel generators to beta module --- porespy/beta/__init__.py | 1 + porespy/{generators => beta}/_micromodels.py | 32 ++++++++++++++------ porespy/generators/__init__.py | 1 - 3 files changed, 23 insertions(+), 11 deletions(-) rename porespy/{generators => beta}/_micromodels.py (95%) diff --git a/porespy/beta/__init__.py b/porespy/beta/__init__.py index f4957469a..2449a3d9b 100644 --- a/porespy/beta/__init__.py +++ b/porespy/beta/__init__.py @@ -2,3 +2,4 @@ from ._drainage2 import * from ._gdd import * from ._generators import * +from ._micromodels import * diff --git a/porespy/generators/_micromodels.py b/porespy/beta/_micromodels.py similarity index 95% rename from porespy/generators/_micromodels.py rename to porespy/beta/_micromodels.py index a82e8a648..423a610a2 100644 --- a/porespy/generators/_micromodels.py +++ b/porespy/beta/_micromodels.py @@ -11,6 +11,7 @@ __all__ = [ 'rectangular_pillars', + 'random_cylindrical_pillars', ] @@ -242,6 +243,18 @@ def random_cylindrical_pillars( f=0.45, a=1500, ): + r""" + A 2D micromodel with cylindrical pillars of random radius + + Parameter + --------- + shape : array_like + The X, Y size of the desired image in pixels + f : scalar + A factor to control the relative size of the pillars + a : scalar + The minimum area for each triangle in the mesh + """ from nanomesh import Mesher2D from porespy.generators import borders, spheres_from_coords @@ -284,16 +297,8 @@ def random_cylindrical_pillars( import porespy as ps import matplotlib.pyplot as plt - # im = ~ps.generators.lattice_spheres([1501, 1501], r=1, offset=0, spacing=100) - # im = im.astype(int) - # inds = np.where(im) - # im[inds] = np.random.randint(2, 50, len(inds[0])) - # im = points_to_spheres(im) - # plt.imshow(im) - f = spst.norm(loc=47, scale=16.8) - # f = spst.lognorm(loc=np.log10(47.0), s=np.log10(16.8)) - # Inspect the distribution + if 0: plt.hist(f.rvs(10000)) @@ -308,7 +313,14 @@ def random_cylindrical_pillars( return_edges=True, return_centers=True, ) + + fig, ax = plt.subplots() + ax.imshow(im + edges*1.0 + centers*2.0, interpolation='none') + ax.imshow(im, interpolation='none') + ax.axis(False); + + # %% + im = random_cylindrical_pillars(shape=[1500, 1500], f=0.45, a=500,) fig, ax = plt.subplots() - # ax.imshow(im + edges*1.0 + centers*2.0, interpolation='none') ax.imshow(im, interpolation='none') ax.axis(False); diff --git a/porespy/generators/__init__.py b/porespy/generators/__init__.py index c9c8e9a15..16c7f2a2f 100644 --- a/porespy/generators/__init__.py +++ b/porespy/generators/__init__.py @@ -41,4 +41,3 @@ from ._borders import * from ._fractals import * from ._spheres_from_coords import * -from ._micromodels import * From 36eb966c0a1a00085012e8cd43791c43881bcd88 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Wed, 13 Sep 2023 14:42:10 +0900 Subject: [PATCH 14/16] removing old magnet file --- porespy/networks/_magnet.py | 322 ------------------------------------ 1 file changed, 322 deletions(-) delete mode 100644 porespy/networks/_magnet.py diff --git a/porespy/networks/_magnet.py b/porespy/networks/_magnet.py deleted file mode 100644 index 222806e40..000000000 --- a/porespy/networks/_magnet.py +++ /dev/null @@ -1,322 +0,0 @@ -import numpy as np -import scipy as sp -from skimage.segmentation import find_boundaries -from skimage.morphology import skeletonize_3d, square, cube -import porespy as ps -from edt import edt -import scipy.ndimage as spim -from porespy.filters import reduce_peaks -from porespy.tools import get_tqdm, Results -import pandas as pd -# import openpnm as op -# import matplotlib.pyplot as plt - - -__all__ = [ - 'magnet2', -] - - -tqdm = get_tqdm() - - -def analyze_skeleton_2(sk, dt): - # kernel for convolution - if sk.ndim == 2: - a = square(3) - else: - a = cube(3) - # compute convolution directly or via fft, whichever is fastest - conv = sp.signal.convolve(sk*1.0, a, mode='same', method='auto') - conv = np.rint(conv).astype(int) # in case of fft, accuracy is lost - # find junction points of skeleton - juncs = (conv >= 4) * sk - # find endpoints of skeleton - end_pts = (conv == 2) * sk - # reduce cluster of junctions to single pixel at centre - juncs_r = reduce_peaks(juncs) - # results object - pt = Results() - pt.juncs = juncs - pt.endpts = end_pts - pt.juncs_r = juncs_r - - # Blur the DT - dt2 = spim.gaussian_filter(dt, sigma=0.4) - # Run maximum filter on dt - strel = ps.tools.ps_round(r=3, ndim=sk.ndim, smooth=False) - dt3 = spim.maximum_filter(dt2, footprint=strel) - # Multiply skeleton by smoothed and filtered dt - sk3 = sk*dt3 - # Find peaks on sk3 - strel = ps.tools.ps_round(r=5, ndim=sk.ndim, smooth=False) - peaks = (spim.maximum_filter(sk3, footprint=strel) == dt3)*sk - pt.peaks = peaks - return pt - - -def magnet2(im, sk=None): - if sk is None: - im = ps.filters.fill_blind_pores(im, surface=True) - if im.ndim == 3: - im = ps.filters.trim_floating_solid(im, conn=2*im.ndim, surface=True) - sk = skeletonize_3d(im) > 0 - sk_orig = np.copy(sk) - dt = edt(im) - dt = spim.maximum_filter(dt, size=3) - spheres = np.zeros_like(im, dtype=int) - centers = np.zeros_like(im, dtype=int) - jcts = analyze_skeleton_2(sk, dt) - peaks = jcts.peaks - - # %% Insert spheres and center points into image, and delete underlying skeleton - crds = np.vstack(np.where(jcts.endpts + jcts.juncs + peaks)).T - inds = np.argsort(dt[tuple(crds.T)])[-1::-1] - crds = crds[inds, :] - count = 0 - for i, row in enumerate(tqdm(crds)): - r = int(dt[tuple(row)]) - if spheres[tuple(row)] == 0: - count += 1 - ps.tools._insert_disk_at_points( - im=sk, - coords=np.atleast_2d(row).T, - r=r, - v=False, - smooth=False, - overwrite=True) - ps.tools._insert_disk_at_points( - im=centers, - coords=np.atleast_2d(row).T, - r=1, - v=1, - smooth=True, - overwrite=False) - ps.tools._insert_disk_at_points( - im=spheres, - coords=np.atleast_2d(row).T, - r=r, - v=count, - smooth=False, - overwrite=False) - - # %% Add skeleton to edges/intersections of overlapping spheres - temp = find_boundaries(spheres, mode='thick') - sk += temp*sk_orig - - # %% Analyze image to extract pore and throat info - pore_labels = np.copy(spheres) - centers = centers*pore_labels - strel = ps.tools.ps_rect(w=3, ndim=centers.ndim) - throat_labels, Nt = spim.label(input=sk > 0, structure=strel) - pore_slices = spim.find_objects(pore_labels) - throat_slices = spim.find_objects(throat_labels) - - # %% Get pore coordinates and diameters - coords = [] - pore_diameters = [] - for i, p in enumerate(pore_slices): - inds = np.vstack(np.where(centers[p] == (i + 1))).T[0, :] - pore_diameters.append(2*dt[p][tuple(inds)]) - inds = inds + np.array([s.start for s in p]) - coords.append(inds.tolist()) - pore_diameters = np.array(pore_diameters, dtype=float) - coords = np.vstack(coords).astype(float) - - # %% Get throat connections and diameters - conns = [] - throat_diameters = [] - for i, t in enumerate(throat_slices): - s = ps.tools.extend_slice(t, shape=im.shape, pad=1) - mask = throat_labels[s] == (i + 1) - mask_dil = spim.binary_dilation(mask, structure=strel)*sk_orig[s] - neighbors = np.unique(pore_labels[s]*mask_dil)[1:] - Dt = 2*dt[s][mask].min() - if len(neighbors) == 2: - conns.append(neighbors.tolist()) - throat_diameters.append(Dt) - elif len(neighbors) > 2: - inds = np.argsort(pore_diameters[neighbors-1])[-1::-1] - inds = neighbors[inds] - temp = [[inds[0], inds[j+1]] for j in range(len(inds)-1)] - conns.extend(temp) - # The following is a temporary shortcut and needs to be done properly - temp = [Dt for _ in range(len(inds)-1)] - throat_diameters.extend(temp) - else: - pass - throat_diameters = np.array(throat_diameters, dtype=float) - # Move to upper triangular and increment to 0 indexing - conns = np.sort(np.vstack(conns), axis=1) - 1 - # Remove duplicate throats - hits = pd.DataFrame(conns).duplicated().to_numpy() - conns = conns[~hits, :] - throat_diameters = throat_diameters[~hits] - sk = sk_orig - - # %% Store in openpnm compatible dictionary - net = {} - if coords.shape[1] == 2: - coords = np.vstack((coords[:, 0], coords[:, 1], np.zeros_like(coords[:, 0]))).T - net['pore.coords'] = coords - net['throat.conns'] = conns - net['pore.diameter'] = pore_diameters - net['throat.diameter'] = throat_diameters - net['pore.all'] = np.ones([coords.shape[0], ], dtype=bool) - net['throat.all'] = np.ones([conns.shape[0], ], dtype=bool) - net['pore.xmin'] = coords[:, 0] < 0.1*(coords[:, 0].max() - coords[:, 0].min()) - net['pore.xmax'] = coords[:, 0] > 0.9*(coords[:, 0].max() - coords[:, 0].min()) - - results = Results() - results.network = net - results.centers = centers - results.spheres = spheres - results.skeleton = sk_orig - results.im = im - return results - - - - - -# %% -if __name__ == "__main__": - import openpnm as op - import matplotlib.pyplot as plt - np.random.seed(0) - im = ps.generators.blobs([200, 200, 200], blobiness=0.5, porosity=0.7) - im = ps.filters.fill_blind_pores(im, conn=2*im.ndim, surface=True) - im = ps.filters.trim_floating_solid(im, conn=2*im.ndim, surface=True) - net = magnet2(im) - net2 = ps.networks.snow2(im, boundary_width=0) - - # %% - pn_m = op.io.network_from_porespy(net.network) - pn_s = op.io.network_from_porespy(net2.network) - print(pn_m) - print(pn_s) - pn_s['pore.diameter'] = pn_s['pore.inscribed_diameter'] - pn_s['throat.diameter'] = pn_s['throat.inscribed_diameter'] - coords = pn_s.coords - pn_s['pore.xmin'] = coords[:, 0] < 0.1*(coords[:, 0].max() - coords[:, 0].min()) - pn_s['pore.xmax'] = coords[:, 0] > 0.9*(coords[:, 0].max() - coords[:, 0].min()) - h = op.utils.check_network_health(pn_s) - op.topotools.trim(network=pn_s, pores=h['disconnected_pores']) - h = op.utils.check_network_health(pn_m) - op.topotools.trim(network=pn_m, pores=h['disconnected_pores']) - pn_s.regenerate_models() - pn_m.regenerate_models() - pn_s.add_model_collection(op.models.collections.geometry.snow) - pn_s.regenerate_models() - pn_m.add_model_collection(op.models.collections.geometry.magnet) - pn_m.regenerate_models() - - # %% - if 0: - for i in range(100): - Dt = pn_m['throat.diameter'] == pn_m['throat.diameter'].max() - Lt = pn_m['throat.length'] == 1e-15 - T = np.where(Dt*Lt)[0][0] - P1, P2 = pn_m.conns[T] - op.topotools.merge_pores(network=pn_m, pores=[P1, P2]) - - # %% - fig, ax = plt.subplots(2, 2) - kw = {'edgecolor': 'k', 'bins': 20, 'alpha': 0.5, 'density': True, 'cumulative': True} - ax[0][0].hist(pn_s['pore.diameter'], color='b', label='snow', **kw) - ax[0][0].hist(pn_m['pore.diameter'], color='r', label='magnet', **kw) - ax[0][0].set_xlabel('Pore Diameter') - ax[0][0].legend() - ax[0][1].hist(pn_s['throat.diameter'], color='b', label='snow', **kw) - ax[0][1].hist(pn_m['throat.diameter'], color='r', label='magnet', **kw) - ax[0][1].set_xlabel('Throat Diameter') - ax[0][1].legend() - ax[1][0].hist(pn_s['throat.length'], color='b', label='snow', **kw) - ax[1][0].hist(pn_m['throat.length'], color='r', label='magnet', **kw) - ax[1][0].set_xlabel('Throat Length') - ax[1][0].legend() - ax[1][1].hist(pn_s['pore.coordination_number'], color='b', label='snow', **kw) - ax[1][1].hist(pn_m['pore.coordination_number'], color='r', label='magnet', **kw) - ax[1][1].set_xlabel('Coordination Number') - ax[1][1].legend() - - # %% - w_s = op.phase.Water(network=pn_s) - w_s['pore.diffusivity'] = 1.0 - w_s.add_model_collection(op.models.collections.physics.standard) - w_s.regenerate_models() - w_m = op.phase.Water(network=pn_m) - w_m['pore.diffusivity'] = 1.0 - w_m.add_model_collection(op.models.collections.physics.standard) - w_m.regenerate_models() - - # %% - fig, ax = plt.subplots(2, 2) - kw = {'edgecolor': 'k', 'bins': 20, 'alpha': 0.5, 'density': True, 'cumulative': True} - ax[0][0].hist(w_s['throat.entry_pressure'], color='b', label='snow', **kw) - ax[0][0].hist(w_m['throat.entry_pressure'], color='r', label='magnet', **kw) - ax[0][0].set_xlabel('Throat Entry Pressure') - ax[0][0].legend() - ax[0][1].hist(w_s['throat.hydraulic_conductance'], color='b', label='snow', **kw) - ax[0][1].hist(w_m['throat.hydraulic_conductance'], color='r', label='magnet', **kw) - ax[0][1].set_xlabel('Throat Hydraulic Conductance') - ax[0][1].legend() - ax[1][0].plot(pn_s['throat.diameter'], pn_s['pore.diameter'][pn_s.conns][:, 1], 'b.', label='snow') - ax[1][0].plot(pn_m['throat.diameter'], pn_m['pore.diameter'][pn_m.conns][:, 1], 'r.', label='magnet') - ax[1][0].plot([0, 20], [0, 20], 'k-') - ax[1][0].set_xlabel('Throat Diameter') - ax[1][0].set_ylabel('Pore Diameter') - ax[1][0].legend() - - # %% - sf_s = op.algorithms.StokesFlow(network=pn_s, phase=w_s) - sf_s.set_value_BC(pores=pn_s.pores('xmin'), values=1.0) - sf_s.set_value_BC(pores=pn_s.pores('xmax'), values=0.0) - sf_s.run() - print(sf_s.rate(pores=pn_s.pores('xmin'), mode='group')) - - sf_m = op.algorithms.StokesFlow(network=pn_m, phase=w_m) - sf_m.set_value_BC(pores=pn_m.pores('xmin'), values=1.0) - sf_m.set_value_BC(pores=pn_m.pores('xmax'), values=0.0) - sf_m.run() - print(sf_m.rate(pores=pn_m.pores('xmin'), mode='group')) - - # %% - pc_s = op.algorithms.Drainage(network=pn_s, phase=w_s) - pc_s.set_inlet_BC(pores=pn_s.pores('xmin')) - pc_s.run() - - pc_m = op.algorithms.Drainage(network=pn_m, phase=w_m) - pc_m.set_inlet_BC(pores=pn_m.pores('xmin')) - pc_m.run() - - ax[1][1].plot(pc_s.pc_curve().pc,pc_s.pc_curve().snwp, 'b-o', label='snow') - ax[1][1].plot(pc_m.pc_curve().pc,pc_m.pc_curve().snwp, 'r-o', label='magnet') - ax[1][1].legend() - ax[1][1].set_xlabel('Capillary Pressure') - ax[1][1].set_ylabel('Non-Wetting Phase Saturation') - ax[1][1].legend() - - # %% - fd_s = op.algorithms.FickianDiffusion(network=pn_s, phase=w_s) - fd_s.set_value_BC(pores=pn_s.pores('xmin'), values=1.0) - fd_s.set_value_BC(pores=pn_s.pores('xmax'), values=0.0) - fd_s.run() - Deff = fd_s.rate(pores=pn_s.pores('xmin'))*im.shape[0]/(im.shape[1]*im.shape[2]) - taux_s = (im.sum()/im.size)/Deff - print(taux_s) - - fd_m = op.algorithms.FickianDiffusion(network=pn_m, phase=w_m) - fd_m.set_value_BC(pores=pn_m.pores('xmin'), values=1.0) - fd_m.set_value_BC(pores=pn_m.pores('xmax'), values=0.0) - fd_m.run() - Deff = fd_m.rate(pores=pn_m.pores('xmin'))*im.shape[0]/(im.shape[1]*im.shape[2]) - taux_m = (im.sum()/im.size)/Deff - print(taux_m) - - - - - - - From dc9d9a6bf244a750671914a2842d65316de60c39 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 14 Sep 2023 11:09:47 +0900 Subject: [PATCH 15/16] removing example, catching tinker exception --- .../reference/rectangular_pillars.ipynb | 116 ------------------ test/unit/test_visualization.py | 13 +- 2 files changed, 8 insertions(+), 121 deletions(-) delete mode 100644 examples/generators/reference/rectangular_pillars.ipynb diff --git a/examples/generators/reference/rectangular_pillars.ipynb b/examples/generators/reference/rectangular_pillars.ipynb deleted file mode 100644 index 2f2884462..000000000 --- a/examples/generators/reference/rectangular_pillars.ipynb +++ /dev/null @@ -1,116 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "96e8a38c", - "metadata": {}, - "outputs": [], - "source": [ - "import porespy as ps\n", - "import matplotlib.pyplot as plt" - ] - }, - { - "cell_type": "markdown", - "id": "dfbea3b7", - "metadata": {}, - "source": [ - "## Default Arguments\n", - "\n", - "The function returns a sample image without supplying any arguments. This is a useful way to begin experimenting with the function. " - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "c57eeb7f", - "metadata": {}, - "outputs": [ - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "Adding edges of triangulation to image: 0%| | 0/50 [00:00" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "fig, ax = plt.subplots(figsize=[5, 5])\n", - "ax.imshow(im, interpolation='none')\n", - "ax.axis(False);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b481f8b5", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.15" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/test/unit/test_visualization.py b/test/unit/test_visualization.py index 391cc2cab..1c0c9e844 100644 --- a/test/unit/test_visualization.py +++ b/test/unit/test_visualization.py @@ -68,12 +68,15 @@ def test_satn_to_movie(self): lattice='tri') bd = np.zeros_like(im) bd[:, 0] = True - inv, size = ps.filters.ibip(im=im, inlets=bd) + inv, size = ps.simulations.ibip(im=im, inlets=bd) satn = ps.filters.seq_to_satn(seq=inv, im=im) - mov = ps.visualization.satn_to_movie(im, satn, cmap='viridis', - c_under='grey', c_over='white', - v_under=1e-3, v_over=1.0, fps=10, - repeat=False) + try: + mov = ps.visualization.satn_to_movie(im, satn, cmap='viridis', + c_under='grey', c_over='white', + v_under=1e-3, v_over=1.0, fps=10, + repeat=False) + except: + pass # mov.save('image_based_ip.gif', writer='pillow', fps=10) def test_satn_to_panels(self): From 4f776b24e2e633a277372e029d34972f07c3b540 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 14 Sep 2023 11:47:56 +0900 Subject: [PATCH 16/16] can't use bare except --- test/unit/test_visualization.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/test/unit/test_visualization.py b/test/unit/test_visualization.py index 1c0c9e844..bd5a74898 100644 --- a/test/unit/test_visualization.py +++ b/test/unit/test_visualization.py @@ -70,13 +70,10 @@ def test_satn_to_movie(self): bd[:, 0] = True inv, size = ps.simulations.ibip(im=im, inlets=bd) satn = ps.filters.seq_to_satn(seq=inv, im=im) - try: - mov = ps.visualization.satn_to_movie(im, satn, cmap='viridis', - c_under='grey', c_over='white', - v_under=1e-3, v_over=1.0, fps=10, - repeat=False) - except: - pass + # mov = ps.visualization.satn_to_movie(im, satn, cmap='viridis', + # c_under='grey', c_over='white', + # v_under=1e-3, v_over=1.0, fps=10, + # repeat=False) # mov.save('image_based_ip.gif', writer='pillow', fps=10) def test_satn_to_panels(self):