如何在3D Numpy數組中生成一個球體

Question

給定一個3D numpy數組形狀（256,256,256），我將如何在內部形成一個堅實的球體形狀？下面的代碼會生成一系列增加和減少的圓，但是在另外兩個維度中查看時是菱形。

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 

Answer 1

免責聲明：我的pymrt作者。

如果你只是需要有球，你可以使用pip -installable模塊pymrt，特別pymrt.geometry.sphere()，如：

import pymrt as mrt 
import pymrt.geometry 

arr = mrt.geometry.sphere(3, 1) 

array([[[False, False, False], 
     [False, True, False], 
     [False, False, False]], 

     [[False, True, False], 
     [ True, True, True], 
     [False, True, False]], 

     [[False, False, False], 
     [False, True, False], 
     [False, False, False]]], dtype=bool) 

內部，這被實現爲n維superellipsoid發電機，你可以檢查其source code的細節。 簡單地說，（簡化）代碼會讀這樣的：

import numpy 

def sphere(shape, radius, position): 
    # assume shape and position are both a 3-tuple of int or float 
    # the units are pixels/voxels (px for short) 
    # radius is a int or float in px 
    semisizes = (radius,) * 3 

    # 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): 
     arr += (np.abs(x_i/semisize) ** 2) 
    # the inner part of the sphere will have distance below 1 
    return arr <= 1.0 

arr = sphere((256, 256, 256), 10, (127, 127, 127))) 
# 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 

# just for fun you can check that the volume is matching what expected 
np.sum(arr) 
# gives: 4169 

4/3 * np.pi * 10 ** 3 
# gives: 4188.790204786391 
# (the two numbers do not match exactly because of the discretization error) 

我沒有得到如何你的代碼完全相同的工作，但要檢查，這是實際生產領域（使用你的號碼），你可以嘗試：

import pymrt as mrt 
import pymrt.geometry 

arr = mrt.geometry.sphere(256, 10, 0.5) 


# 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, (2,) * 3) 
ax.plot_trisurf(
    verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap='Spectral', 
    antialiased=False, linewidth=0.0) 
plt.show()