points_to_spheres not adding spheres correctly for 3D images

[heinsimon]

The function points_to_spheres is not correctly adding spheres for a 3D image.

radius=10
im=np.full((41,41,41),0)
im[20,20,20]=radius

res_ps=ps.tools.points_to_spheres(im=im)
res_fixed = points_to_spheres_fixed(im=im)
print(f'Number voxels porespy: {np.sum(res_ps)}') # Number voxels porespy: 1
print(f'Number voxels fixed: {np.sum(res_fixed)}') # Number voxels fixed: 4169
print(f'Volume analytic: {4/3*radius**3 *np.pi}') # Volume analytic: 4188.790204786391

fig, ax = plt.subplots(3, 2, figsize=[12, 6])
ax[0,0].imshow(ps.visualization.sem(~res_ps, axis=0), interpolation='none', origin='lower')
ax[0,0].axis(False)
ax[1,0].imshow(ps.visualization.sem(~res_ps, axis=2).T, interpolation='none', origin='lower')
ax[1,0].axis(False);
ax[2,0].imshow(ps.visualization.sem(~res_ps, axis=1).T, interpolation='none', origin='lower')
ax[2,0].axis(False);

ax[0,1].imshow(ps.visualization.sem(~res_fixed, axis=0), interpolation='none', origin='lower')
ax[0,1].axis(False)
ax[1,1].imshow(ps.visualization.sem(~res_fixed, axis=2).T, interpolation='none', origin='lower')
ax[1,1].axis(False);
ax[2,1].imshow(ps.visualization.sem(~res_fixed, axis=1).T, interpolation='none', origin='lower')
ax[2,1].axis(False);

Potential fix:

• remove '.T' in _sphere_insertions.py#L38
• Add check for 'im.ndim==3' and 'r = im[x, y, z].flatten()' in _sphere_insertions.py#L48

Fix used for example above

def points_to_spheres_fixed(im):
    r"""
    Inserts disks/spheres into an image at locations indicated by non-zero values

    Parameters
    ----------
    im : ndarray
        The image containing nonzeros indicating the locations to insert spheres.
        If the non-zero values are `bool`, then the maximal size is found and used;
        if the non-zeros are `int` then these values are used as the radii.

    Returns
    -------
    spheres : ndarray
        A `bool` array with disks/spheres inserted at each nonzero location in `im`.
    """
    from scipy.spatial import distance_matrix
    if im.ndim == 3:
        x, y, z = np.where(im > 0)
        coords = np.vstack((x, y, z))
    else:
        x, y = np.where(im > 0)
        coords = np.vstack((x, y))
    if im.dtype == bool:
        dmap = distance_matrix(coords.T, coords.T)
        mask = dmap < 1
        dmap[mask] = np.inf
        r = np.around(dmap.min(axis=0)/2, decimals=0).astype(int)
    else:
        if im.ndim == 3:
            r = im[x, y, z].flatten()
        else:
            r = im[x, y].flatten()
    im_spheres = np.zeros_like(im, dtype=bool)
    im_spheres = ps.tools._insert_disks_at_points_parallel(
        im_spheres,
        coords=coords,
        radii=r,
        v=True,
        smooth=False,
    )
    return im_spheres
    ```

[jgostick]

Thanks for the detailed bug report and solution. Much appreciated. Would you be willing to do a pull request? Those changes are pretty minor so it should be a quick and easy.