1 1
import sys
2 1
import numpy as np
3 1
import openpnm as op
4 1
from tqdm import tqdm
5 1
import scipy.ndimage as spim
6 1
from porespy.tools import extend_slice
7 1
import openpnm.models.geometry as op_gm
8

9

10 1
def regions_to_network(im, dt=None, voxel_size=1):
11
    r"""
12
    Analyzes an image that has been partitioned into pore regions and extracts
13
    the pore and throat geometry as well as network connectivity.
14

15
    Parameters
16
    ----------
17
    im : ND-array
18
        An image of the pore space partitioned into individual pore regions.
19
        Note that this image must have zeros indicating the solid phase.
20

21
    dt : ND-array
22
        The distance transform of the pore space.  If not given it will be
23
        calculated, but it can save time to provide one if available.
24

25
    voxel_size : scalar
26
        The resolution of the image, expressed as the length of one side of a
27
        voxel, so the volume of a voxel would be **voxel_size**-cubed.  The
28
        default is 1, which is useful when overlaying the PNM on the original
29
        image since the scale of the image is alway 1 unit lenth per voxel.
30

31
    Returns
32
    -------
33
    A dictionary containing all the pore and throat size data, as well as the
34
    network topological information.  The dictionary names use the OpenPNM
35
    convention (i.e. 'pore.coords', 'throat.conns') so it may be converted
36
    directly to an OpenPNM network object using the ``update`` command.
37

38
    """
39 1
    print('-'*60)
40 1
    print('Extracting pore and throat information from image')
41 1
    from skimage.morphology import disk, ball
42 1
    struc_elem = disk if im.ndim == 2 else ball
43

44
    # if ~np.any(im == 0):
45
    #     raise Exception('The received image has no solid phase (0\'s)')
46

47 1
    if dt is None:
48 1
        dt = spim.distance_transform_edt(im > 0)
49 1
        dt = spim.gaussian_filter(input=dt, sigma=0.5)
50

51
    # Get 'slices' into im for each pore region
52 1
    slices = spim.find_objects(im)
53

54
    # Initialize arrays
55 1
    Ps = np.arange(1, np.amax(im)+1)
56 1
    Np = np.size(Ps)
57 1
    p_coords = np.zeros((Np, im.ndim), dtype=float)
58 1
    p_volume = np.zeros((Np, ), dtype=float)
59 1
    p_dia_local = np.zeros((Np, ), dtype=float)
60 1
    p_dia_global = np.zeros((Np, ), dtype=float)
61 1
    p_label = np.zeros((Np, ), dtype=int)
62 1
    p_area_surf = np.zeros((Np, ), dtype=int)
63 1
    t_conns = []
64 1
    t_dia_inscribed = []
65 1
    t_area = []
66 1
    t_perimeter = []
67 1
    t_coords = []
68
    # dt_shape = np.array(dt.shape)
69

70
    # Start extracting size information for pores and throats
71 1
    for i in tqdm(Ps, file=sys.stdout):
72 1
        pore = i - 1
73 1
        if slices[pore] is None:
74 0
            continue
75 1
        s = extend_slice(slices[pore], im.shape)
76 1
        sub_im = im[s]
77 1
        sub_dt = dt[s]
78 1
        pore_im = sub_im == i
79 1
        padded_mask = np.pad(pore_im, pad_width=1, mode='constant')
80 1
        pore_dt = spim.distance_transform_edt(padded_mask)
81 1
        s_offset = np.array([i.start for i in s])
82 1
        p_label[pore] = i
83 1
        p_coords[pore, :] = spim.center_of_mass(pore_im) + s_offset
84 1
        p_volume[pore] = np.sum(pore_im)
85 1
        p_dia_local[pore] = (2*np.amax(pore_dt)) - np.sqrt(3)
86 1
        p_dia_global[pore] = 2*np.amax(sub_dt)
87 1
        p_area_surf[pore] = np.sum(pore_dt == 1)
88 1
        im_w_throats = spim.binary_dilation(input=pore_im, structure=struc_elem(1))
89 1
        im_w_throats = im_w_throats*sub_im
90 1
        Pn = np.unique(im_w_throats)[1:] - 1
91 1
        for j in Pn:
92 1
            if j > pore:
93 1
                t_conns.append([pore, j])
94 1
                vx = np.where(im_w_throats == (j + 1))
95 1
                t_dia_inscribed.append(2*np.amax(sub_dt[vx]))
96 1
                t_perimeter.append(np.sum(sub_dt[vx] < 2))
97 1
                t_area.append(np.size(vx[0]))
98 1
                t_inds = tuple([i+j for i, j in zip(vx, s_offset)])
99 1
                temp = np.where(dt[t_inds] == np.amax(dt[t_inds]))[0][0]
