Extracting Pores from a CT Image and Visualizing in ParaView

[je-ha-jung]

Dear Professor,

I have learned from various experiments and modeling using CT images with OPENPNM, and I realized that its utility is very high. I am currently studying OPENPNM and would like to extract pores from my CT image and visualize them using ParaView.

However, as a beginner with OPENPNM, I used the snow extraction advanced method from Porespy as it is. The tutorial used artificial images, but I have converted my CT scan stack (515 TIFF images) into a numpy.ndarray format. To verify if ParaView can load 3D images, I converted the numpy.ndarray file into a VTK file, and upon opening it in ParaView, I was able to see it in 3D, as shown in the image below. Since I need to extract pores and create a pore network model from this 3D image, I used the same code as in the tutorial.

However, unlike the tutorial, I am encountering an issue where all peak values are coming out as 1, which prevents the pore network model from being constructed.

Could you please help me understand what steps I might be missing? I look forward to your detailed guidance.

Thank you.

This image does not include pore extraction; it only represents the 3D visualization of the image.

import imageio
import numpy as np
import porespy as ps
import openpnm as op
import scipy.ndimage as spim
import matplotlib.pyplot as plt
from porespy.filters import find_peaks, trim_saddle_points, trim_nearby_peaks
from porespy.tools import randomize_colors
from skimage.segmentation import watershed

ps.settings.tqdm['disable'] = True
ps.visualization.set_mpl_style()
np.random.seed(10)

file_path = r"C:\Users\user\Desktop\lab\cube_515\3d_image.npy"
im = np.load(file_path)
im.shape, type(im)
fig, ax = plt.subplots()
ax.imshow(ps.visualization.show_planes(im), cmap=plt.cm.bone);

sigma = 0.4
dt = spim.distance_transform_edt(input=im)
dt1 = spim.gaussian_filter(input=dt, sigma=sigma)
peaks = find_peaks(dt=dt)

print('Initial number of peaks: ', spim.label(peaks)[1])
peaks = trim_saddle_points(peaks=peaks, dt=dt1)
print('Peaks after trimming saddle points: ', spim.label(peaks)[1])
peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
peaks, N = spim.label(peaks)
print('Peaks after trimming nearby peaks: ', N)

regions = watershed(image=-dt, markers=peaks, mask=dt > 0)
regions = randomize_colors(regions)

fig, ax = plt.subplots()
ax.imshow((regions*im), cmap=plt.cm.nipy_spectral);

net = ps.networks.regions_to_network(regions*im, voxel_size=1)
pn = op.io.network_from_porespy(net)
print(pn)

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)

[jgostick]

I think the solution is simple: You need to convert your image from greyscale to binary. The numpy array needs to be True for voxels which are void and False for voxels which are solid. This is called segmentation. It's easy in principle, but takes some effort. Hopefully this tutorial helps.

[je-ha-jung]

Dear Professor,

Following your advice on binarizing the images, I was able to visualize them using ParaView. Thank you so much for your help! However, I’ve encountered some issues. The sizes of the pores appear to be uniformly set, and the throats seem insufficient to connect all the pores. Since I am using CT images of basalt, most of the pores should be interconnected rather than isolated, but in the visualized image, it looks like many pores are not connected by throats. I'm not sure what I might be missing.

I thought it might be due to using a voxel size of 1, so I calculated the actual voxel size of the TIF image as 0.055 and visualized it again, but the overall diameter was still visualized as uniform.

This phrase should clearly convey your intention.

In the previous code, I only modified the binarization process and adjusted the voxel size from 1 to 0.055. Since there seem to be more pores compared to throats, I also adjusted the sigma from 0.4 to 0.8. The binarized image was generated correctly, and the watershed process also seems to have proceeded well.

I look forward to your response. Thank you.

[jgostick]

Have you looked at the paraview tutorial to get the sizes of the pores and throats represented properly?

https://openpnm.org/examples/tutorials/10_visualization_options.html

[je-ha-jung]

