Generating an image of a sphere packing from a list of center points and radii

[jgostick]

There have been several requests for this lately, so I cooked up a prototype function. We'll be adding something like to PoreSpy in the next release I suppose but for now you can use this:

import numpy as np
from porespy.tools import insert_sphere

def custom_sphere_packing(d, voxel_size, shape=None):
    r"""
    Generates a sphere packing given a list of centers and radii

    Parameters
    ----------
    df : DataFrame or dict
        The X, Y, Z center coordinates and R of each sphere in the packing
        should be stored in their own column or key. See ``Notes`` for more
        detail.
    voxel_size : float
        The resolution of the image to generate in terms of how many 'units'
        in the supplied packing data correspond to the side-length of a voxel.

    Returns
    -------
    spheres : ndarray
        A numpy ndarray of ``True`` values indicating the spheres and
        ``False`` elsewhere.

    Notes
    -----
    The input data should be in column format, either as a dictionary of 1D
    *numpy* arrays like this:

    .. code::

        d['X'] = np.array([1, 1, 2])
        d['Y'] = np.array([1, 3, 2])
        d['Z'] = np.array([1, 1, 1])
        d['R'] = np.array([0.5, 0.7, 0.6])

    Or a *pandas* ``DataFrame` like this:

    ==== ==== ==== ==== ====
    ID   'X'  'Y'  'Z'  'R'
    ==== ==== ==== ==== ====
    0    1    1    1    0.5
    1    1    3    1    0.7
    2    2    2    1    0.6
    ==== ==== ==== ==== ====

    """
    vx = voxel_size
    x = ((d['X'] - d['X'].min() + d['R'].max())/vx).astype(int)
    y = ((d['Y'] - d['Y'].min() + d['R'].max())/vx).astype(int)
    z = ((d['Z'] - d['Z'].min() + d['R'].max())/vx).astype(int)
    try:  # Add a try/except to catch files with no radius to use r=1 voxel
        r = (d['R']/vx).astype(int)
    except KeyError:
        r = np.ones_like(x)
    R = r.max()
    shape = np.ceil([x.max()+R+1, y.max()+R+1, z.max()+R+1]).astype(int)
    im = np.zeros(shape, dtype=bool)
    for i in range(len(x)):
        im = insert_sphere(im, c=[x[i], y[i], z[i]], r=r[i], v=True)
    return im

It would be used like this:

import pandas as pd
import matplotlib.pyplot as plt
import porespy as ps

df = pd.read_csv('finney_packing.csv')
df['R'] = 1.5
im = custom_sphere_packing(df, voxel_size=0.18)
plt.imshow(ps.visualization.sem(~im, axis=0))

Where the finney_packing.csv file was obtained from the digitalrocksportal.org but is also uploaded here: finney_packing.csv

[jgostick]

Note that the above function needs a bit more work. For instance, it can't handle 2D images. I also has no check to enure the generated image is not too large. It might make more sense if the shape argument was passed and voxel_size was computed internally...then the image would never get larer than desired.

[Soheil1993]

Dear Dr. @jgostick

Many thanks for your kind consideration and valuable time.

If I am not mistaken, we should replace im = insert_sphere(...) by im = ps.tools.insert_sphere(...) in the code to get the result.

Best regards,
Soheil

[jgostick]

Yes, good point! I was running this script in the same file as the insert_sphere function so didn't need to import it. I will update the function.

[Soheil1993]

Dear Dr. [@jgostick]

I used this code for finney_packing.csv and it worked properly. Aslo I could postprocess the results using Paraview that is my desired goal. (attached image)

I obtained the position (X,Y,Z) and radii(R) of particles from Yade in a text file and then I converted it to a csv file.

But it did not worked and I could not obtain any result.I attached the text file and csv file and the code for your consideration.

I think the problem is it can not recognize my packing. Because I gained this results:
‘‘ Initial number of peaks: 1
Peaks after trimming saddle points: 1
Peaks after trimming nearby peaks: 1 ‘’

It should be mentioned that I also tried sphere packing with the same and different sizes of particles. Then I removed the R column and it did not work as well, because there is no R (radii) data in the finney_packing.csv.