100 1
                if im.ndim == 2:
101 1
                    t_coords.append(tuple((t_inds[0][temp],
102
                                           t_inds[1][temp])))
103
                else:
104 1
                    t_coords.append(tuple((t_inds[0][temp],
105
                                           t_inds[1][temp],
106
                                           t_inds[2][temp])))
107
    # Clean up values
108 1
    Nt = len(t_dia_inscribed)  # Get number of throats
109 1
    if im.ndim == 2:  # If 2D, add 0's in 3rd dimension
110 1
        p_coords = np.vstack((p_coords.T, np.zeros((Np, )))).T
111 1
        t_coords = np.vstack((np.array(t_coords).T, np.zeros((Nt, )))).T
112

113 1
    net = {}
114 1
    net['pore.all'] = np.ones((Np, ), dtype=bool)
115 1
    net['throat.all'] = np.ones((Nt, ), dtype=bool)
116 1
    net['pore.coords'] = np.copy(p_coords)*voxel_size
117 1
    net['pore.centroid'] = np.copy(p_coords)*voxel_size
118 1
    net['throat.centroid'] = np.array(t_coords)*voxel_size
119 1
    net['throat.conns'] = np.array(t_conns)
120 1
    net['pore.label'] = np.array(p_label)
121 1
    net['pore.volume'] = np.copy(p_volume)*(voxel_size**3)
122 1
    net['throat.volume'] = np.zeros((Nt, ), dtype=float)
123 1
    net['pore.diameter'] = np.copy(p_dia_local)*voxel_size
124 1
    net['pore.inscribed_diameter'] = np.copy(p_dia_local)*voxel_size
125 1
    net['pore.equivalent_diameter'] = 2*((3/4*net['pore.volume']/np.pi)**(1/3))
126 1
    net['pore.extended_diameter'] = np.copy(p_dia_global)*voxel_size
127 1
    net['pore.surface_area'] = np.copy(p_area_surf)*(voxel_size)**2
128 1
    net['throat.diameter'] = np.array(t_dia_inscribed)*voxel_size
129 1
    net['throat.inscribed_diameter'] = np.array(t_dia_inscribed)*voxel_size
130 1
    net['throat.area'] = np.array(t_area)*(voxel_size**2)
131 1
    net['throat.perimeter'] = np.array(t_perimeter)*voxel_size
132 1
    net['throat.equivalent_diameter'] = (np.array(t_area) * (voxel_size**2))**0.5
133 1
    P12 = net['throat.conns']
134 1
    PT1 = np.sqrt(np.sum(((p_coords[P12[:, 0]]-t_coords) * voxel_size)**2, axis=1))
135 1
    PT2 = np.sqrt(np.sum(((p_coords[P12[:, 1]]-t_coords) * voxel_size)**2, axis=1))
136 1
    net['throat.total_length'] = PT1 + PT2
137 1
    PT1 = PT1-p_dia_local[P12[:, 0]]/2*voxel_size
138 1
    PT2 = PT2-p_dia_local[P12[:, 1]]/2*voxel_size
139 1
    net['throat.length'] = PT1 + PT2
140 1
    dist = (p_coords[P12[:, 0]]-p_coords[P12[:, 1]])*voxel_size
141 1
    net['throat.direct_length'] = np.sqrt(np.sum(dist**2, axis=1))
142
    # Make a dummy openpnm network to get the conduit lengths
143 1
    pn = op.network.GenericNetwork()
144 1
    pn.update(net)
145 1
    pn.add_model(propname='throat.endpoints',
146
                 model=op_gm.throat_endpoints.spherical_pores,
147
                 pore_diameter='pore.inscribed_diameter',
148
                 throat_diameter='throat.inscribed_diameter')
149 1
    pn.add_model(propname='throat.conduit_lengths',
150
                 model=op_gm.throat_length.conduit_lengths)
151 1
    pn.add_model(propname='pore.area',
152
                 model=op_gm.pore_area.sphere)
153 1
    net['throat.endpoints.head'] = pn['throat.endpoints.head']
154 1
    net['throat.endpoints.tail'] = pn['throat.endpoints.tail']
155 1
    net['throat.conduit_lengths.pore1'] = pn['throat.conduit_lengths.pore1']
156 1
    net['throat.conduit_lengths.pore2'] = pn['throat.conduit_lengths.pore2']
157 1
    net['throat.conduit_lengths.throat'] = pn['throat.conduit_lengths.throat']
158 1
    net['pore.area'] = pn['pore.area']
159 1
    prj = pn.project
160 1
    prj.clear()
161 1
    wrk = op.Workspace()
162 1
    wrk.close_project(prj)
163

164 1
    return net

Read our documentation on viewing source code .

Loading