Despite visualizing based on the tutorial, I was unable to resolve my concerns. I am still struggling with understanding why the throats do not connect all the pores and why, even though I input the actual voxel size of 0.055 instead of 1 in the scaling adjustment section, the pore sizes are visualized as if they are the same.

[jgostick]

I'm not sure what else to recommend. Is it possible that the image you're seeing is correct and that your image has a lot of isolated pores? Also, did you ensure that the "Masking/Glyph Mode" is set to 'all points' in paraview?

[je-ha-jung]

Dear Professor,

I hope this message finds you well. I greatly appreciate the assistance you’ve provided so far. However, I am still encountering some difficulties.

Although I have managed to adjust the pore sizes successfully and converted them to equivalent diameters, the pore colors remain uniform. I intended to color the pores according to their diameters, but it seems that the coloring might have been reversed, as it appears to be applied to the throats instead.

Additionally, contrary to the tutorial instructions, many of my pores still seem not to be connected properly.

Could you please review the attached code and help me address these issues? I am still a beginner and am finding it challenging, but I am working diligently to overcome these obstacles. Your guidance would be greatly appreciated.

Thank you very much for your time and support.

import os
import imageio
import numpy as np
import porespy as ps
import openpnm as op
import scipy.ndimage as spim
import matplotlib.pyplot as plt
from porespy.filters import find_peaks, trim_saddle_points, trim_nearby_peaks
from porespy.tools import randomize_colors
from skimage.segmentation import watershed

input_image_path = r'C:\Users\user\Desktop\lab\basalt_images'

print(f"Checking file at path: {input_image_path}")
print(f"File exists: {os.path.isfile(input_image_path)}")

folder_path = r"C:\Users\user\Desktop\lab\basalt_images"
file_names = sorted([f for f in os.listdir(folder_path) if f.endswith('.tif')])

output_folder = r"C:\Users\user\Desktop\lab\basalt_images_binary"
if not os.path.exists(output_folder):
os.makedirs(output_folder)

for file_name in file_names:
image_path = os.path.join(folder_path, file_name)
im = imageio.imread(image_path)

smoothed_im = filters.gaussian(im, sigma=1)

inverted_im = np.max(smoothed_im) - smoothed_im

thresh = np.percentile(inverted_im, 75.37)
binary_im = inverted_im > thresh

fig, ax = subplots(figsize=[5, 5])
ax.imshow(binary_im, cmap='gray')
plt.show()

output_path = os.path.join(output_folder, f"binary_{file_name}")
imageio.imwrite(output_path, binary_im.astype(np.uint8) * 255)
tif_files = glob.glob(os.path.join(folder_path, "*.tif"))

tif_files.sort()

stack = np.zeros((len(tif_files), 105, 105), dtype=np.uint16)

for i, tif_file in enumerate(tif_files):
img = Image.open(tif_file)
stack[i, :, :] = np.array(img)

print("3차원 배열의 shape:", stack.shape)

np.save(os.path.join(folder_path, "3d_image.npy"), stack)

fig, ax = plt.subplots()
ax.imshow(ps.visualization.show_planes(im), cmap=plt.cm.bone);

sigma = 1
dt = spim.distance_transform_edt(input=im)
dt1 = spim.gaussian_filter(input=dt, sigma=sigma)
peaks = find_peaks(dt=dt)

print('Initial number of peaks: ', spim.label(peaks)[1])
peaks = trim_saddle_points(peaks=peaks, dt=dt1)
print('Peaks after trimming saddle points: ', spim.label(peaks)[1])
peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
peaks, N = spim.label(peaks)
print('Peaks after trimming nearby peaks: ', N)

regions = watershed(image=-dt, markers=peaks, mask=dt > 0)
regions = randomize_colors(regions)
fig, ax = plt.subplots()
ax.imshow((regions*im)[:, :, 100], cmap=plt.cm.nipy_spectral);

net = ps.networks.regions_to_network(regions*im, voxel_size=0.055)
pn = op.io.network_from_porespy(net)
print(pn)

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\basalt_images_binary\final_network.vtk"
op.io.project_to_vtk(project=pn.project, filename=output_vtk_file)

