PMEAL / porespy

A set of tools for characterizing and analying 3D images of porous materials
https://porespy.org
MIT License
290 stars 97 forks source link

pseudo_gravity_packing: for clearance=0 and smooth=True: empty voxel between spheres in some siuations #957

Open heinsimon opened 2 months ago

heinsimon commented 2 months ago

Situation

With 58656a5 in dev: The generator pseudo_gravity_packing with smooth enabled and clearance=0 will leave one voxel distance between spheres in some situations. The spheres touch with smooth=False

Expected

I would expect no voxel distance between spheres. Hence, my question: Is this a bug?

Example


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

plt.close('all')

im_box=np.full((65,31,31),False) # shows bug
im_box=np.full((65,34,34),False) # shows bug
# im_box=np.full((65,39,38),False) # does not show bug

ps.settings.loglevel = "DEBUG"
#%%
numspheres=2;

packing_SmoothOff=ps.generators.pseudo_gravity_packing(im=im_box,r=15,clearance=0,maxiter=numspheres,seed=42,smooth=False)

# Bug in pseudo_gravity_packing: Spheres with smooth=True do not touch
packing_SmoothOn=ps.generators.pseudo_gravity_packing(im=im_box,r=15,clearance=0,maxiter=numspheres,seed=42,smooth=True)

#%%
fig,ax=plt.subplots(2,3, figsize=[12, 10])
ax[0,0].imshow(np.sum(packing_SmoothOn,axis=1))
ax[0,0].set_title('packing SmoothOn')
ax[0,1].imshow(np.sum(packing_SmoothOn,axis=1)>0)
ax[0,1].set_title('packing SmoothOn sum axis=1 >0')
ax[0,2].plot(np.int16(np.sum(np.sum(packing_SmoothOn,axis=2),axis=1)==0)*500)
ax[0,2].plot(np.sum(np.sum(packing_SmoothOn,axis=2),axis=1))
ax[0,2].set_title('packing SmoothOn sum axis1 and 2')

ax[1,0].imshow(np.sum(packing_SmoothOff,axis=1))
ax[1,0].set_title('packing SmoothOff')
ax[1,1].imshow(np.sum(packing_SmoothOff,axis=1)>0)
ax[1,1].set_title('packing SmoothOff sum axis=1 >0')
ax[1,2].plot(np.int16(np.sum(np.sum(packing_SmoothOff,axis=2),axis=1)==0)*500)
ax[1,2].plot(np.sum(np.sum(packing_SmoothOff,axis=2),axis=1))
ax[1,2].set_title('packing SmoothOff sum axis1 and 2')

#%%
padded_packing_SmoothOff=np.pad(packing_SmoothOff,pad_width=((0,0),(1,1),(1,1)), mode='constant',constant_values=0)
padded_packing_SmoothOn=np.pad(packing_SmoothOn,pad_width=((0,0),(1,1),(1,1)), mode='constant',constant_values=0)

discon_voxels_SmoothOff = ps.filters.find_disconnected_voxels(im=padded_packing_SmoothOff,conn=6)
discon_voxels_SmoothOn = ps.filters.find_disconnected_voxels(im=padded_packing_SmoothOn,conn=6)
print(f'Unique off: {np.unique(discon_voxels_SmoothOff)}')
print(f'Unique on: {np.unique(discon_voxels_SmoothOn)}')# will show [False True] if bug is present

discon_voxels_SmoothOn_unpad=ps.tools.unpad(discon_voxels_SmoothOn,((0,0),(1,1),(1,1)))

fig,ax=plt.subplots(1,2, figsize=[12, 10])
ax[0].imshow(np.sum(packing_SmoothOn,axis=1))
ax[0].set_title('packing SmoothOn')
ax[1].imshow(np.sum(discon_voxels_SmoothOn_unpad,axis=1))
ax[1].set_title('packing SmoothOn: Disconnected voxels')

Figures for image size (65,34,34)

Comparision smooth on and off

Figure_1_65_34_34_comparision_smoothOn_Off

unconnected voxels for smooth on

Figure_2_65_34_34_smoothOn_unconnected

jgostick commented 2 months ago

It might be a bug. The sphere touches the edge of the image image, which is correct, so perhaps I need to check the logic for the clearance. Thanks for reporting this.