Visualization of Pore Network Models Using ParaView

[je-ha-jung]

Hello,
I am currently using CT images to create pore network models and calculate absolute permeability.
For the images, I used 1000 slices of the Bentheimer sandstone dataset from the Digital Rocks Portal. Initially, I used the SNOW2 algorithm, but due to memory limitations during repeated runs, I switched to using the parallel implementation of snow partitioning in PoreSpy. I set the sigma value to 0.4 initially, but it seemed to overly segment the pores, so I changed it to 1.

Here are the challenges I am facing:
The shape of the pore network model seems to be well-constructed, but the equivalent diameter values are extremely small. The pores and throats appear to be formed very narrowly, resulting in very low absolute permeability values. How can I adjust the model to better reflect the actual pore sizes of the image?
I have tried slightly overestimating the binarized images and tweaking the sigma parameter, but the pore sizes in the model still seem much smaller than the actual pore sizes in the rock.
Since the pores and throats were constructed based on the watershed algorithm, I’m not sure if this is the reason they appear so small. Could you please help me resolve this issue? Below is my code.

import os
import glob
import imageio.v2 as imageio
import porespy as ps
import numpy as np
import openpnm as op
import matplotlib.pyplot as plt
from skimage.segmentation import watershed
from scipy import ndimage as spim
from porespy.tools import randomize_colors
from porespy.filters import find_peaks, trim_saddle_points, trim_nearby_peaks
from PIL import Image
from vtk.util import numpy_support
import vtk
folder_path = r"C:\Users\user\Desktop\lab\포어 네트워크 모델\bentheimer_otsubinary_image_1000_tif"

tif_files = glob.glob(os.path.join(folder_path, "*.tif"))
tif_files.sort()

sample_img = Image.open(tif_files[0])
width, height = sample_img.size

stack = np.zeros((len(tif_files), height, width), dtype=np.uint8)
for i, tif_file in enumerate(tif_files):
img = Image.open(tif_file)
stack[i, :, :] = np.array(img)
print("3차원 배열의 shape:", stack.shape)

im = stack
resolution = 0.00000225 # 2.25 μm = 0.00225 mm
print(f"Resolution (mm/voxel): {resolution}")

snow_out = ps.filters.snow_partitioning_parallel(
im=im, # 이진화된 CT 이미지 (3D numpy 배열)
divs=2, # x, y, z 방향으로 2개씩 분할
r_max=10, # 최대 공극 반지름
sigma=1 # Gaussian 필터의 sigma 값
)
fig, ax = plt.subplots(1, 1, figsize=[6, 6])
ax.imshow(randomize_colors(snow_out.regions[:, :, snow_out.regions.shape[2] // 2]), cmap='jet') # 중간 슬라이스
ax.set_title('Segmented Image')
plt.show()

fig, ax = plt.subplots(1, 2, figsize=(12, 6))
ax[0].imshow(im[:, :, im.shape[2] // 2], cmap='gray', origin='lower')
ax[0].set_title("Original Image (Central Slice)")
ax[1].imshow(randomize_colors(snow_out.regions[:, :, snow_out.regions.shape[2] // 2]), cmap='jet', origin='lower')
ax[1].set_title("Segmented Image (Central Slice)")
plt.tight_layout()
plt.show()

network_dict = ps.networks.regions_to_network(
regions=snow_out.regions,
voxel_size=resolution
)
pn = op.io.network_from_porespy(network_dict)
print("OpenPNM Network Summary:")
print(pn)
print("\nNetwork Keys:")
print(pn.keys())
print("\nExample Pore Coordinates:")
print(pn['pore.coords'][:5])
print("\nExample Throat Connections:")
print(pn['throat.conns'][:5])

im = ps.tools.align_image_with_openpnm(im)
imageio.volsave('image.tif', np.array(im, dtype=np.int8))
op.io.project_to_vtk(project=pn.project)
output_vtk_file = r"C:\Users\user\Desktop\lab\포어 네트워크 모델\bentheimer_otsubinary_image_1000_tif\sigma1_partitioning_network.vtk"
op.io.project_to_vtk(project=pn.project, filename=output_vtk_file)

max_diameter = pn['pore.equivalent_diameter'].max()
min_diameter = pn['pore.equivalent_diameter'].min()

print(f"Maximum pore equivalent diameter: {max_diameter}")
print(f"Minimum pore equivalent diameter: {min_diameter}")

Maximum pore equivalent diameter: 0.00027545349464088
Minimum pore equivalent diameter: 2.7915772090473024e-06