[jgostick]

How about you send me your image and I'll see what's wrong with them. Put the images on dropbox (or equivalent) and send a link to my email, either uwaterloo or gmail works.

[je-ha-jung]

Thank you very much for your assistance. I used basalt CT images as input data. I have created two folders: one contains the original TIFF images, and the other contains the binarized images. Inside the folder with the binarized images, there is also the VTK file I used for visualization, named final_network.vtp.

Could you please check if there are any steps I missed during the code implementation or the visualization modeling process? I want to ensure that the visualization I conducted was done correctly. Thank you, Professor.

https://www.dropbox.com/scl/fo/b04fniteea6is7s5ar2hz/ABSv4o3-8fER_Y7JD2DBuVg?rlkey=cth423wrewud3igx9y889rubb&e=1&st=nsdyayib&dl=0

[je-ha-jung]

Thank you very much for your assistance. I used basalt CT images as input data. I have created two folders: one contains the original TIFF images, and the other contains the binarized images. Inside the folder with the binarized images, there is also the VTK file I used for visualization, named final_network.vtp. Could you please check if there are any steps I missed during the code implementation or the visualization modeling process? I want to ensure that the visualization I conducted was done correctly. Thank you, Professor. https://www.dropbox.com/scl/fo/b04fniteea6is7s5ar2hz/ABSv4o3-8fER_Y7JD2DBuVg?rlkey=cth423wrewud3igx9y889rubb&e=1&st=nsdyayib&dl=0 2024년 8월 28일 (수) 오전 6:39, jgostick ***@***.***>님이 작성:

…

How about you send me your image and I'll see what's wrong with them. Put the images on dropbox (or equivalent) and send a link to my email, either uwaterloo or gmail works. — Reply to this email directly, view it on GitHub <#2936 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BKVCKR3U6H5O4X3F3ZXYOP3ZTTWZRAVCNFSM6AAAAABM2GHNCOVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANBWG44DQOI> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[jgostick]

I have taken a look at your images, and I think the problem starts there. They are too small and their resolution is too low. Maybe you have only uploaded 'low res' version, but if these are the only images you have then they are not sufficient. Because the resolution is so low, you cannot see any connections between the pore bodies. This is why you are getting only isolated pores...because this is what your images has. Assuming you can get higher resolution images, I am still not sure that you will have success. It looks to me like these pore bodies are almost completely isolated from their neighbors. I have never worked with Basalt, but perhaps this is one of their main characteristics? What is the permeability of these materials? I am guess it is quite low, meaning that these pore bodies are not well connected.

[je-ha-jung]

I understood which part I was missing. I thought I had to zoom in to see the connectivity of the pore closely, so I think I used low quality and a small number of images. I think I'll have to re-implement it with the number of photos and high quality images. I was wondering if shape 300, 300, 300 would be okay, so I'm going to re-produce it and try it, but I'm glad if there's nothing I missed in the code I implemented. When I only looked at the pore network model, the connectivity seemed to be low, so I didn't get the transmittance. 2024년 8월 29일 (목) 오전 2:40, jgostick ***@***.***>님이 작성:

…

I have taken a look at your images, and I think the problem starts there. They are too small and their resolution is too low. Maybe you have only uploaded 'low res' version, but if these are the only images you have then they are not sufficient. Because the resolution is so low, you cannot see any connections between the pore bodies. This is why you are getting only isolated pores...because this is what your images has. Assuming you can get higher resolution images, I am still not sure that you will have success. It looks to me like these pore bodies are almost completely isolated from their neighbors. I have never worked with Basalt, but perhaps this is one of their main characteristics? What is the permeability of these materials? I am guess it is quite low, meaning that these pore bodies are not well connected. — Reply to this email directly, view it on GitHub <#2936 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BKVCKR2AKD6BL2P2KGDDKV3ZTYDRNAVCNFSM6AAAAABM2GHNCOVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANBXHAZDOMA> . You are receiving this because you authored the thread.Message ID: ***@***.***>