如何用Python的NumPy获取3D图像中球体内的所有像素

核心思路是通过计算每个像素到球心的欧几里得距离，筛选出距离小于等于半径的像素。为了提升计算效率，我们可以直接比较距离的平方（省去开根号的浮点运算，结果判断逻辑完全一致）。

具体步骤如下：

1. 计算各像素坐标与球心的差值
   利用NumPy的广播特性，直接对整个坐标数组做减法：

   dx = indicesX - centerX
   dy = indicesY - centerY
   dz = indicesZ - centerZ
   

2. 计算距离的平方
   基于三维欧几里得距离公式，计算每个像素到球心的距离平方：

   distance_squared = dx ** 2 + dy ** 2 + dz ** 2
   

3. 生成筛选掩码
   创建布尔数组，标记出距离平方小于等于半径平方的像素：

   radius_squared = r ** 2
   mask = distance_squared <= radius_squared
   

4. 提取球体内的像素
   用掩码直接索引原坐标数组，得到球体内的所有像素坐标；如果需要提取对应图像的像素值，同样用掩码索引图像数组即可：

   # 获取球体内的像素坐标
   inside_indices_x = indicesX[mask]
   inside_indices_y = indicesY[mask]
   inside_indices_z = indicesZ[mask]
   
   # 假设你的3D图像数组为img，提取球体内的像素值
   # inside_pixels = img[mask]
   

完整示例 #

如果你的indicesX、indicesY、indicesZ是通过np.mgrid生成的（这是3D图像坐标常见的生成方式），完整代码示例如下：

import numpy as np

# 假设3D图像尺寸为(64, 64, 64)
img_shape = (64, 64, 64)
indicesZ, indicesY, indicesX = np.mgrid[0:img_shape[0], 0:img_shape[1], 0:img_shape[2]]

# 球心坐标和半径
centerX, centerY, centerZ = 32, 32, 32
r = 10

# 计算差值
dx = indicesX - centerX
dy = indicesY - centerY
dz = indicesZ - centerZ

# 计算距离平方并生成掩码
distance_squared = dx**2 + dy**2 + dz**2
mask = distance_squared <= r**2

# 提取球体内的像素坐标和值
inside_x = indicesX[mask]
inside_y = indicesY[mask]
inside_z = indicesZ[mask]
# inside_pixels = img[mask]  # 如果有图像数组的话

这种方法完全利用NumPy的向量化运算，避免了循环遍历，效率远高于逐像素判断。

内容的提问来源于stack exchange，提问作者Mary