How to extract a network from a very large image?

[JoHommel]

I have trouble to extract a single PN from my µCT-data of column experiments, likely due to the size of the data set and memory limitations, but I would like to have a single PN for each experiment. I searched the porespy documentation, but I could not find any function that looks as if it would do the job.
In principal, it should be simply searching for pore bodies with identical/similar coordinates, merging them, and checking and updating the pore throats, if necessary.
Doing this by hand seems tedious and prone to errors. I am grateful for any hints or insights into this issue.

[jgostick]

I am not sure I understand your problem. I assume the pore structure is not changing, hence you expect a single PNM to fit all the images. So why can't you just extract the PNM from the first experiment and use that for all the others?

[jgostick]

wow, 2000^3 is pretty big! But you may be in luck. We have recently developed a parallelized watershed algorithm that can actually process the image in serial chunks. I am just realizing now that we don't have an example showing its use, but it actually runs automatically when you use ps.networks.snow2. We do have an example on using the snow_parallel function which finds the watershed, but does not return a network. Once you have the watershed you can pass it to porespy.networks.regions_to_network to get a network. This code is pretty new and has not been stress tested, so please keep me updated with any bumps you experience.

[JoHommel]

I tried using snow_parallel: (snow_out = ps.filters.snow_partitioning_parallel(binaryframes, divs=2, r_max=5, sigma=0.4)) and for a resized version of the image set (approx. 1000^3), the partitioning seems to work.
Also creating a PNM out of the regions seems to work:
net = ps.networks.regions_to_network(snow_out.regions, voxel_size=voxel_size)
...
pn = op.network.GenericNetwork()
pn.update(net)

op.io.VTK.save(network = pn, filename=(outputfile +'.vtk'))

The PNM is also successfully saved as a vtk file for visualization on Paraview, but somehow, writing the PNM to a csv file fails:

df = Pandas.to_dataframe(network=pn, phases=[], join=True, delim =',')
df.to_csv((outputfile +'.csv'))

With the erroro message:
File "PNGtoPN_3D_SNOW_2021_11_3_parallel.py", line 171, in
df = Pandas.to_dataframe(network=pn, phases=[], join=True, delim =',')
File "/.../anaconda3/lib/python3.7/site-packages/openpnm/io/Pandas.py", line 76, in to_dataframe
tmp = np.split(arr, arr.shape[1], axis=1)
IndexError: tuple index out of range

And this error message is not that helpful so that at the moment I do not even understand which index is supposedly out of bounds.

[jgostick]

So this sounds like mostly good news! For Paraview, I recommend the XDMF format, which is WAY faster. It creates two files: XDMF and HDF5. The XDMF file is similar to the VTP/VTK file, but instead of containing data it points to data in the accompanying HDF5 file.

For the CSV error, I'l lhave to investigate. It looks like you're creating a dataframe, then converting to CSV, so maybe try using the openpnm CSV class instead?

[JoHommel]

Thank you for the reply. I am still a bi new to python and mainly reusing scripts created by other colleagues. I never questioned their choice of how to write files and such basic things. I'll try your suggestions instead and see how that works.

[JoHommel]

I finally got some time to trying the openpnm CSV class and it gives me the same error (IndexError: tuple index out of range).

[JoHommel]

Further, looking at the network in paraview, it seems as if there is some kind of boundary or discontinuity within the PNM. Is there maybe something gone wrong with stitching together the two PNM generated in snow_parallel: (snow_out = ps.filters.snow_partitioning_parallel(binaryframes, divs=2, r_max=5, sigma=0.4)) ?

However, if this due to stiching, I would expect the discontinuity to be rather somewhat in the middle and not at roughly one third to two thirds along the length of the image stack.

[jgostick]

are you sure that discontinuity is not in the image? Did some slices get lost somehow? I'm not saying porespy/openpnm is infallible, but the location of that discontinuity seems odd. You used div=2, so that that would put the divisions at 1/3 and 2/3 along each axis. That discontinuity might be at 1/3, but why isn't there one at 2/3, and why not the other 2 axes?

[JoHommel]

It is quite difficult to check in the images, as I do not have a 3D rendering of the images and it is in the direction of the stack. But it is good to know that divs is the number of divisions, not the number of domains as I understood it from the documentation https://porespy.org/modules/generated/porespy.filters.snow_partitioning_parallel.html:
divs (list or int) –
Number of domains each axis will be divided.
I checked the images yesterday and nothing really caught my eye. I hope to find more time for a thorough check today and I agree that something in the images is probably more likely for this singular discontinuity. I really could not see any in any other direction, so I guess that this is the more likely case.
I just thought maybe it has something to do with the problem of writing out the csv file, that we discussed above.

[JoHommel]

Update: I just found that somehow the image indices got messed up during the resizing from 2000³ to 1000³, which is why the inlet 1/3 of the images got somehow placed at the outlet. So, this was a false alarm and has nothing to do with porespy/openpnm.

[jgostick]

Thanks for letting me know about the documentation error. In fact, I may be wrong about the meaning of divs. I'll have to double check the code. Glad you found the problem though.

[falvlima]

Hi, Jeff.

I am trying to process some huge data as well.
Actually, I have reduced it to a 1677, 1799, 500 image.
By using the parallel processing, I could generate the dic file and also the prj, but when I am trying to export the data to read in paraview no file comes to the io.
I have tried also to export as a HDF5 file, but, even though the following message appears in spyder, no file came.

op.io.HDF5.export_data(prj)
Out[31]: <HDF5 file "proj_06.hdf" (mode r+)>

Thanks in advance

Also the following is the code I am using...

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
import json
import cv2
from skimage import io

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

im = io.imread('E:\ILLDATA\Filezilla\Prescan\FLVIG16bin3.tif')
print(im.shape)

#im = cv2.bitwise_not(im)
#im = im.astype(np.bool)
#fig, ax = plt.subplots()
#ax.imshow(ps.visualization.show_planes(im), cmap=plt.cm.bone);

sigma = 4

regions = ps.filters.snow_partitioning_parallel(
im=im, divs=4, r_max=5, sigma=4)

net = ps.networks.regions_to_network(regions.regions, voxel_size=15)

prj = op.io.PoreSpy.import_data(net)

proj = op.io.PoreSpy.load(net)
op.io.VTK.save(proj.network)

op.io.HDF5.export_data(prj)

[jgostick]

Hi @falvlima
Skip the last line above and open the VTK (actually it's a 'vtp' file). Once you open paraview you can follow the instruction which we accidentally deleted and have re-added here: https://github.com/PMEAL/OpenPNM/blob/readd_paraview/examples/notebooks/io/view_vtp_in_paraview.md

This line: op.io.VTK.save(proj.network), should write a file to your local working directory, which you can find with os.get_cwd() (import os first). Hopefully that helps, but if not let me know.

BTW, we know that our docs are lacking, but working on the code always seems to take all our time. Sorry.