pseudo_gravity_packing: periodic lateral boundaries and returning centers of added spheres

[heinsimon]

Hi,

I am thinking of using the function pseudo_gravity_packing to create some electrode structures.
I just saw the recent changes in fixing_pseudo_packing and was wondering if there is time/room for some more changes:

• Is it possible to allow for a periodic lateral boundary/edge for adding the spheres?
• Is it possible to return the centers and radius (maybe in the format to be used directly with points_to_spheres) of the added spheres?

Thank you!

Simon

[heinsimon]

Additional idea:

• Allow for a certain value (2*value?) to be set at overlaps of added spheres (only relevant in the case of clearance<0)

[jgostick]

1. Doing periodic edges seems interesting but periodicity is always such a pain. I will give it some thought.
2. Returning centers and radii is totally doable and quite easy. I will for sure add this to the PR before it gets merged.
3. This could be done but would require some changes deep in the _insert_disk_at_point function. At present there is a flag for overwrite=True/False but this could be changed to mode='keep|overwrite|add'.

[jgostick]

Out of curiosity, does this function work for you?:

from skimage.segmentation import watershed
from edt import edt
import pandas as pd
import scipy.ndimage as spim
import porespy as ps
import numpy as np

def spheres_to_coords(im):
    dt = edt(im)
    dt = spim.gaussian_filter(dt, sigma=0.4)
    pk = (spim.maximum_filter(dt, size=3) == dt)*im
    ws = watershed(-dt, spim.label(pk)[0])*im
    lt = ps.filters.local_thickness(im)
    xyz = np.vstack(np.where(pk)).T
    r = lt[pk.astype(bool)]
    df = pd.DataFrame()
    df['r'] = r
    df['x'] = xyz[:, 0]
    df['y'] = xyz[:, 1]
    if im.ndim == 3:
        df['z'] = xyz[:, 2]
    return df

[heinsimon]

1. Hopefully it is not to much hassle and will be added at some point. Thank you.
2. That would be nice.
3. Okay.

To your code example: It identifies the centers, but the radius is off:

from skimage.segmentation import watershed
from edt import edt
import pandas as pd
import scipy.ndimage as spim
import porespy as ps
import numpy as np
import copy

def spheres_to_coords(im):
    dt = edt(im)
    dt = spim.gaussian_filter(dt, sigma=0.4)
    pk = (spim.maximum_filter(dt, size=3) == dt)*im
    ws = watershed(-dt, spim.label(pk)[0])*im
    lt = ps.filters.local_thickness(im)
    xyz = np.vstack(np.where(pk)).T
    r = lt[pk.astype(bool)]
    df = pd.DataFrame()
    df['r'] = r
    df['x'] = xyz[:, 0]
    df['y'] = xyz[:, 1]
    if im.ndim == 3:
        df['z'] = xyz[:, 2]
    return df


imbox=np.full((100,50,50),True)

packing=ps.generators.pseudo_gravity_packing(imbox,r=15,clearance=0,seed=42,maxiter=10)

df=spheres_to_coords(packing)

packing2 = ps.generators.spheres_from_coords(df,mode='contained',smooth=True)

df_fix=copy.deepcopy(df)
df_fix['r']=14 # see issue 935
packing2_fix = ps.generators.spheres_from_coords(df_fix,mode='contained',smooth=True)

print(f'np.sum(packing): {np.sum(packing)}')
print(f'np.sum(packing2): {np.sum(packing2)}')
print(f'np.sum(packing2_fix): {np.sum(packing2_fix)}')