@@ -1298,66 +1298,66 @@
Loading
1298 1298
        inlets = np.pad(inlets, mode="symmetric", pad_width=pw)
1299 1299
        # sizes = np.unique(np.around(sizes, decimals=0).astype(int))[-1::-1]
1300 1300
        imresults = np.zeros(np.shape(impad))
1301 -
        with tqdm(sizes, file=sys.stdout) as pbar:
1302 -
            for r in sizes:
1303 -
                pbar.update()
1304 -
                if parallel:
1305 -
                    imtemp = chunked_func(func=spim.binary_erosion,
1306 -
                                          input=impad, structure=strel(r),
1307 -
                                          overlap=int(2*r) + 1,
1308 -
                                          cores=cores, divs=divs)
1309 -
                elif fft:
1310 -
                    imtemp = fftmorphology(impad, strel(r), mode="erosion")
1311 -
                else:
1312 -
                    imtemp = spim.binary_erosion(input=impad,
1313 -
                                                 structure=strel(r))
1314 -
                if access_limited:
1315 -
                    imtemp = trim_disconnected_blobs(imtemp, inlets)
1301 +
        pbar = tqdm(sizes, file=sys.stdout)
1302 +
        for r in sizes:
1303 +
            pbar.update()
1304 +
            if parallel:
1305 +
                imtemp = chunked_func(func=spim.binary_erosion,
1306 +
                                      input=impad, structure=strel(r),
1307 +
                                      overlap=int(2*r) + 1,
1308 +
                                      cores=cores, divs=divs)
1309 +
            elif fft:
1310 +
                imtemp = fftmorphology(impad, strel(r), mode="erosion")
1311 +
            else:
1312 +
                imtemp = spim.binary_erosion(input=impad,
1313 +
                                             structure=strel(r))
1314 +
            if access_limited:
1315 +
                imtemp = trim_disconnected_blobs(imtemp, inlets)
1316 +
            if parallel:
1317 +
                imtemp = chunked_func(func=spim.binary_dilation,
1318 +
                                      input=imtemp, structure=strel(r),
1319 +
                                      overlap=int(2*r) + 1,
1320 +
                                      cores=cores, divs=divs)
1321 +
            elif fft:
1322 +
                imtemp = fftmorphology(imtemp, strel(r), mode="dilation")
1323 +
            else:
1324 +
                imtemp = spim.binary_dilation(input=imtemp,
1325 +
                                              structure=strel(r))
1326 +
            if np.any(imtemp):
1327 +
                imresults[(imresults == 0) * imtemp] = r
1328 +
        imresults = extract_subsection(imresults, shape=im.shape)
1329 +
    elif mode == "dt":
1330 +
        imresults = np.zeros(np.shape(im))
1331 +
        pbar = tqdm(sizes, file=sys.stdout)
1332 +
        for r in sizes:
1333 +
            pbar.update()
1334 +
            imtemp = dt >= r
1335 +
            if access_limited:
1336 +
                imtemp = trim_disconnected_blobs(imtemp, inlets)
1337 +
            if np.any(imtemp):
1338 +
                imtemp = edt(~imtemp) < r
1339 +
                imresults[(imresults == 0) * imtemp] = r
1340 +
    elif mode == "hybrid":
1341 +
        imresults = np.zeros(np.shape(im))
1342 +
        pbar = tqdm(sizes, file=sys.stdout)
1343 +
        for r in sizes:
1344 +
            pbar.update()
1345 +
            imtemp = dt >= r
1346 +
            if access_limited:
1347 +
                imtemp = trim_disconnected_blobs(imtemp, inlets)
1348 +
            if np.any(imtemp):
1316 1349
                if parallel:
1317 1350
                    imtemp = chunked_func(func=spim.binary_dilation,
1318 1351
                                          input=imtemp, structure=strel(r),
1319 1352
                                          overlap=int(2*r) + 1,
1320 1353
                                          cores=cores, divs=divs)
1321 1354
                elif fft:
1322 -
                    imtemp = fftmorphology(imtemp, strel(r), mode="dilation")
1355 +
                    imtemp = fftmorphology(imtemp, strel(r),
1356 +
                                           mode="dilation")
1323 1357
                else:
1324 1358
                    imtemp = spim.binary_dilation(input=imtemp,
1325 1359
                                                  structure=strel(r))
