1 1
import sys
2 1
import numpy as np
3 1
from edt import edt
4 1
import porespy as ps
5 1
from numba import njit
6 1
import scipy.spatial as sptl
7 1
import scipy.ndimage as spim
8 1
from porespy.tools import norm_to_uniform, ps_ball, ps_disk, get_border
9 1
from typing import List
10 1
from numpy import array
11 1
from tqdm import tqdm
12

13

14 1
def insert_shape(im, element, center=None, corner=None, value=1, mode="overwrite"):
15
    r"""
16
    Inserts sub-image into a larger image at the specified location.
17

18
    If the inserted image extends beyond the boundaries of the image it will
19
    be cropped accordingly.
20

21
    Parameters
22
    ----------
23
    im : ND-array
24
        The image into which the sub-image will be inserted
25
    element : ND-array
26
        The sub-image to insert
27
    center : tuple
28
        Coordinates indicating the position in the main image where the
29
        inserted imaged will be centered.  If ``center`` is given then
30
        ``corner`` cannot be specified.  Note that ``center`` can only be
31
        used if all dimensions of ``element`` are odd, otherwise the meaning
32
        of center is not defined.
33
    corner : tuple
34
        Coordinates indicating the position in the main image where the
35
        lower corner (i.e. [0, 0, 0]) of the inserted image should be anchored.
36
        If ``corner`` is given then ``corner`` cannot be specified.
37
    value : scalar
38
        A scalar value to apply to the sub-image.  The default is 1.
39
    mode : string
40
        If 'overwrite' (default) the inserted image replaces the values in the
41
        main image.  If 'overlay' the inserted image is added to the main
42
        image.  In both cases the inserted image is multiplied by ``value``
43
        first.
44

45
    Returns
46
    -------
47
    im : ND-array
48
        A copy of ``im`` with the supplied element inserted.
49

50
    """
51 1
    im = im.copy()
52 1
    if im.ndim != element.ndim:
53 0
        raise Exception(
54
            f"Image shape {im.shape} and element shape {element.shape} do not match"
55
        )
56 1
    s_im = []
57 1
    s_el = []
58 1
    if (center is not None) and (corner is None):
59 1
        for dim in range(im.ndim):
60 1
            r, d = np.divmod(element.shape[dim], 2)
61 1
            if d == 0:
62 1
                raise Exception(
63
                    "Cannot specify center point when element "
64
                    + "has one or more even dimension"
65
                )
66 1
            lower_im = np.amax((center[dim] - r, 0))
67 1
            upper_im = np.amin((center[dim] + r + 1, im.shape[dim]))
68 1
            s_im.append(slice(lower_im, upper_im))
69 1
            lower_el = np.amax((lower_im - center[dim] + r, 0))
70 1
            upper_el = np.amin((upper_im - center[dim] + r, element.shape[dim]))
71 1
            s_el.append(slice(lower_el, upper_el))
72 1
    elif (corner is not None) and (center is None):
73 1
        for dim in range(im.ndim):
74 1
            L = int(element.shape[dim])
75 1
            lower_im = np.amax((corner[dim], 0))
76 1
            upper_im = np.amin((corner[dim] + L, im.shape[dim]))
77 1
            s_im.append(slice(lower_im, upper_im))
78 1
            lower_el = np.amax((lower_im - corner[dim], 0))
79 1
            upper_el = np.amin((upper_im - corner[dim], element.shape[dim]))
80 1
            s_el.append(slice(min(lower_el, upper_el), upper_el))
81
    else:
82 0
        raise Exception("Cannot specify both corner and center")
83

84 1
    if mode == "overlay":
85 1
        im[tuple(s_im)] = im[tuple(s_im)] + element[tuple(s_el)] * value
86 1
    elif mode == "overwrite":
87 1
        im[tuple(s_im)] = element[tuple(s_el)] * value
88
    else:
89 0
        raise Exception("Invalid mode " + mode)
90 1
    return im
91

92

93 1
def RSA(im: array,
94
        radius: int,
95
        volume_fraction: int = 1,
96
        n_max: int = None,
97
        mode: str = "contained"):
98
    r"""
99
    Generates a sphere or disk packing using Random Sequential Addition
100

101
    This algorithm ensures that spheres do not overlap but does not
102
    guarantee they are tightly packed.
103

104
    This function adds spheres to the background of the received ``im``, which
105
    allows iteratively adding spheres of different radii to the unfilled space,
106
    be repeatedly passing in the result of previous calls to RSA.
107

108
    Parameters
109
    ----------
110
    im : ND-array
111
        The image into which the spheres should be inserted.  By accepting an
112
        image rather than a shape, it allows users to insert spheres into an
113
        already existing image.  To begin the process, start with an array of
114
        zeros such as ``im = np.zeros([200, 200, 200], dtype=bool)``.
115
    radius : int
116
        The radius of the disk or sphere to insert.
117
    volume_fraction : scalar (default is 1.0)
118
        The fraction of the image that should be filled with spheres.  The
119
        spheres are added as 1's, so each sphere addition increases the
120
        ``volume_fraction`` until the specified limit is reach.  Note that if
121
        ``n_max`` is reached first, then ``volume_fraction`` will not be
122
        acheived.
123
    n_max : int (default is 10,000)
124
        The maximum number of spheres to add.  By default the value of
125
        ``n_max`` is high so that the addition of spheres will go indefinately
126
        until ``volume_fraction`` is met, but specifying a smaller value
127
        will halt addition after the given number of spheres are added.
128
    mode : string (default is 'contained')
129
        Controls how the edges of the image are handled.  Options are:
130

131
        'contained' - Spheres are all completely within the image
132

133
        'extended' - Spheres are allowed to extend beyond the edge of the
134
        image.  In this mode the volume fraction will be less that requested
135
        since some spheres extend beyond the image, but their entire volume
136
        is counted as added for computational efficiency.
137

138
    Returns
139
    -------
140
    image : ND-array
141
        A handle to the input ``im`` with spheres of specified radius
142
        *added* to the background.
143

144
    Notes
145
    -----
146
    This function uses Numba to speed up the search for valid sphere insertion
147
    points.  It seems that Numba does not look at the state of the scipy
148
    random number generator, so setting the seed to a known value has no
149
    effect on the output of this function. Each call to this function will
150
    produce a unique image.  If you wish to use the same realization multiple
151
    times you must save the array (e.g. ``numpy.save``).
152

153
    References
154
    ----------
155
    [1] Random Heterogeneous Materials, S. Torquato (2001)
156

157
    """
158 1
    print(80 * "-")
159 1
    print(f"RSA: Adding spheres of size {radius}")