I am grateful for any hints or insights into this issue.

Thank you so much in advance and I look forward to hearing from you.

Best regards,
Soheil

Soheil.zip

[jgostick]

Try putting voxel_size = 0.018 or even 0.01. This will make an image of size ~300 cubed, which will presumably create enough room for the particles to not overlap.

[jgostick]

Also, try without the ID column in your spreadsheet. I don't think my prototype function accounts for that

[Soheil1993]

Dear Dr.@jgostick

Thank you so much for your kind reply and valuable information.

The finney_packing.csv file includes ID,X,Y and Z columns and there is no R data. But I removed the ID column in spreadsheet according to your recommendation.

I also used the (#678 (comment)) text file and received the same results as well.

It should be mentioned that I am using SNOW network extraction (advanced) code. (https://porespy.org/examples/networks/tutorials/snow_advanced.html) it works for finney_packing.csv .

I attached the results and errors for your consideration.

Best regards,
Soheil

[jgostick]

I noticed the ['R'] problem and updated the example function just recently. If R is not found it will assume 1 pixel. The problem you are getting is that you put R as 1.5, in the same units as X, Y, and Z. If you look at the code you can see that it normalizes the values by the voxel_resolution, so the R of 1.5 results in a sphere of radius 1.5/0.001 or 150 voxels, which is the whole image. Try R=0.01, and voxel_size of 0.001, and you'll get spheres with radius of 10 pixels.

BTW, posting the code above was meant as an example for you to learn from while you implement your own function. I cannot keep giving you this basic tech support while you learn the ins-and-outs of volumetric image analysis. Sorry to be grumpy, but you could have answered your own questions if you just stopped to look at what the code is actually doing.

[Soheil1993]

Dear Dr. @jgostick

Thank you so much for your detailed description.

I am so sorry for late reply and these simple questions. I appreciate your patience and kind consideration :)

I tried R=0.01, and voxel_size of 0.001 but I am faced with this error: “Kernel Restarting: The kernel appears to have died.It will restart automatically.”
I waited one hour for that but It did not work. Is there another way to do it?

Would you please let me know if there is any option to extract pore networks based on particle locations and radii?

[Soheil1993]

Dear Dr. @jgostick

I tried R=0.03, and voxel_size of 0.003 and fortunately it worked :)

But in paraview the pore diameters are equal. I do not know what the mistake is.

The throat diameters are different

I attached the particle packing for your consideration.

I am grateful for any hints or insights into this issue.

[jgostick]

Great! The kernel dying was no doubt because the image was too big. The function should probably take a shape argument, which the desired size of the final image, then all the values get scaled accordingly. Anyway, good luck with this now that you have what you need.

[Soheil1993]

Dear Dr. @jgostick

Thank you very much for your valuable time spent on solving my issue.

Do you have any recommendations for solving the equal pore diameter sizes?

Should I change any other value in SNOW network extraction code?

One more thing, Is it necessary to multiply all the parameters (X,Y,Z,R) by 1e-6 for changing m to um?

[Soheil1993]

Dear Dr. @jgostick

So sorry for bothering you again. I tried several times with different packing and different R and voxel sizes, and in all cases the pore diameters “network/net 01/properties/pore.diameter” are equal (In paraview). I could not find the issue.

I gain different values for pores just when I select “network/net 01/properties/pore._id." What is "pore._id"?

Would you please let me know what is the difference between the pore network items?

For obtaining the pore size distribution, which one should I select?

[jgostick]

Perhaps I am mis-understading, but in one of your posts above, you typed:

I tried R=0.03, and voxel_size of 0.003 and fortunately it worked :)

And now you're saying that all the particles are the same size...presumably becuase you set all the R values to 0.03? You need to put a column of R values into the csv file, with the unique R value for each sphere, which also has it's unique x, y, z position.

[Soheil1993]

Dear Dr. @jgostick

Sorry for asking too many questions and bothering you with simple questions.
I said that when I use particles with different radii values ‘R’ I gain equal “pore” diameters in paraview. Fortunately, I found the issue and it works properly.

Thank you very much for your kind reply and valuable time.
Wish you all the best.