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/beta/_micromodels.py b/porespy/beta/_micromodels.py new file mode 100644 index 000000000..423a610a2 --- /dev/null +++ b/porespy/beta/_micromodels.py @@ -0,0 +1,326 @@ +# 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, get_tqdm, Results +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__ = [ + 'rectangular_pillars', + 'random_cylindrical_pillars', +] + + +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 + 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, + dist=None, + Rmin=5, + Rmax=None, + lattice='sc', + return_edges=False, + return_centers=False +): + 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 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. 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. 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: + + ======== =================================================================== + 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. + ======== =================================================================== + + return_edges : boolean, optional, default is ``False`` + 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`` an image with markers located at each pore center is also + returned along with the micromodel + + Returns + ------- + 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 + ========== ================================================================= + 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. + """ + # 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) + # 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): + 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 + # 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)): + # 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: + ims.edges = edges + if return_centers: + ims.centers = centers + 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=[1500, 1500], + 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 + + 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 + + mesher = Mesher2D(im) + mesher.generate_contour(max_edge_dist=50, level=0.999) + + mesh = mesher.triangulate(opts=f'q0a{a}ne') + # 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 = tri['vertices'][mask] + 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, smooth=True, mode='contained') + return im_w_spheres + + +if __name__ == '__main__': + import porespy as ps + import matplotlib.pyplot as plt + + f = spst.norm(loc=47, scale=16.8) + + 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); + + # %% + im = random_cylindrical_pillars(shape=[1500, 1500], f=0.45, a=500,) + fig, ax = plt.subplots() + ax.imshow(im, interpolation='none') + ax.axis(False); diff --git a/test/unit/test_visualization.py b/test/unit/test_visualization.py index 391cc2cab..bd5a74898 100644 --- a/test/unit/test_visualization.py +++ b/test/unit/test_visualization.py @@ -68,12 +68,12 @@ 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) + # 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):