160 1
    im = im.astype(bool)
161 1
    if n_max is None:
162 1
        n_max = 10000
163 1
    vf_final = volume_fraction
164 1
    vf_start = im.sum() / im.size
165 1
    print("Initial volume fraction:", vf_start)
166 1
    if im.ndim == 2:
167 1
        template_lg = ps_disk(radius * 2)
168 1
        template_sm = ps_disk(radius)
169
    else:
170 1
        template_lg = ps_ball(radius * 2)
171 1
        template_sm = ps_ball(radius)
172 1
    vf_template = template_sm.sum() / im.size
173
    # Pad image by the radius of large template to enable insertion near edges
174 1
    im = np.pad(im, pad_width=2 * radius, mode="edge")
175
    # Depending on mode, adjust mask to remove options around edge
176 1
    if mode == "contained":
177 1
        border = get_border(im.shape, thickness=2 * radius, mode="faces")
178 1
    elif mode == "extended":
179 1
        border = get_border(im.shape, thickness=radius + 1, mode="faces")
180
    else:
181 0
        raise Exception("Unrecognized mode: ", mode)
182
    # Remove border pixels
183 1
    im[border] = True
184
    # Dilate existing objects by strel to remove pixels near them
185
    # from consideration for sphere placement
186 1
    print("Dilating foreground features by sphere radius")
187 1
    dt = edt(im == 0)
188 1
    options_im = dt >= radius
189
    # ------------------------------------------------------------------------
190
    # Begin inserting the spheres
191 1
    vf = vf_start
192 1
    free_sites = np.flatnonzero(options_im)
193 1
    i = 0
194 1
    while (vf <= vf_final) and (i < n_max) and (len(free_sites) > 0):
195 1
        c, count = _make_choice(options_im, free_sites=free_sites)
196
        # The 100 below is arbitrary and may change performance
197 1
        if count > 100:
198
            # Regenerate list of free_sites
199 1
            print("Regenerating free_sites after", i, "iterations")
200 1
            free_sites = np.flatnonzero(options_im)
201 1
        if all(np.array(c) == -1):
202 1
            break
203 1
        s_sm = tuple([slice(x - radius, x + radius + 1, None) for x in c])
204 1
        s_lg = tuple([slice(x - 2 * radius, x + 2 * radius + 1, None) for x in c])
205 1
        im[s_sm] += template_sm  # Add ball to image
206 1
        options_im[s_lg][template_lg] = False  # Update extended region
207 1
        vf += vf_template
208 1
        i += 1
209 1
    print("Number of spheres inserted:", i)
210
    # ------------------------------------------------------------------------
211
    # Get slice into returned image to retain original size
212 1
    s = tuple([slice(2 * radius, d - 2 * radius, None) for d in im.shape])
213 1
    im = im[s]
214 1
    vf = im.sum() / im.size
215 1
    print("Final volume fraction:", vf)
216 1
    return im
217

218

219 1
@njit
220
def _make_choice(options_im, free_sites):
221
    r"""
222
    This function is called by _begin_inserting to find valid insertion points
223

224
    Parameters
225
    ----------
226
    options_im : ND-array
227
        An array with ``True`` at all valid locations and ``False`` at all
228
        locations where a sphere already exists PLUS a region of radius R
229
        around each sphere since these points are also invalid insertion
230
        points.
231
    free_sites : array_like
232
        A 1D array containing valid insertion indices.  This list is used to
233
        select insertion points from a limited which occasionally gets
234
        smaller.
235

236
    Returns
237
    -------
238
    coords : list
239
        The XY or XYZ coordinates of the next insertion point
240
    count : int
241
        The number of attempts that were needed to find valid point.  If
242
        this value gets too high, a short list of ``free_sites`` should be
243
        generated in the calling function.
244

245
    """
246 0
    choice = False
247 0
    count = 0
248 0
    upper_limit = len(free_sites)
249 0
    max_iters = upper_limit * 20
250 0
    if options_im.ndim == 2:
251 0
        coords = [-1, -1]
252 0
        Nx, Ny = options_im.shape
253 0
        while not choice:
254 0
            if count >= max_iters:
255 0
                coords = [-1, -1]
256 0
                break
257 0
            ind = np.random.randint(0, upper_limit)
258
            # This numpy function is not supported by numba yet
259
            # c1, c2 = np.unravel_index(free_sites[ind], options_im.shape)
260
            # So using manual unraveling
261 0
            coords[1] = free_sites[ind] % Ny