1326 -
                if np.any(imtemp):
1327 -
                    imresults[(imresults == 0) * imtemp] = r
1328 -
        imresults = extract_subsection(imresults, shape=im.shape)
1329 -
    elif mode == "dt":
1330 -
        imresults = np.zeros(np.shape(im))
1331 -
        with tqdm(sizes, file=sys.stdout) as pbar:
1332 -
            for r in sizes:
1333 -
                pbar.update()
1334 -
                imtemp = dt >= r
1335 -
                if access_limited:
1336 -
                    imtemp = trim_disconnected_blobs(imtemp, inlets)
1337 -
                if np.any(imtemp):
1338 -
                    imtemp = edt(~imtemp) < r
1339 -
                    imresults[(imresults == 0) * imtemp] = r
1340 -
    elif mode == "hybrid":
1341 -
        imresults = np.zeros(np.shape(im))
1342 -
        with tqdm(sizes, file=sys.stdout) as pbar:
1343 -
            for r in sizes:
1344 -
                pbar.update()
1345 -
                imtemp = dt >= r
1346 -
                if access_limited:
1347 -
                    imtemp = trim_disconnected_blobs(imtemp, inlets)
1348 -
                if np.any(imtemp):
1349 -
                    if parallel:
1350 -
                        imtemp = chunked_func(func=spim.binary_dilation,
1351 -
                                              input=imtemp, structure=strel(r),
1352 -
                                              overlap=int(2*r) + 1,
1353 -
                                              cores=cores, divs=divs)
1354 -
                    elif fft:
1355 -
                        imtemp = fftmorphology(imtemp, strel(r),
1356 -
                                               mode="dilation")
1357 -
                    else:
1358 -
                        imtemp = spim.binary_dilation(input=imtemp,
1359 -
                                                      structure=strel(r))
1360 -
                    imresults[(imresults == 0) * imtemp] = r
1360 +
                imresults[(imresults == 0) * imtemp] = r
1361 1361
    else:
1362 1362
        raise Exception("Unrecognized mode " + mode)
1363 1363
    return imresults

@@ -914,8 +914,13 @@
Loading
914 914
    return im
915 915
916 916
917 -
def cylinders(shape: List[int], radius: int, ncylinders: int,
918 -
              phi_max: float = 0, theta_max: float = 90, length: float = None):
917 +
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):
919 924
    r"""
920 925
    Generates a binary image of overlapping cylinders.
921 926
@@ -974,30 +979,170 @@
Loading
974 979
    im = np.zeros(shape, dtype=bool)
975 980
    n = 0
976 981
    L = min(H, R)
977 -
    with tqdm(range(1, ncylinders), file=sys.stdout) as pbar:
978 -
        while n < ncylinders:
979 -
            # Choose a random starting point in domain
980 -
            x = np.random.rand(3) * (shape + 2 * L)
981 -
            # Chose a random phi and theta within given ranges
982 -
            phi = (np.pi / 2 - np.pi * np.random.rand()) * phi_max / 90
983 -
            theta = (np.pi / 2 - np.pi * np.random.rand()) * theta_max / 90
984 -
            X0 = R * np.array([np.cos(phi) * np.cos(theta),
985 -
                               np.cos(phi) * np.sin(theta),
986 -
                               np.sin(phi)])
987 -
            [X0, X1] = [x + X0, x - X0]
988 -
            crds = line_segment(X0, X1)
989 -
            lower = ~np.any(np.vstack(crds).T < [L, L, L], axis=1)
990 -
            upper = ~np.any(np.vstack(crds).T >= shape + L, axis=1)
991 -
            valid = upper * lower
992 -
            if np.any(valid):
993 -
                im[crds[0][valid] - L, crds[1][valid] - L, crds[2][valid] - L] = 1
994 -
                n += 1
995 -
                pbar.update()
982 +
    pbar = tqdm(total=ncylinders, file=sys.stdout, disable=not verbose)
983 +
    while n < ncylinders:
984 +
        # Choose a random starting point in domain
985 +
        x = np.random.rand(3) * (shape + 2 * L)
986 +
        # Chose a random phi and theta within given ranges
987 +
        phi = (np.pi / 2 - np.pi * np.random.rand()) * phi_max / 90
988 +
        theta = (np.pi / 2 - np.pi * np.random.rand()) * theta_max / 90
989 +
        X0 = R * np.array([np.cos(phi) * np.cos(theta),
990 +
                           np.cos(phi) * np.sin(theta),
991 +
                           np.sin(phi)])
992 +
        [X0, X1] = [x + X0, x - X0]
993 +
        crds = line_segment(X0, X1)
994 +
        lower = ~np.any(np.vstack(crds).T < [L, L, L], axis=1)
995 +
        upper = ~np.any(np.vstack(crds).T >= shape + L, axis=1)
996 +
        valid = upper * lower
997 +
        if np.any(valid):
998 +
            im[crds[0][valid] - L, crds[1][valid] - L, crds[2][valid] - L] = 1
999 +
            n += 1
1000 +
            pbar.update()
996 1001
    im = np.array(im, dtype=bool)
997 1002
    dt = edt(~im) < radius
998 1003
    return ~dt
999 1004
1000 1005
1006 +
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 +
    if ncylinders is not None:
1088 +
        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 +
        return im
1097 +
1098 +
    if porosity is None:
1099 +
        raise Exception("'ncylinders' and 'porosity' can't be both None")
1100 +
1101 +
    if max_iter < 3:
1102 +
        raise Exception("Iterations must be greater than or equal to 3")
1103 +
1104 +
    vol_total = float(np.prod(shape))
1105 +
1106 +
    def get_num_pixels(porosity):
1107 +
        r"""
