作者热门文章
- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
给定一个形状为 (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
最佳答案
编辑:pymrt.geometry
已被删除以支持 raster_geometry
.
免责声明:我是pymrt
和raster_geometry
的作者。
如果您只需要球体,您可以使用 pip
可安装模块 raster_geometry
,尤其是 raster_geometry.sphere()
,例如:
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 维超椭圆体生成器实现的,你可以查看它的 source code了解详情。简而言之,(简化的)代码如下所示:
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()
和 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
是上面的 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
关于python - 如何在 3D Numpy 数组中生成球体,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/46626267/
我是一名优秀的程序员,十分优秀!