262 0
            coords[0] = (free_sites[ind] // Ny) % Nx
263 0
            choice = options_im[coords[0], coords[1]]
264 0
            count += 1
265 0
    if options_im.ndim == 3:
266 0
        coords = [-1, -1, -1]
267 0
        Nx, Ny, Nz = options_im.shape
268 0
        while not choice:
269 0
            if count >= max_iters:
270 0
                coords = [-1, -1, -1]
271 0
                break
272 0
            ind = np.random.randint(0, upper_limit)
273
            # This numpy function is not supported by numba yet
274
            # c1, c2, c3 = np.unravel_index(free_sites[ind], options_im.shape)
275
            # So using manual unraveling
276 0
            coords[2] = free_sites[ind] % Nz
277 0
            coords[1] = (free_sites[ind] // Nz) % Ny
278 0
            coords[0] = (free_sites[ind] // (Nz * Ny)) % Nx
279 0
            choice = options_im[coords[0], coords[1], coords[2]]
280 0
            count += 1
281 0
    return coords, count
282

283

284 1
def bundle_of_tubes(shape: List[int], spacing: int):
285
    r"""
286
    Create a 3D image of a bundle of tubes, in the form of a rectangular
287
    plate with randomly sized holes through it.
288

289
    Parameters
290
    ----------
291
    shape : list
292
        The size the image, with the 3rd dimension indicating the plate
293
        thickness.  If the 3rd dimension is not given then a thickness of
294
        1 voxel is assumed.
295

296
    spacing : scalar
297
        The center to center distance of the holes.  The hole sizes will be
298
        randomly distributed between this values down to 3 voxels.
299

300
    Returns
301
    -------
302
    image : ND-array
303
        A boolean array with ``True`` values denoting the pore space
304
    """
305 1
    shape = np.array(shape)
306 1
    if np.size(shape) == 1:
307 0
        shape = np.full((3,), int(shape))
308 1
    if np.size(shape) == 2:
309 0
        shape = np.hstack((shape, [1]))
310 1
    temp = np.zeros(shape=shape[:2])
311 1
    Xi = np.ceil(
312
        np.linspace(spacing / 2, shape[0] - (spacing / 2) - 1, int(shape[0] / spacing))
313
    )
314 1
    Xi = np.array(Xi, dtype=int)
315 1
    Yi = np.ceil(
316
        np.linspace(spacing / 2, shape[1] - (spacing / 2) - 1, int(shape[1] / spacing))
317
    )
318 1
    Yi = np.array(Yi, dtype=int)
319 1
    temp[tuple(np.meshgrid(Xi, Yi))] = 1
320 1
    inds = np.where(temp)
321 1
    for i in range(len(inds[0])):
322 1
        r = np.random.randint(1, (spacing / 2))
323 1
        try:
324 1
            s1 = slice(inds[0][i] - r, inds[0][i] + r + 1)
325 1
            s2 = slice(inds[1][i] - r, inds[1][i] + r + 1)
326 1
            temp[s1, s2] = ps_disk(r)
327 0
        except ValueError:
328 0
            odd_shape = np.shape(temp[s1, s2])
329 0
            temp[s1, s2] = ps_disk(r)[: odd_shape[0], : odd_shape[1]]
330 1
    im = np.broadcast_to(array=np.atleast_3d(temp), shape=shape)
331 1
    return im
332

333

334 1
def polydisperse_spheres(
335
    shape: List[int], porosity: float, dist, nbins: int = 5, r_min: int = 5
336
):
337
    r"""
338
    Create an image of randomly place, overlapping spheres with a distribution
339
    of radii.
340

341
    Parameters
342
    ----------
343
    shape : list
344
        The size of the image to generate in [Nx, Ny, Nz] where Ni is the
345
        number of voxels in each direction.  If shape is only 2D, then an
346
        image of polydisperse disks is returns
347

348
    porosity : scalar
349
        The porosity of the image, defined as the number of void voxels
350
        divided by the number of voxels in the image. The specified value
351
        is only matched approximately, so it's suggested to check this value
352
        after the image is generated.
353

354
    dist : scipy.stats distribution object
355
        This should be an initialized distribution chosen from the large number
356
        of options in the ``scipy.stats`` submodule.  For instance, a normal
357
        distribution with a mean of 20 and a standard deviation of 10 can be
358
        obtained with ``dist = scipy.stats.norm(loc=20, scale=10)``
359

360
    nbins : scalar
361
        The number of discrete sphere sizes that will be used to generate the
362
        image.  This function generates  ``nbins`` images of monodisperse
363
        spheres that span 0.05 and 0.95 of the possible values produced by the
364
        provided distribution, then overlays them to get polydispersivity.
365

366
    Returns
367
    -------
368
    image : ND-array
369
        A boolean array with ``True`` values denoting the pore space
370
    """
371 1
    shape = np.array(shape)
372 1
    if np.size(shape) == 1:
373 0
        shape = np.full((3,), int(shape))
374 1
    Rs = dist.interval(np.linspace(0.05, 0.95, nbins))
375 1
    Rs = np.vstack(Rs).T
376 1
    Rs = (Rs[:-1] + Rs[1:]) / 2
377 1
    Rs = np.clip(Rs.flatten(), a_min=r_min, a_max=None)
378 1
    phi_desired = 1 - (1 - porosity) / (len(Rs))
379 1
    im = np.ones(shape, dtype=bool)
380 1
    for r in Rs:
381 1
        phi_im = im.sum() / np.prod(shape)
382 1
        phi_corrected = 1 - (1 - phi_desired) / phi_im
383 1
        temp = overlapping_spheres(shape=shape, radius=r, porosity=phi_corrected)
384 1
        im = im * temp
385 1
    return im
386

387

388 1
def voronoi_edges(shape: List[int], radius: int, ncells: int, flat_faces: bool = True):
389
    r"""
390
    Create an image of the edges in a Voronoi tessellation
391

392
    Parameters
393
    ----------
394
    shape : array_like
395
        The size of the image to generate in [Nx, Ny, Nz] where Ni is the
396
        number of voxels in each direction.
397

398
    radius : scalar
399
        The radius to which Voronoi edges should be dilated in the final image.
400

401
    ncells : scalar
402
        The number of Voronoi cells to include in the tesselation.
403

404
    flat_faces : Boolean
405
        Whether the Voronoi edges should lie on the boundary of the
406
        image (True), or if edges outside the image should be removed (False).
407

408
    Returns
409
    -------
410
    image : ND-array
411
        A boolean array with ``True`` values denoting the pore space
412

413
    """
414 1
    print(60 * "-")
415 1
    print("voronoi_edges: Generating", ncells, "cells")
416 1
    shape = np.array(shape)
417 1
    if np.size(shape) == 1:
418 0
        shape = np.full((3,), int(shape))
419 1
    im = np.zeros(shape, dtype=bool)
420 1
    base_pts = np.random.rand(ncells, 3) * shape
421 1
    if flat_faces:
422
        # Reflect base points
423 1
        Nx, Ny, Nz = shape
424 1
        orig_pts = base_pts
425 1
        base_pts = np.vstack((base_pts, [-1, 1, 1] * orig_pts + [2.0 * Nx, 0, 0]))
426 1
        base_pts = np.vstack((base_pts, [1, -1, 1] * orig_pts + [0, 2.0 * Ny, 0]))
427 1
        base_pts = np.vstack((base_pts, [1, 1, -1] * orig_pts + [0, 0, 2.0 * Nz]))
428 1
        base_pts = np.vstack((base_pts, [-1, 1, 1] * orig_pts))
429 1
        base_pts = np.vstack((base_pts, [1, -1, 1] * orig_pts))
430 1
        base_pts = np.vstack((base_pts, [1, 1, -1] * orig_pts))
431 1
    vor = sptl.Voronoi(points=base_pts)
432 1
    vor.vertices = np.around(vor.vertices)
433 1
    vor.vertices *= (np.array(im.shape) - 1) / np.array(im.shape)
434 1
    vor.edges = _get_Voronoi_edges(vor)
435 1
    for row in vor.edges:
436 1
        pts = vor.vertices[row].astype(int)
437 1
        if np.all(pts >= 0) and np.all(pts < im.shape):
438 1
            line_pts = line_segment(pts[0], pts[1])
439 1
            im[tuple(line_pts)] = True
440 1
    im = edt(~im) > radius
441 1
    return im
442

443

444 1
def _get_Voronoi_edges(vor):
445
    r"""
446
    Given a Voronoi object as produced by the scipy.spatial.Voronoi class,
447
    this function calculates the start and end points of eeach edge in the
448
    Voronoi diagram, in terms of the vertex indices used by the received
449
    Voronoi object.
450

451
    Parameters
452
    ----------
453
    vor : scipy.spatial.Voronoi object
454

455
    Returns
456
    -------
457
    A 2-by-N array of vertex indices, indicating the start and end points of
458
    each vertex in the Voronoi diagram.  These vertex indices can be used to
459
    index straight into the ``vor.vertices`` array to get spatial positions.
460
    """
461 1
    edges = [[], []]
462 1
    for facet in vor.ridge_vertices:
463
        # Create a closed cycle of vertices that define the facet
464 1
        edges[0].extend(facet[:-1] + [facet[-1]])
465 1
        edges[1].extend(facet[1:] + [facet[0]])
466 1
    edges = np.vstack(edges).T  # Convert to scipy-friendly format
467 1
    mask = np.any(edges == -1, axis=1)  # Identify edges at infinity
468 1
    edges = edges[~mask]  # Remove edges at infinity
469 1
    edges = np.sort(edges, axis=1)  # Move all points to upper triangle
470
    # Remove duplicate pairs
471 1
    edges = edges[:, 0] + 1j * edges[:, 1]  # Convert to imaginary
472 1
    edges = np.unique(edges)  # Remove duplicates
473 1
    edges = np.vstack((np.real(edges), np.imag(edges))).T  # Back to real
474 1
    edges = np.array(edges, dtype=int)
475 1
    return edges
476

477

478 1
def lattice_spheres(
479
    shape: List[int], radius: int, offset: int = 0, lattice: str = "sc"
480
):
481
    r"""
482
    Generates a cubic packing of spheres in a specified lattice arrangement
483

484
    Parameters
485
    ----------
486
    shape : list
487
        The size of the image to generate in [Nx, Ny, Nz] where N is the
488
        number of voxels in each direction.  For a 2D image, use [Nx, Ny].
489

490
    radius : scalar
491
        The radius of spheres (circles) in the packing
492

493
    offset : scalar
494
        The amount offset (+ or -) to add between sphere centers.
495

496
    lattice : string
497
        Specifies the type of lattice to create.  Options are:
498

499
        'sc' - Simple Cubic (default)
500

501
        'fcc' - Face Centered Cubic
502

503
        'bcc' - Body Centered Cubic
504

505
        For 2D images, 'sc' gives a square lattice and both 'fcc' and 'bcc'
506
        give a triangular lattice.
507

508
    Returns
509
    -------
510
    image : ND-array
511
        A boolean array with ``True`` values denoting the pore space
512
    """
513 1
    print(60 * "-")
514 1
    print("lattice_spheres: Generating " + lattice + " lattice")
515 1
    r = radius
516 1
    shape = np.array(shape)
517 1
    if np.size(shape) == 1:
518 0
        shape = np.full((3,), int(shape))
519 1
    im = np.zeros(shape, dtype=bool)
520 1
    im = im.squeeze()
521

522
    # Parse lattice type
523 1
    lattice = lattice.lower()
524 1
    if im.ndim == 2:
525 1
        if lattice in ["sc"]:
526 1
            lattice = "sq"
527 1
        if lattice in ["bcc", "fcc"]:
528 0
            lattice = "tri"
529

530 1
    if lattice in ["sq", "square"]:
531 1
        spacing = 2 * r
532 1
        s = int(spacing / 2) + np.array(offset)
533 1
        coords = np.mgrid[r : im.shape[0] - r : 2 * s, r : im.shape[1] - r : 2 * s]
534 1
        im[coords[0], coords[1]] = 1
535 1
    elif lattice in ["tri", "triangular"]:
536 1
        spacing = 2 * np.floor(np.sqrt(2 * (r ** 2))).astype(int)
537 1
        s = int(spacing / 2) + offset
538 1
        coords = np.mgrid[r : im.shape[0] - r : 2 * s, r : im.shape[1] - r : 2 * s]
539 1
        im[coords[0], coords[1]] = 1
540 1
        coords = np.mgrid[
541
            s + r : im.shape[0] - r : 2 * s, s + r : im.shape[1] - r : 2 * s
542
        ]
543 1
        im[coords[0], coords[1]] = 1
544 1
    elif lattice in ["sc", "simple cubic", "cubic"]:
545 1
        spacing = 2 * r
546 1
        s = int(spacing / 2) + np.array(offset)
547 1
        coords = np.mgrid[
548
            r : im.shape[0] - r : 2 * s,
549
            r : im.shape[1] - r : 2 * s,
550
            r : im.shape[2] - r : 2 * s,
551
        ]
552 1
        im[coords[0], coords[1], coords[2]] = 1
553 1
    elif lattice in ["bcc", "body cenetered cubic"]:
554 1
        spacing = 2 * np.floor(np.sqrt(4 / 3 * (r ** 2))).astype(int)
555 1
        s = int(spacing / 2) + offset
556 1
        coords = np.mgrid[
557
            r : im.shape[0] - r : 2 * s,
558
            r : im.shape[1] - r : 2 * s,
559
            r : im.shape[2] - r : 2 * s,
560
        ]
561 1
        im[coords[0], coords[1], coords[2]] = 1
562 1
        coords = np.mgrid[
563
            s + r : im.shape[0] - r : 2 * s,
564
            s + r : im.shape[1] - r : 2 * s,
565
            s + r : im.shape[2] - r : 2 * s,
566
        ]
567 1
        im[coords[0], coords[1], coords[2]] = 1
568 1
    elif lattice in ["fcc", "face centered cubic"]:
569 1
        spacing = 2 * np.floor(np.sqrt(2 * (r ** 2))).astype(int)
570 1
        s = int(spacing / 2) + offset
571 1
        coords = np.mgrid[
572
            r : im.shape[0] - r : 2 * s,
573
            r : im.shape[1] - r : 2 * s,
574
            r : im.shape[2] - r : 2 * s,
575
        ]
576 1
        im[coords[0], coords[1], coords[2]] = 1
577 1
        coords = np.mgrid[
578
            r : im.shape[0] - r : 2 * s,
579
            s + r : im.shape[1] - r : 2 * s,
580
            s + r : im.shape[2] - r : 2 * s,
581
        ]
582 1
        im[coords[0], coords[1], coords[2]] = 1
583 1
        coords = np.mgrid[
584
            s + r : im.shape[0] - r : 2 * s,
585
            s : im.shape[1] - r : 2 * s,
586
            s + r : im.shape[2] - r : 2 * s,
587
        ]
588 1
        im[coords[0], coords[1], coords[2]] = 1
589 1
        coords = np.mgrid[
590
            s + r : im.shape[0] - r : 2 * s,
591
            s + r : im.shape[1] - r : 2 * s,
592
            s : im.shape[2] - r : 2 * s,
593
        ]
594 1
        im[coords[0], coords[1], coords[2]] = 1
595 1
    im = ~(edt(~im) < r)
596 1
    return im
597

598

599 1
def overlapping_spheres(shape: List[int],
600
                        radius: int,
601
                        porosity: float,
602
                        iter_max: int = 10,
603
                        tol: float = 0.01):
604
    r"""
605
    Generate a packing of overlapping mono-disperse spheres
606

607
    Parameters
608
    ----------
609
    shape : list
610
        The size of the image to generate in [Nx, Ny, Nz] where Ni is the
611
        number of voxels in the i-th direction.
612

613
    radius : scalar
614
        The radius of spheres in the packing.
615

616
    porosity : scalar
617
        The porosity of the final image, accurate to the given tolerance.
618

619
    iter_max : int
620
        Maximum number of iterations for the iterative algorithm that improves
621
        the porosity of the final image to match the given value.
622

623
    tol : float
624
        Tolerance for porosity of the final image compared to the given value.
625

626
    Returns
627
    -------
628
    image : ND-array
629
        A boolean array with ``True`` values denoting the pore space
630

631
    Notes
632
    -----
633
    This method can also be used to generate a dispersion of hollows by
634
    treating ``porosity`` as solid volume fraction and inverting the
635
    returned image.
636

637
    """
638 1
    shape = np.array(shape)
639 1
    if np.size(shape) == 1:
640 0
        shape = np.full((3, ), int(shape))
641 1
    ndim = (shape != 1).sum()
642 1
    s_vol = ps_disk(radius).sum() if ndim == 2 else ps_ball(radius).sum()
643

644 1
    bulk_vol = np.prod(shape)
645 1
    N = int(np.ceil((1 - porosity) * bulk_vol / s_vol))
646 1
    im = np.random.random(size=shape)
647

648
    # Helper functions for calculating porosity: phi = g(f(N))
649 1
    def f(N):
650 1
        return edt(im > N / bulk_vol) < radius
651

652 1
    def g(im):
653
        r"""Returns fraction of 0s, given a binary image"""
654 1
        return 1 - im.sum() / np.prod(shape)
655

656
    # # Newton's method for getting image porosity match the given
657
    # w = 1.0                         # Damping factor
658
    # dN = 5 if ndim == 2 else 25     # Perturbation
659
    # for i in range(iter_max):
660
    #     err = g(f(N)) - porosity
661
    #     d_err = (g(f(N+dN)) - g(f(N))) / dN
662
    #     if d_err == 0:
663
    #         break
664
    #     if abs(err) <= tol:
665
    #         break
666
    #     N2 = N - int(err/d_err)   # xnew = xold - f/df
667
    #     N = w * N2 + (1-w) * N
668

669
    # Bisection search: N is always undershoot (bc. of overlaps)
670 1
    N_low, N_high = N, 4 * N
671 1
    for i in range(iter_max):
672 1
        N = np.mean([N_high, N_low], dtype=int)
673 1
        err = g(f(N)) - porosity
674 1
        if err > 0:
675 1
            N_low = N
676
        else:
677 1
            N_high = N
678 1
        if abs(err) <= tol:
679 1
            break
680

681 1
    return ~f(N)
682

683

684 1
def perlin_noise(shape: List[int], porosity=None, octaves: int = 3,
685
                 frequency: List[int] = 2, persistence: float = 0.5):
686
    r"""
687
    Generate a Perlin noise field
688

689
    Parameters
690
    ----------
691
    shape : array_like
692
        The shape of the desired image
693
    frequncy : array_like
694
        Controls the frequency of the noise, with higher values leading to
695
        smaller features or more tightly spaced undulations in the brightness.
696
    porosity : float
697
        If specified, the returned image will be thresholded to the specified
698
        porosity.  If not provided, the greyscale noise is returned (default).
699
    octaves : int
700
        Controls the texture of the noise, with higher values giving more
701
        comlex features of larger length scales.
702
    persistence : float
703
        Controls how prominent each successive octave is.  Shoul be a number
704
        less than 1.
705

706
    Returns
707
    -------
708
    An ND-array of the specified ``shape``.  If ``porosity`` is not given
709
    then the array contains greyscale values distributed normally about 0.
710
    Use ``porespy.tools.norm_to_uniform`` to create an well-scale image for
711
    thresholding.  If ``porosity`` is given then these steps are done
712
    internally and a boolean image is returned.
713

714
    Notes
715
    -----
716
    The implementation used here is a bit fussy about the values of
717
    ``frequency`` and ``octaves``.  (1) the image ``shape`` must an integer
718
    multiple of ``frequency`` in each direction, and (2) ``frequency`` to the
719
    power of ``octaves`` must be less than or equal the``shape`` in each
720
    direction.  Exceptions are thrown if these conditions are not met.
721

722
    References
723
    ----------
724
    This implementation is taken from Pierre Vigier's
725
    `Github repo <https://github.com/pvigier/perlin-numpy>`_
726

727
    """
728
    # Parse args
729 1
    shape = np.array(shape)
730 1
    if shape.size == 1:  # Assume 3D
731 0
        shape = np.ones(3, dtype=int) * shape
732 1
    res = np.array(frequency)
733 1
    if res.size == 1:  # Assume shape as shape
734 1
        res = np.ones(shape.size, dtype=int) * res
735

736
    # Check inputs for various sins
737 1
    if res.size != shape.size:
738 0
        raise Exception('shape and res must have same dimensions')
739 1
    if np.any(np.mod(shape, res) > 0):
740 0
        raise Exception('res must be a multiple of shape along each axis')
741 1
    if np.any(shape / res**octaves < 1):
742 0
        raise Exception('(res[i])**octaves must be <= shape[i]')
743 1
    check = shape / (res**octaves)
744 1
    if np.any(check % 1):
745 0
        raise Exception("Image size must be factor of res**octaves")
746

747
    # Generate noise
748 1
    noise = np.zeros(shape)
749 1
    frequency = 1
750 1
    amplitude = 1
751 1
    for _ in tqdm(range(octaves), file=sys.stdout):
752 1
        if noise.ndim == 2:
753 1
            noise += amplitude * _perlin_noise_2D(shape, frequency * res)
754 1
        elif noise.ndim == 3:
755 1
            noise += amplitude * _perlin_noise_3D(shape, frequency * res)
756 1
        frequency *= 2
757 1
        amplitude *= persistence
758

759 1
    if porosity is not None:
760 0
        noise = norm_to_uniform(noise, scale=[0, 1])
761 0
        noise = noise > porosity
762

763 1
    return noise
764

765

766 1
def _perlin_noise_3D(shape, res):
767 1
    def f(t):
768 1
        return 6 * t**5 - 15 * t**4 + 10 * t**3
769

770 1
    delta = res / shape
771 1
    d = shape // res
772 1
    grid = np.mgrid[0:res[0]:delta[0], 0:res[1]:delta[1], 0:res[2]:delta[2]]
773 1
    grid = grid.transpose(1, 2, 3, 0) % 1
774
    # Gradients
775 1
    theta = 2 * np.pi * np.random.rand(*(res + 1))
776 1
    phi = 2 * np.pi * np.random.rand(*(res + 1))
777 1
    gradients = np.stack((np.sin(phi) * np.cos(theta),
778
                          np.sin(phi) * np.sin(theta),
779
                          np.cos(phi)), axis=3)
780 1
    g000 = gradients[0:-1, 0:-1, 0:-1]
781 1
    g000 = g000.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
782 1
    g100 = gradients[1:, 0:-1, 0:-1]
783 1
    g100 = g100.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
784 1
    g010 = gradients[0:-1, 1:, 0:-1]
785 1
    g010 = g010.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
786 1
    g110 = gradients[1:, 1:, 0:-1]
787 1
    g110 = g110.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
788 1
    g001 = gradients[0:-1, 0:-1, 1:]
789 1
    g001 = g001.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
790 1
    g101 = gradients[1:, 0:-1, 1:]
791 1
    g101 = g101.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
792 1
    g011 = gradients[0:-1, 1:, 1:]
793 1
    g011 = g011.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
794 1
    g111 = gradients[1:, 1:, 1:]
795 1
    g111 = g111.repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
796
    # Ramps
797 1
    n000 = np.sum(np.stack((grid[..., 0],
798
                            grid[..., 1],
799
                            grid[..., 2]), axis=3) * g000, 3)
800 1
    n100 = np.sum(np.stack((grid[..., 0] - 1,
801
                            grid[..., 1],
802
                            grid[..., 2]), axis=3) * g100, 3)
803 1
    n010 = np.sum(np.stack((grid[..., 0],
804
                            grid[..., 1] - 1,
805
                            grid[..., 2]), axis=3) * g010, 3)
806 1
    n110 = np.sum(np.stack((grid[..., 0] - 1,
807
                            grid[..., 1] - 1,
808
                            grid[..., 2]), axis=3) * g110, 3)
809 1
    n001 = np.sum(np.stack((grid[..., 0],
810
                            grid[..., 1],
811
                            grid[..., 2] - 1), axis=3) * g001, 3)
812 1
    n101 = np.sum(np.stack((grid[..., 0] - 1,
813
                            grid[..., 1],
814
                            grid[..., 2] - 1), axis=3) * g101, 3)
815 1
    n011 = np.sum(np.stack((grid[..., 0],
816
                            grid[..., 1] - 1,
817
                            grid[..., 2] - 1), axis=3) * g011, 3)
818 1
    n111 = np.sum(np.stack((grid[..., 0] - 1,
819
                            grid[..., 1] - 1,
820
                            grid[..., 2] - 1), axis=3) * g111, 3)
821
    # Interpolation
822 1
    t = f(grid)
823 1
    n00 = n000 * (1 - t[..., 0]) + t[..., 0] * n100
824 1
    n10 = n010 * (1 - t[..., 0]) + t[..., 0] * n110
825 1
    n01 = n001 * (1 - t[..., 0]) + t[..., 0] * n101
826 1
    n11 = n011 * (1 - t[..., 0]) + t[..., 0] * n111
827 1
    n0 = (1 - t[..., 1]) * n00 + t[..., 1] * n10
828 1
    n1 = (1 - t[..., 1]) * n01 + t[..., 1] * n11
829 1
    return ((1 - t[..., 2]) * n0 + t[..., 2] * n1)
830

831

832 1
def _perlin_noise_2D(shape, res):
833 1
    def f(t):
834 1
        return 6 * t**5 - 15 * t**4 + 10 * t**3
835

836 1
    delta = res / shape
837 1
    d = shape // res
838 1
    grid = np.mgrid[0:res[0]:delta[0],
839
                    0:res[1]:delta[1]].transpose(1, 2, 0) % 1
840

841
    # Gradients
842 1
    angles = 2 * np.pi * np.random.rand(res[0] + 1, res[1] + 1)
843 1
    gradients = np.dstack((np.cos(angles), np.sin(angles)))
844 1
    g00 = gradients[0:-1, 0:-1].repeat(d[0], 0).repeat(d[1], 1)
845 1
    g10 = gradients[1:, 0:-1].repeat(d[0], 0).repeat(d[1], 1)
846 1
    g01 = gradients[0:-1, 1:].repeat(d[0], 0).repeat(d[1], 1)
847 1
    g11 = gradients[1:, 1:].repeat(d[0], 0).repeat(d[1], 1)
848

849
    # Ramps
850 1
    n00 = np.sum(np.dstack((grid[..., 0], grid[..., 1])) * g00, 2)
851 1
    n10 = np.sum(np.dstack((grid[..., 0] - 1, grid[..., 1])) * g10, 2)
852 1
    n01 = np.sum(np.dstack((grid[..., 0], grid[..., 1] - 1)) * g01, 2)
853 1
    n11 = np.sum(np.dstack((grid[..., 0] - 1, grid[..., 1] - 1)) * g11, 2)
854

855
    # Interpolation
856 1
    t = f(grid)
857 1
    n0 = n00 * (1 - t[:, :, 0]) + t[:, :, 0] * n10
858 1
    n1 = n01 * (1 - t[:, :, 0]) + t[:, :, 0] * n11
859

860 1
    return np.sqrt(2) * ((1 - t[:, :, 1]) * n0 + t[:, :, 1] * n1)
861

862

863 1
def blobs(shape: List[int], porosity: float = 0.5, blobiness: int = 1,
864
          **kwargs):
865
    """
866
    Generates an image containing amorphous blobs
867

868
    Parameters
869
    ----------
870
    shape : list
871
        The size of the image to generate in [Nx, Ny, Nz] where N is the
872
        number of voxels
873

874
    porosity : float
875
        If specified, this will threshold the image to the specified value
876
        prior to returning.  If ``None`` is specified, then the scalar noise
877
        field is converted to a uniform distribution and returned without
878
        thresholding.
879

880
    blobiness : int or list of ints(default = 1)
881
        Controls the morphology of the blobs.  A higher number results in
882
        a larger number of small blobs.  If a list is supplied then the blobs
883
        are anisotropic.
884

885
    Returns
886
    -------
887
    image : ND-array
888
        A boolean array with ``True`` values denoting the pore space
889

890
    See Also
891
    --------
892
    norm_to_uniform
893

894
    """
895 1
    blobiness = np.array(blobiness)
896 1
    shape = np.array(shape)
897 1
    parallel = kwargs.pop('parallel', False)
898 1
    divs = kwargs.pop('divs', 2)
899 1
    cores = kwargs.pop('cores', None)
900 1
    if np.size(shape) == 1:
901 1
        shape = np.full((3, ), int(shape))
902 1
    sigma = np.mean(shape) / (40 * blobiness)
903 1
    im = np.random.random(shape)
904 1
    if parallel:
905
        # TODO: The determination of the overlap should be done rigorously
906 1
        im = ps.filters.chunked_func(func=spim.gaussian_filter,
907
                                     input=im, sigma=sigma,
908
                                     divs=divs, cores=cores, overlap=10)
909
    else:
910 1
        im = spim.gaussian_filter(im, sigma=sigma)
911 1
    im = norm_to_uniform(im, scale=[0, 1])
912 1
    if porosity:
913 1
        im = im < porosity
914 1
    return im
915

916

917 1
def _cylinders(shape: List[int],
918
               radius: int,
919
               ncylinders: int,
920
               phi_max: float = 0,
921
               theta_max: float = 90,
922
               length: float = None,
923
               verbose: bool = True):
924
    r"""
925
    Generates a binary image of overlapping cylinders.
926

927
    This is a good approximation of a fibrous mat.
928

929
    Parameters
930
    ----------
931
    shape : list
932
        The size of the image to generate in [Nx, Ny, Nz] where N is the
933
        number of voxels. 2D images are not permitted.
934
    radius : scalar
935
        The radius of the cylinders in voxels
936
    ncylinders : scalar
937
        The number of cylinders to add to the domain. Adjust this value to
938
        control the final porosity, which is not easily specified since
939
        cylinders overlap and intersect different fractions of the domain.
940
    phi_max : scalar
941
        A value between 0 and 90 that controls the amount that the cylinders
942
        lie *out of* the XY plane, with 0 meaning all cylinders lie in the XY
943
        plane, and 90 meaning that cylinders are randomly oriented out of the
944
        plane by as much as +/- 90 degrees.
945
    theta_max : scalar
946
        A value between 0 and 90 that controls the amount of rotation *in the*
947
        XY plane, with 0 meaning all cylinders point in the X-direction, and
948
        90 meaning they are randomly rotated about the Z axis by as much
949
        as +/- 90 degrees.
950
    length : scalar
951
        The length of the cylinders to add.  If ``None`` (default) then the
952
        cylinders will extend beyond the domain in both directions so no ends
953
        will exist. If a scalar value is given it will be interpreted as the
954
        Euclidean distance between the two ends of the cylinder.  Note that
955
        one or both of the ends *may* still lie outside the domain, depending
956
        on the randomly chosen center point of the cylinder.
957

958
    Returns
959
    -------
960
    image : ND-array
961
        A boolean array with ``True`` values denoting the pore space
962
    """
963 1
    shape = np.array(shape)
964 1
    if np.size(shape) == 1:
965 0
        shape = np.full((3, ), int(shape))
966 1
    elif np.size(shape) == 2:
967 1
        raise Exception("2D cylinders don't make sense")
968
    # Find hypotenuse of domain from [0,0,0] to [Nx,Ny,Nz]
969 1
    H = np.sqrt(np.sum(np.square(shape))).astype(int)
970 1
    if length is None:  # Assume cylinders span domain if length not given
971 1
        length = 2 * H
972 1
    R = min(int(length / 2), 2 * H)  # Trim given length to 2H if too long
973
    # Adjust max angles to be between 0 and 90
974 1
    if (phi_max > 90) or (phi_max < 0):
975 0
        raise Exception('phi_max must be betwen 0 and 90')
976 1
    if (theta_max > 90) or (theta_max < 0):
977 0
        raise Exception('theta_max must be betwen 0 and 90')
978
    # Create empty image for inserting into
979 1
    im = np.zeros(shape, dtype=bool)
980 1
    n = 0
981 1
    L = min(H, R)
982 1
    pbar = tqdm(total=ncylinders, file=sys.stdout, disable=not verbose)
983 1
    while n < ncylinders:
984
        # Choose a random starting point in domain
985 1
        x = np.random.rand(3) * (shape + 2 * L)
986
        # Chose a random phi and theta within given ranges
987 1
        phi = (np.pi / 2 - np.pi * np.random.rand()) * phi_max / 90
988 1
        theta = (np.pi / 2 - np.pi * np.random.rand()) * theta_max / 90
989 1
        X0 = R * np.array([np.cos(phi) * np.cos(theta),
990
                           np.cos(phi) * np.sin(theta),
991
                           np.sin(phi)])
992 1
        [X0, X1] = [x + X0, x - X0]
993 1
        crds = line_segment(X0, X1)
994 1
        lower = ~np.any(np.vstack(crds).T < [L, L, L], axis=1)
995 1
        upper = ~np.any(np.vstack(crds).T >= shape + L, axis=1)
996 1
        valid = upper * lower
997 1
        if np.any(valid):
998 1
            im[crds[0][valid] - L, crds[1][valid] - L, crds[2][valid] - L] = 1
999 1
            n += 1
1000 1
            pbar.update()
1001 1
    im = np.array(im, dtype=bool)
1002 1
    dt = edt(~im) < radius
1003 1
    return ~dt
1004

1005

1006 1
def cylinders(shape: List[int],
1007
              radius: int,
1008
              ncylinders: int = None,
1009
              porosity: float = None,
1010
              phi_max: float = 0,
1011
              theta_max: float = 90,
1012
              length: float = None,
1013
              max_iter: int = 3):
1014
    r"""
1015
    Generates a binary image of overlapping cylinders given porosity OR number
1016
    of cylinders.
1017

1018
    This is a good approximation of a fibrous mat.
1019

1020
    Parameters
1021
    ----------
1022
    shape : list
1023
        The size of the image to generate in [Nx, Ny, Nz] where N is the
1024
        number of voxels. 2D images are not permitted.
1025
    radius : scalar
1026
        The radius of the cylinders in voxels
1027
    ncylinders : scalar
1028
        The number of cylinders to add to the domain. Adjust this value to
1029
        control the final porosity, which is not easily specified since
1030
        cylinders overlap and intersect different fractions of the domain.
1031
    porosity : scalar
1032
        The targeted value for the porosity of the generated mat. The
1033
        function uses an algorithm for predicted the number of required
1034
        number of cylinder, and refines this over a certain number of
1035
        fractional insertions (according to the 'iterations' input).
1036
    phi_max : scalar
1037
        A value between 0 and 90 that controls the amount that the cylinders
1038
        lie *out of* the XY plane, with 0 meaning all cylinders lie in the XY
1039
        plane, and 90 meaning that cylinders are randomly oriented out of the
1040
        plane by as much as +/- 90 degrees.
1041
    theta_max : scalar
1042
        A value between 0 and 90 that controls the amount of rotation *in the*
1043
        XY plane, with 0 meaning all cylinders point in the X-direction, and
1044
        90 meaning they are randomly rotated about the Z axis by as much
1045
        as +/- 90 degrees.
1046
    length : scalar
1047
        The length of the cylinders to add.  If ``None`` (default) then the
1048
        cylinders will extend beyond the domain in both directions so no ends
1049
        will exist. If a scalar value is given it will be interpreted as the
1050
        Euclidean distance between the two ends of the cylinder.  Note that
1051
        one or both of the ends *may* still lie outside the domain, depending
1052
        on the randomly chosen center point of the cylinder.
1053
    max_iter : scalar
1054
        The number of fractional fiber insertions used to target the requested
1055
        porosity. By default a value of 3 is used (and this is typically
1056
        effective in getting very close to the targeted porosity), but a
1057
        greater number can be input to improve the achieved porosity.
1058
    return_fiber_number : bool
1059
        Determines whether the function will return the number of fibers
1060
        along with the image
1061

1062
    Returns
1063
    -------
1064
    image : ND-array
1065
        A boolean array with ``True`` values denoting the pore space
1066

1067
    Notes
1068
    -----
1069
    The cylinders_porosity function works by estimating the number of
1070
    cylinders needed to be inserted into the domain by estimating
1071
    cylinder length, and exploiting the fact that, when inserting any
1072
    potentially overlapping objects randomly into a volume v_total (which
1073
    has units of pixels and is equal to dimx x dimy x dimz, for example),
1074
    such that the total volume of objects added to the volume is v_added
1075
    (and includes any volume that was inserted but overlapped with already
1076
    occupied space), the resulting porosity will be equal to
1077
    exp(-v_added/v_total).
1078

1079
    After intially estimating the cylinder number and inserting a small
1080
    fraction of the estimated number, the true cylinder volume is
1081
    calculated, the estimate refined, and a larger fraction of cylinders
1082
    inserted. This is repeated a number of times according to the
1083
    'max_iter' argument, yielding an image with a porosity close to
1084
    the goal.
1085

1086
    """
1087 1
    if ncylinders is not None:
1088 1
        im = _cylinders(
1089
            shape=shape,
1090
            radius=radius,
1091
            ncylinders=ncylinders,
1092
            phi_max=phi_max,
1093
            theta_max=theta_max,
1094
            length=length,
1095
        )
1096 1
        return im
1097

1098 1
    if porosity is None:
1099 0
        raise Exception("'ncylinders' and 'porosity' can't be both None")
1100

1101 1
    if max_iter < 3:
1102 1
        raise Exception("Iterations must be greater than or equal to 3")
1103

1104 1
    vol_total = float(np.prod(shape))
1105

1106 1
    def get_num_pixels(porosity):
1107
        r"""
1108
        Helper method to calculate number of pixels given a porosity
1109
        """
1110 1
        return -np.log(porosity) * vol_total
1111

1112
    # Crudely estimate fiber length as cube root of product of dims
1113 1
    length_estimate = vol_total ** (1 / 3) if length is None else length
1114

1115
    # Rough fiber volume estimate
1116 1
    vol_fiber = length_estimate * np.pi * radius * radius
1117 1
    n_pixels_to_add = get_num_pixels(porosity)
1118

1119
    # Rough estimate of n_fibers
1120 1
    n_fibers_added = 0
1121
    # Calculate fraction of fibers to be added in each iteration.
1122 1
    subdif = 0.8 / np.sum(np.arange(1, max_iter) ** 2)
1123 1
    fractions = [0.2]
1124 1
    for i in range(1, max_iter):
1125 1
        fractions.append(fractions[i - 1] + (max_iter - i) ** 2 * subdif)
1126

1127 1
    im = np.ones(shape, dtype=bool)
1128 1
    for frac in tqdm(fractions, file=sys.stdout, desc="Adding fibers"):
1129 1
        n_fibers_total = n_pixels_to_add / vol_fiber
1130 1
        n_fibers = int(np.ceil(frac * n_fibers_total) - n_fibers_added)
1131 1
        if n_fibers > 0:
1132 1
            im = im & _cylinders(
1133
                shape, radius, n_fibers, phi_max, theta_max, length, verbose=False
1134
            )
1135 1
        n_fibers_added += n_fibers
1136
        # Update parameters for next iteration
1137 1
        porosity = ps.metrics.porosity(im)
1138 1
        vol_added = get_num_pixels(porosity)
1139 1
        vol_fiber = vol_added / n_fibers_added
1140

1141 1
    print(f"{n_fibers_added} fibers were added to reach the target porosity.\n")
1142

1143 1
    return im
1144

1145

1146 1
def line_segment(X0, X1):
1147
    r"""
1148
    Calculate the voxel coordinates of a straight line between the two given
1149
    end points
1150

1151
    Parameters
1152
    ----------
1153
    X0 and X1 : array_like
1154
        The [x, y] or [x, y, z] coordinates of the start and end points of
1155
        the line.
1156

1157
    Returns
1158
    -------
1159
    coords : list of lists
1160
        A list of lists containing the X, Y, and Z coordinates of all voxels
1161
        that should be drawn between the start and end points to create a solid
1162
        line.
1163
    """
1164 1
    X0 = np.around(X0).astype(int)
1165 1
    X1 = np.around(X1).astype(int)
1166 1
    if len(X0) == 3:
1167 1
        L = np.amax(np.absolute([[X1[0] - X0[0]], [X1[1] - X0[1]], [X1[2] - X0[2]]])) + 1
1168 1
        x = np.rint(np.linspace(X0[0], X1[0], L)).astype(int)
1169 1
        y = np.rint(np.linspace(X0[1], X1[1], L)).astype(int)
1170 1
        z = np.rint(np.linspace(X0[2], X1[2], L)).astype(int)
1171 1
        return [x, y, z]
1172
    else:
1173 1
        L = np.amax(np.absolute([[X1[0] - X0[0]], [X1[1] - X0[1]]])) + 1
1174 1
        x = np.rint(np.linspace(X0[0], X1[0], L)).astype(int)
1175 1
        y = np.rint(np.linspace(X0[1], X1[1], L)).astype(int)
1176 1
        return [x, y]

Read our documentation on viewing source code .

Loading