EDIT: pymrt.geometry
has been removed in favor of raster_geometry
.
DISCLAIMER: I am the author of both pymrt
and raster_geometry
.
If you just need to have the sphere, you can use the pip
-installable module raster_geometry
, and particularly raster_geometry.sphere()
, e.g:
import raster_geometry as rg
arr = rg.sphere(3, 1)
print(arr.astype(np.int_))
# [[[0 0 0]
# [0 1 0]
# [0 0 0]]
# [[0 1 0]
# [1 1 1]
# [0 1 0]]
# [[0 0 0]
# [0 1 0]
# [0 0 0]]]
internally, this is implemented as an n-dimensional superellipsoid generator, you can check its source code for details.
Briefly, the (simplified) code would reads like this:
import numpy as np
def sphere(shape, radius, position):
"""Generate an n-dimensional spherical mask."""
# assume shape and position have the same length and contain ints
# the units are pixels / voxels (px for short)
# radius is a int or float in px
assert len(position) == len(shape)
n = len(shape)
semisizes = (radius,) * len(shape)
# genereate the grid for the support points
# centered at the position indicated by position
grid = [slice(-x0, dim - x0) for x0, dim in zip(position, shape)]
position = np.ogrid[grid]
# calculate the distance of all points from `position` center
# scaled by the radius
arr = np.zeros(shape, dtype=float)
for x_i, semisize in zip(position, semisizes):
# this can be generalized for exponent != 2
# in which case `(x_i / semisize)`
# would become `np.abs(x_i / semisize)`
arr += (x_i / semisize) ** 2
# the inner part of the sphere will have distance below or equal to 1
return arr <= 1.0
and testing it:
# this will save a sphere in a boolean array
# the shape of the containing array is: (256, 256, 256)
# the position of the center is: (127, 127, 127)
# if you want is 0 and 1 just use .astype(int)
# for plotting it is likely that you want that
arr = sphere((256, 256, 256), 10, (127, 127, 127))
# just for fun you can check that the volume is matching what expected
# (the two numbers do not match exactly because of the discretization error)
print(np.sum(arr))
# 4169
print(4 / 3 * np.pi * 10 ** 3)
# 4188.790204786391
I am failing to get how your code exactly works, but to check that this is actually producing spheres (using your numbers) you could try:
arr = sphere((256, 256, 256), 10, (127, 127, 127))
# plot in 3D
import matplotlib.pyplot as plt
from skimage import measure
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
verts, faces, normals, values = measure.marching_cubes(arr, 0.5)
ax.plot_trisurf(
verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap='Spectral',
antialiased=False, linewidth=0.0)
plt.show()
Other approaches
One could implement essentially the same with a combination of np.linalg.norm()
and np.indices()
:
import numpy as np
def sphere_idx(shape, radius, position):
"""Generate an n-dimensional spherical mask."""
assert len(position) == len(shape)
n = len(shape)
position = np.array(position).reshape((-1,) + (1,) * n)
arr = np.linalg.norm(np.indices(shape) - position, axis=0)
return arr <= radius
producing the same results (sphere_ogrid
is sphere
from above):
import matplotlib.pyplot as plt
funcs = sphere_ogrid, sphere_idx
fig, axs = plt.subplots(1, len(funcs), squeeze=False, figsize=(4 * len(funcs), 4))
d = 500
n = 2
shape = (d,) * n
position = (d // 2,) * n
size = (d // 8)
base = sphere_ogrid(shape, size, position)
for i, func in enumerate(funcs):
arr = func(shape, size, position)
axs[0, i].imshow(arr)
However, this is going to be substantially slower and requires much more temporary memory n_dim * shape
of the output.
The benchmarks below seems to support the speed assessment:
base = sphere_ogrid(shape, size, position)
for func in funcs:
print(f"{func.__name__:20s}", np.allclose(base, arr), end=" ")
%timeit -o func(shape, size, position)
# sphere_ogrid True 1000 loops, best of 5: 866 µs per loop
# sphere_idx True 100 loops, best of 5: 4.15 ms per loop