How to extract a pore network from a rock CT image

[MingyueTianya1]

Hello，

I would like to ask how to extract holes and throat data from the CT image of the rock and import them into openpnm to generate a network model.

Thank you,
mingyuetianya

[ZeliangSu]

I dug it for a while too. To me, the doc wasn't too clear for this, but it is as simple as the following:

import porespy as ps
import numpy as np
import matplotlib.pyplot as plt

# create a random porous structure
stack = np.random.randn(500, 500, 500)
stack[stack>0.5] = 1
stack[stack != 1] = 0

# extract the PN
# if I understand well the doc, it floods a voxel until it touches the wall
local_thickness = ps.filters.local_thickness(stack, sizes=10) # pls see below y-axis
data = ps.metrics.pore_size_distribution(local_thickness, bins=100, log=False)

# plot it
voxel_size = 1 # unit: um
plt.bar(data.R*voxel_size, data.pdf, width=data.bin_widths) # pore size distribution
plt.plot(data.R*voxel_size, 1 - data.cdf, 'red') # accumulative
plt.xlabel('invasion size [um]')
plt.ylabel('pore size smaller than stated volume') # check this link: https://porespy.org/examples/filters/reference/porosimetry.html

[ZeliangSu]

Well, I'm afraid that what I posted is not what you search as I'm in the same case as yours. I just found another solution:

import porespy as ps
import openpnm as op
import numpy as np
import matplotlib.pyplot as plt

# create a random porous structure
stack = np.random.randn(500, 500, 500)
stack[stack>0.5] = 1
stack[stack != 1] = 0

# compute pore network
snow_output = ps.networks.snow2(stack, voxel_size=1)
try:
    pore_net, geo = op.io.PoreSpy.import_data(snow_output.network)
except AttributeError:
    pore_net, geo = op.io.from_porespy(snow_output.network)

voxel_size = 1 # unit: um
values, base = np.histogram(geo['pore.volume']*(4/3*np.pi*voxel_size**3), bins=100)
cumulative = np.cumsum(values)

fig, ax1 = plt.subplots()
ax1.bar(base[:-1], values, width=base[2]-base[1])
ax1.set_xlabel('pore volume [um^3]')
ax1.set_ylabel('counted')
ax2 = ax1.twinx()
ax2.plot(base[:-1], cumulative/cumulative.max()*100, 'red')
ax2.set_ylabel('accumulative volume fraction')
fig.show()

Note that if you use the geo.keys() method in the dictionary geo above, you can find all the information you need including: throat.diameter, pore.volume, pore.equivalent_diameter, etc. I will invite you to search the properties that you want to do the stat, accessing the dictionary by doing geo['some key word here']. Hope it helps!

[MingyueTianya1]

I am very sorry to give you a reply so late. Thank you very much for your reply to the problem, which benefit me a lot. The pore network has been successfully extracted.

…

------------------ 原始邮件 ------------------ 发件人: "PMEAL/porespy" ***@***.***>; 发送时间: 2022年5月10日(星期二) 晚上11:01 ***@***.***>; ***@***.******@***.***>; 主题: Re: [PMEAL/porespy] How to extract a pore network from a rock CT image (Discussion #688) Well, I'm afraid that what I posted is not what you search as I'm in the same case like yours. I just found another solution: import porespy as ps import openpnm as op import numpy as np import matplotlib.pyplot as plt # create a random porous structure stack = np.random.randn(500, 500, 500) stack[stack>0.5] = 1 stack[stack != 1] = 0 # compute pore network snow_output = ps.networks.snow2(stack, voxel_size=1) try: pore_net, geo = op.io.PoreSpy.import_data(snow_output.network) except AttributeError: pore_net, geo = op.io.from_porespy(snow_output.network) voxel_size = 1 # unit: um values, base = np.histogram(geo['pore.volume']*(4/3*np.pi*voxel_size**3), bins=100) cumulative = np.cumsum(values) fig, ax1 = plt.subplots() ax1.bar(base[:-1], values, width=base[2]-base[1]) ax1.set_xlabel('pore volume [um^3]') ax1.set_ylabel('counted') ax2 = ax1.twinx() ax2.plot(base[:-1], cumulative/cumulative.max()*100, 'red') ax2.set_ylabel('accumulative volume fraction') fig.show() Note that if you use the geo.keys() method in the dictionary geo above, you can find all the information you need including: throat.diameter, pore.volume, pore.equivalent_diameter, etc. I will invite you to search the properties that you want to do the stat, accessing the dictionary by doing geo['some key word here']. Hope it helps! — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>

[wangjiayun258138]

Can we discuss it, I want to learn from you

[MingyueTianya1]

Hello If you want to extract the pore network information, you can use the method of snow2 operation, and there are other examples on the official website to learn in many aspects, you are also welcome to join the nine five two seven six three four eight seven group to learn and exchange.

…

------------------ 原始邮件 ------------------ 发件人: "PMEAL/porespy" ***@***.***>; 发送时间: 2023年2月18日(星期六) 晚上8:04 ***@***.***>; ***@***.******@***.***>; 主题: Re: [PMEAL/porespy] How to extract a pore network from a rock CT image (Discussion #688) Can we discuss it, I want to learn from you — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>

[pengzhulin]

Hello, how can I make my CT results (such as .mesh or .stl files) to the 3D numpy array? Thank you!

[pengzhulin]

Hello，

I would like to ask how to extract holes and throat data from the CT image of the rock and import them into openpnm to generate a network model.

Thank you, mingyuetianya

Hello, how can I make my CT results (such as .mesh or .stl files) to the 3D numpy array? Thank you!

[jgostick]

that is actually impossible to answer without knowing what is inside your files. In fact, knowing what STL and MESH files typically describe, I don't like your chances of getting a numpy array out of them. You'll have to google and stackoverflow for a package that you can try.

[pengzhulin]

that is actually impossible to answer without knowing what is inside your files. In fact, knowing what STL and MESH files typically describe, I don't like your chances of getting a numpy array out of them. You'll have to google and stackoverflow for a package that you can try.

Thank you for your reply. I got my CT data for porous carbon materials. Actually, it is a series of tiff images. I use the software Dragonfly to segment and reconstruct the 3D object. Currently, I want to export this 3D object to my PC so I can use my OpenPNM scripts to conduct pore-scale simulation. However, the option I can choose is 'vtk', '.mesh' or '.stl'. These two is not easily used for OpenPNM. While, in OpenPNM, a 3D numpy array is always used. I am trying to convert this.

Another method I am trying is to use Dragonfly OpenPNM plugin, I export the networkx.classes.graph.Graph object by using the function named getNetworkX(), and then save it to my PC from the Dragonfly. Then, I load the object to openpnm, op.io.network_from_networkx(object_file). However, it fails. The network is different from what we normally use from the demos or generate from porespy, such as (1)op.network.Demo(shape=[25, 25, 1], spacing=1e-4) (2)pn = op.io.network_from_porespy(snow_output.network). I am currently trying to figure this out.

Do you have any suggestions for this question? I would be really appreciated. Thank you, Prof. Gostick.

[jgostick]

Sorry, I am not able to do Dragonfly tech support.

[pengzhulin]

Sorry, I am not able to do Dragonfly tech support.

I am sorry. I don't intend for that. I just use that software to process my CT images and get the geometry. Would you mind sharing any experience about creating Porespy array based on a series of CT images?

Thank you very much.

[jgostick]

Here is a recent tutorial on using snow2. You can also find some tutorials on actual materials here.

[pengzhulin]

Thank you very much. Your tutorials help a lot! Thank you!