1108 +
        Helper method to calculate number of pixels given a porosity
1109 +
        """
1110 +
        return -np.log(porosity) * vol_total
1111 +
1112 +
    # Crudely estimate fiber length as cube root of product of dims
1113 +
    length_estimate = vol_total ** (1 / 3) if length is None else length
1114 +
1115 +
    # Rough fiber volume estimate
1116 +
    vol_fiber = length_estimate * np.pi * radius * radius
1117 +
    n_pixels_to_add = get_num_pixels(porosity)
1118 +
1119 +
    # Rough estimate of n_fibers
1120 +
    n_fibers_added = 0
1121 +
    # Calculate fraction of fibers to be added in each iteration.
1122 +
    subdif = 0.8 / np.sum(np.arange(1, max_iter) ** 2)
1123 +
    fractions = [0.2]
1124 +
    for i in range(1, max_iter):
1125 +
        fractions.append(fractions[i - 1] + (max_iter - i) ** 2 * subdif)
1126 +
1127 +
    im = np.ones(shape, dtype=bool)
1128 +
    for frac in tqdm(fractions, file=sys.stdout, desc="Adding fibers"):
1129 +
        n_fibers_total = n_pixels_to_add / vol_fiber
1130 +
        n_fibers = int(np.ceil(frac * n_fibers_total) - n_fibers_added)
1131 +
        if n_fibers > 0:
1132 +
            im = im & _cylinders(
1133 +
                shape, radius, n_fibers, phi_max, theta_max, length, verbose=False
1134 +
            )
1135 +
        n_fibers_added += n_fibers
1136 +
        # Update parameters for next iteration
1137 +
        porosity = ps.metrics.porosity(im)
1138 +
        vol_added = get_num_pixels(porosity)
1139 +
        vol_fiber = vol_added / n_fibers_added
1140 +
1141 +
    print(f"{n_fibers_added} fibers were added to reach the target porosity.\n")
1142 +
1143 +
    return im
1144 +
1145 +
1001 1146
def line_segment(X0, X1):
1002 1147
    r"""
1003 1148
    Calculate the voxel coordinates of a straight line between the two given

@@ -181,51 +181,51 @@
Loading
181 181
    print('Calculating regionprops')
182 182
183 183
    results = regionprops(im)
184 -
    with tqdm(range(len(results)), file=sys.stdout) as pbar:
185 -
        for i in range(len(results)):
186 -
            pbar.update()
187 -
            mask = results[i].image
188 -
            mask_padded = np.pad(mask, pad_width=1, mode='constant')
189 -
            temp = spim.distance_transform_edt(mask_padded)
190 -
            dt = extract_subsection(temp, shape=mask.shape)
191 -
192 -
            # Slice indices
193 -
            results[i].slice = results[i]._slice
194 -
195 -
            # Volume of regions in voxels
196 -
            results[i].volume = results[i].area
197 -
198 -
            # Volume of bounding box, in voxels
199 -
            results[i].bbox_volume = np.prod(mask.shape)
