给定一个形状为 (256, 256, 256) 的 3D numpy 数组，我如何在里面制作一个实心球体形状？下面的代码生成一系列递增和递减的圆圈，但在其他两个维度中查看时呈菱形。

def make_sphere(arr, x_pos, y_pos, z_pos, radius=10, size=256, plot=False):

    val = 255            
    for r in range(radius):
        y, x = np.ogrid[-x_pos:n-x_pos, -y_pos:size-y_pos]
        mask = x*x + y*y <= r*r 
        top_half = arr[z_pos+r]
        top_half[mask] = val #+ np.random.randint(val)
        arr[z_pos+r] = top_half

    for r in range(radius, 0, -1):
        y, x = np.ogrid[-x_pos:size-x_pos, -y_pos:size-y_pos]
        mask = x*x + y*y <= r*r 
        bottom_half = arr[z_pos+r]
        bottom_half[mask] = val#+ np.random.randint(val)
        arr[z_pos+2*radius-r] = bottom_half

    if plot:
        for i in range(2*radius):
            if arr[z_pos+i].max() != 0:
                print(z_pos+i)
                plt.imshow(arr[z_pos+i])
                plt.show()

    return arr

EDIT: pymrt.geometry已被删除以支持raster_geometry https://pypi.org/project/raster-geometry/.

免责声明: 我是这两本书的作者pymrt and raster_geometry.

如果您只需要球体，则可以使用pip-可安装模块raster_geometry，特别是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]]]

在内部，这是作为一个n维超椭球生成器实现的，你可以检查它源代码 https://github.com/norok2/raster_geometry/blob/e6e6d2fadca13bf3568da7fe29d59fb84d073934/raster_geometry/raster.py#L1156了解详情。 简而言之，（简化的）代码将如下所示：

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

并测试它：

# 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

我无法了解您的代码到底是如何工作的，但要检查这是否实际上产生了球体（使用您的数字），您可以尝试：

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()

其他方法

人们可以通过组合来实现本质上相同的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

产生相同的结果（sphere_ogrid is sphere从上面）：

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)

然而，这会慢得多并且需要更多的临时内存n_dim * shape的输出。 以下基准似乎支持速度评估：

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