Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
p_dia_global = sp.zeros((Np, ), dtype=float)
p_label = sp.zeros((Np, ), dtype=int)
p_area_surf = sp.zeros((Np, ), dtype=int)
t_conns = []
t_dia_inscribed = []
t_area = []
t_perimeter = []
t_coords = []
# dt_shape = sp.array(dt.shape)
# Start extracting size information for pores and throats
for i in tqdm(Ps):
pore = i - 1
if slices[pore] is None:
continue
s = extend_slice(slices[pore], im.shape)
sub_im = im[s]
sub_dt = dt[s]
pore_im = sub_im == i
padded_mask = sp.pad(pore_im, pad_width=1, mode='constant')
pore_dt = spim.distance_transform_edt(padded_mask)
s_offset = sp.array([i.start for i in s])
p_label[pore] = i
p_coords[pore, :] = spim.center_of_mass(pore_im) + s_offset
p_volume[pore] = sp.sum(pore_im)
p_dia_local[pore] = (2*sp.amax(pore_dt)) - sp.sqrt(3)
p_dia_global[pore] = 2*sp.amax(sub_dt)
p_area_surf[pore] = sp.sum(pore_dt == 1)
im_w_throats = spim.binary_dilation(input=pore_im, structure=struc_elem(1))
im_w_throats = im_w_throats*sub_im
Pn = sp.unique(im_w_throats)[1:] - 1
for j in Pn:
References
----------
[1] Gostick, J. "A versatile and efficient network extraction algorithm
using marker-based watershed segmenation". Physical Review E. (2017)
"""
peaks = sp.copy(peaks)
if dt.ndim == 2:
from skimage.morphology import square as cube
else:
from skimage.morphology import cube
labels, N = spim.label(peaks)
slices = spim.find_objects(labels)
for i in range(N):
s = extend_slice(s=slices[i], shape=peaks.shape, pad=10)
peaks_i = labels[s] == i+1
dt_i = dt[s]
im_i = dt_i > 0
iters = 0
peaks_dil = sp.copy(peaks_i)
while iters < max_iters:
iters += 1
peaks_dil = spim.binary_dilation(input=peaks_dil,
structure=cube(3))
peaks_max = peaks_dil*sp.amax(dt_i*peaks_dil)
peaks_extended = (peaks_max == dt_i)*im_i
if sp.all(peaks_extended == peaks_i):
break # Found a true peak
elif sp.sum(peaks_extended*peaks_i) == 0:
peaks_i = False
break # Found a saddle point
that the surface area of region 1 is stored in element 0 of the list.
"""
print('_'*60)
print('Finding surface area of each region')
im = regions.copy()
# Get 'slices' into im for each pore region
slices = spim.find_objects(im)
# Initialize arrays
Ps = sp.arange(1, sp.amax(im)+1)
sa = sp.zeros_like(Ps, dtype=float)
# Start extracting marching cube area from im
for i in tqdm(Ps):
reg = i - 1
if slices[reg] is not None:
s = extend_slice(slices[reg], im.shape)
sub_im = im[s]
mask_im = sub_im == i
mesh = mesh_region(region=mask_im, strel=strel)
sa[reg] = mesh_surface_area(mesh)
result = sa * voxel_size**2
return result
p_area_surf = sp.zeros((Np, ), dtype=int)
mc_sa = sp.zeros((Np, ), dtype=int)
t_area_mc = []
t_conns = []
t_dia_inscribed = []
t_area = []
t_perimeter = []
t_coords = []
dt_shape = sp.array(dt.shape)
# Start extracting size information for pores and throats
for i in tqdm(Ps):
pore = i - 1
# if slices[pore] is None:
# continue
s = extend_slice(slices[pore], im.shape)
sub_im = im[s]
sub_dt = dt[s]
pore_im = sub_im == i
# ---------------------------------------------------------------------
padded_mask = sp.pad(pore_im, pad_width=1, mode='constant')
if sp.any(dt_shape == 1):
mask_2d = sp.pad(pore_im.squeeze(), pad_width=1, mode='constant')
ax = sp.where(dt_shape == 1)[0][0]
pore_dt = spim.distance_transform_edt(input=mask_2d.squeeze())
pore_dt = sp.expand_dims(pore_dt, ax)
else:
pore_dt = spim.distance_transform_edt(padded_mask)
if padded_mask.ndim == 3:
filter_mask = spim.convolve(padded_mask*1.0,
weights=ball(1))/sp.sum(ball(1))
verts, faces, norm, val = measure.marching_cubes_lewiner(filter_mask)
Elementary Volume (Rev), Advances in Transport Phenomena in Porous Media
(1987)
"""
im_temp = sp.zeros_like(im)
crds = sp.array(sp.rand(npoints, im.ndim)*im.shape, dtype=int)
pads = sp.array(sp.rand(npoints)*sp.amin(im.shape)/2+10, dtype=int)
im_temp[tuple(crds.T)] = True
labels, N = spim.label(input=im_temp)
slices = spim.find_objects(input=labels)
porosity = sp.zeros(shape=(N,), dtype=float)
volume = sp.zeros(shape=(N,), dtype=int)
for i in tqdm(sp.arange(0, N)):
s = slices[i]
p = pads[i]
new_s = extend_slice(s, shape=im.shape, pad=p)
temp = im[new_s]
Vp = sp.sum(temp)
Vt = sp.size(temp)
porosity[i] = Vp/Vt
volume[i] = Vt
profile = namedtuple('profile', ('volume', 'porosity'))
profile.volume = volume
profile.porosity = porosity
return profile
s_coords = sp.zeros((Ns, s_im.ndim), dtype=float)
s_volume = sp.zeros((Ns, ), dtype=float)
s_dia_local = sp.zeros((Ns, ), dtype=float)
s_dia_global = sp.zeros((Ns, ), dtype=float)
s_label = sp.zeros((Ns, ), dtype=int)
s_area_surf = sp.zeros((Ns, ), dtype=int)
st_conns = []
st_dia_inscribed = []
st_area = []
st_perimeter = []
st_coords = []
solid_pore_area_surf = sp.zeros((Ns, ), dtype=float)
for k in tqdm(Ss):
solid = k - 1
ext_s_slice = extend_slice(s_slice[solid], s_im.shape)
sub_im_s = s_im[ext_s_slice] # Sub image of solid extended slice
sub_dt_s = s_dt[ext_s_slice] # Sub dt of solid ex_slice
# Surface Area Calculation
sub_im_s_dt = spim.distance_transform_edt(sub_im_s)
sub_im_s_mask = sub_im_s == k
sub_im_s_solid_labels = sub_im_s_mask*sub_im_s_dt
solid_pore_area_surf[solid] = sp.sum(sub_im_s_solid_labels == 1)
im_solid_s = sub_im_s == k
dt_solid_s = spim.distance_transform_edt(sp.pad(im_solid_s,
pad_width=1,
mode='constant'))
solid_offset = sp.array([kx.start for kx in ext_s_slice])
s_label[solid] = k
s_coords[solid, :] = spim.center_of_mass(im_solid_s) + solid_offset
p_solid_area_surf = sp.zeros((Np, ), dtype=int)
p_solid_volume = sp.zeros((Np, ), dtype=int)
pore_solid_conns = []
# Start extracting size information for pores and throats
for i in tqdm(Ps):
pore = i - 1
# if slices[pore] is None:
# continue
s = extend_slice(slices[pore], p_im.shape)
sub_im = p_im[s]
sub_dt = dt[s]
if dual is True:
# This calculates chunk of solid properties connected with pore
ext_p_on_s_slice = extend_slice(p_on_s_slice[pore], p_im.shape)
# Sub image of solid extended slice
sub_sim = p_on_s[ext_p_on_s_slice]
# Sub distance transform of solid extended slide
sub_sdt = spim.distance_transform_edt(sub_sim)
sub_sdt = sub_sdt == 1 # Finding only surface region of solid
# Detecting chunk of solid connected with ith pore
solid_im = (sub_sdt*sub_sim) == i
p_solid_area_surf[pore] = sp.sum(solid_im)
p_solid_volume[pore] = sp.sum(sub_sim == i)
# Finding Pore Solid Neighbours
pore_regions_full = pore_regions[ext_p_on_s_slice]
pore_region_mask = pore_regions[ext_p_on_s_slice] == i
solid_regions_full = solid_regions[ext_p_on_s_slice]
solid_labels_on_pore = ((pore_region_mask * pore_regions_full !=
0) * solid_regions_full)
p_label = sp.zeros((Np, ), dtype=int)
p_area_surf = sp.zeros((Np, ), dtype=int)
t_conns = []
t_dia_inscribed = []
t_area = []
t_perimeter = []
t_coords = []
p_solid_area_surf = sp.zeros((Np, ), dtype=int)
p_solid_volume = sp.zeros((Np, ), dtype=int)
pore_solid_conns = []
# Start extracting size information for pores and throats
for i in tqdm(Ps):
pore = i - 1
# if slices[pore] is None:
# continue
s = extend_slice(slices[pore], p_im.shape)
sub_im = p_im[s]
sub_dt = dt[s]
if dual is True:
# This calculates chunk of solid properties connected with pore
ext_p_on_s_slice = extend_slice(p_on_s_slice[pore], p_im.shape)
# Sub image of solid extended slice
sub_sim = p_on_s[ext_p_on_s_slice]
# Sub distance transform of solid extended slide
sub_sdt = spim.distance_transform_edt(sub_sim)
sub_sdt = sub_sdt == 1 # Finding only surface region of solid
# Detecting chunk of solid connected with ith pore
solid_im = (sub_sdt*sub_sim) == i
p_solid_area_surf[pore] = sp.sum(solid_im)
p_solid_volume[pore] = sp.sum(sub_sim == i)
' unexpected behavior.')
if im.ndim == 2:
cube = square
ball = disk
# Get 'slices' into im for each region
slices = spim.find_objects(im)
# Initialize arrays
Ps = sp.arange(1, sp.amax(im)+1)
sa = sp.zeros_like(Ps, dtype=float)
sa_combined = [] # Difficult to preallocate since number of conns unknown
cn = []
# Start extracting area from im
for i in tqdm(Ps):
reg = i - 1
if slices[reg] is not None:
s = extend_slice(slices[reg], im.shape)
sub_im = im[s]
mask_im = sub_im == i
sa[reg] = areas[reg]
im_w_throats = spim.binary_dilation(input=mask_im,
structure=ball(1))
im_w_throats = im_w_throats*sub_im
Pn = sp.unique(im_w_throats)[1:] - 1
for j in Pn:
if j > reg:
cn.append([reg, j])
merged_region = im[(min(slices[reg][0].start,
slices[j][0].start)):
max(slices[reg][0].stop,
slices[j][0].stop),
(min(slices[reg][1].start,
slices[j][1].start)):