200 -
201 -
            # Create an image of the border
202 -
            results[i].border = dt == 1
203 -
204 -
            # Create an image of the maximal inscribed sphere
205 -
            r = dt.max()
206 -
            inv_dt = spim.distance_transform_edt(dt < r)
207 -
            results[i].inscribed_sphere = inv_dt < r
208 -
209 -
            # Find surface area using marching cubes and analyze the mesh
210 -
            tmp = np.pad(np.atleast_3d(mask), pad_width=1, mode='constant')
211 -
            tmp = spim.convolve(tmp, weights=ball(1)) / 5
212 -
            verts, faces, norms, vals = marching_cubes(volume=tmp, level=0)
213 -
            results[i].surface_mesh_vertices = verts
214 -
            results[i].surface_mesh_simplices = faces
215 -
            area = mesh_surface_area(verts, faces)
216 -
            results[i].surface_area = area
217 -
218 -
            # Find sphericity
219 -
            vol = results[i].volume
220 -
            r = (3 / 4 / np.pi * vol)**(1 / 3)
221 -
            a_equiv = 4 * np.pi * r**2
222 -
            a_region = results[i].surface_area
223 -
            results[i].sphericity = a_equiv / a_region
224 -
225 -
            # Find skeleton of region
226 -
            results[i].skeleton = skeletonize_3d(mask)
227 -
228 -
            # Volume of convex image, equal to area in 2D, so just translating
229 -
            results[i].convex_volume = results[i].convex_area
184 +
    pbar = tqdm(range(len(results)), file=sys.stdout)
185 +
    for i in range(len(results)):
186 +
        pbar.update()
187 +
        mask = results[i].image
188 +
        mask_padded = np.pad(mask, pad_width=1, mode='constant')
189 +
        temp = spim.distance_transform_edt(mask_padded)
190 +
        dt = extract_subsection(temp, shape=mask.shape)
191 +
192 +
        # Slice indices
193 +
        results[i].slice = results[i]._slice
194 +
195 +
        # Volume of regions in voxels
196 +
        results[i].volume = results[i].area
197 +
198 +
        # Volume of bounding box, in voxels
199 +
        results[i].bbox_volume = np.prod(mask.shape)
200 +
201 +
        # Create an image of the border
202 +
        results[i].border = dt == 1
203 +
204 +
        # Create an image of the maximal inscribed sphere
205 +
        r = dt.max()
206 +
        inv_dt = spim.distance_transform_edt(dt < r)
207 +
        results[i].inscribed_sphere = inv_dt < r
208 +
209 +
        # Find surface area using marching cubes and analyze the mesh
210 +
        tmp = np.pad(np.atleast_3d(mask), pad_width=1, mode='constant')
211 +
        tmp = spim.convolve(tmp, weights=ball(1)) / 5
212 +
        verts, faces, norms, vals = marching_cubes(volume=tmp, level=0)
213 +
        results[i].surface_mesh_vertices = verts
214 +
        results[i].surface_mesh_simplices = faces
215 +
        area = mesh_surface_area(verts, faces)
216 +
        results[i].surface_area = area
217 +
218 +
        # Find sphericity
219 +
        vol = results[i].volume
220 +
        r = (3 / 4 / np.pi * vol)**(1 / 3)
221 +
        a_equiv = 4 * np.pi * r**2
222 +
        a_region = results[i].surface_area
223 +
        results[i].sphericity = a_equiv / a_region
224 +
225 +
        # Find skeleton of region
226 +
        results[i].skeleton = skeletonize_3d(mask)
227 +
228 +
        # Volume of convex image, equal to area in 2D, so just translating
229 +
        results[i].convex_volume = results[i].convex_area
230 230
231 231
    return results
Files Coverage
porespy 88.1%
Project Totals (16 files) 88.1%
codecov-umbrella
Build #311800943 -
unittests
1
codecov:
2
  branch: dev
3

4
coverage:
5
  precision: 1
6
  round: down
7
  range: "50...100"
8

9
  status:
10
    project:
11
      default:
12
        target: auto
13
        threshold: 0.5%
14
        branches: null
15

16
    patch:
17
      default:
18
        target: auto
19
        threshold: 0.5%
20
        branches: null
21

22
comment:
23
  layout: "header, diff, changes, sunburst, uncovered"
24
  branches: null
25
  behavior: default
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading