How to generate a sphere in 3D Numpy array

Question

Given a 3D numpy array of shape (256, 256, 256), how would I make a solid sphere shape inside? The code below generates a series of increasing and decreasing circles but is diamond shaped when viewed in the two other dimensions.

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

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

Answer 2

size = 100
radius = 10

x0, y0, z0 = (50, 50, 50)

x, y, z = np.mgrid[0:size:1, 0:size:1, 0:size:1]
r = np.sqrt((x - x0)**2 + (y - y0)**2 + (z - z0)**2)
r[r > radius] = 0

Answer 3

Nice question. My answer to a similar question would be applicable here also.

You can try the following code. In the below mentioned code AA is the matrix that you want.

import numpy as np
from copy import deepcopy

''' size : size of original 3D numpy matrix A.
    radius : radius of circle inside A which will be filled with ones. 
'''
size, radius = 5, 2

''' A : numpy.ndarray of shape size*size*size. '''
A = np.zeros((size,size, size)) 

''' AA : copy of A (you don't want the original copy of A to be overwritten.) '''
AA = deepcopy(A) 

''' (x0, y0, z0) : coordinates of center of circle inside A. '''
x0, y0, z0 = int(np.floor(A.shape[0]/2)), \
        int(np.floor(A.shape[1]/2)), int(np.floor(A.shape[2]/2))


for x in range(x0-radius, x0+radius+1):
    for y in range(y0-radius, y0+radius+1):
        for z in range(z0-radius, z0+radius+1):
            ''' deb: measures how far a coordinate in A is far from the center. 
                    deb>=0: inside the sphere.
                    deb<0: outside the sphere.'''   
            deb = radius - abs(x0-x) - abs(y0-y) - abs(z0-z) 
            if (deb)>=0: AA[x,y,z] = 1

Following is an example of the output for size=5 and radius=2 (a sphere of radius 2 pixels inside a numpy array of shape 5*5*5):

[[[0. 0. 0. 0. 0.]
  [0. 0. 0. 0. 0.]
  [0. 0. 1. 0. 0.]
  [0. 0. 0. 0. 0.]
  [0. 0. 0. 0. 0.]]

 [[0. 0. 0. 0. 0.]
  [0. 0. 1. 0. 0.]
  [0. 1. 1. 1. 0.]
  [0. 0. 1. 0. 0.]
  [0. 0. 0. 0. 0.]]

 [[0. 0. 1. 0. 0.]
  [0. 1. 1. 1. 0.]
  [1. 1. 1. 1. 1.]
  [0. 1. 1. 1. 0.]
  [0. 0. 1. 0. 0.]]

 [[0. 0. 0. 0. 0.]
  [0. 0. 1. 0. 0.]
  [0. 1. 1. 1. 0.]
  [0. 0. 1. 0. 0.]
  [0. 0. 0. 0. 0.]]

 [[0. 0. 0. 0. 0.]
  [0. 0. 0. 0. 0.]
  [0. 0. 1. 0. 0.]
  [0. 0. 0. 0. 0.]
  [0. 0. 0. 0. 0.]]]

I haven't printed the output for the size and radius that you had asked for (size=32 and radius=4), as the output will be very long.

Answer 4

Here is how to create voxels space without numpy, the main idea that you calculate distance between center and voxel and if voxel in radius you will create.

from math import sqrt    

def distance_dimension(xyz0 = [], xyz1 = []):
    delta_OX = pow(xyz0[0] - xyz1[0], 2)
    delta_OY = pow(xyz0[1] - xyz1[1], 2)
    delta_OZ = pow(xyz0[2] - xyz1[2], 2)
    return sqrt(delta_OX+delta_OY+delta_OZ)    



def voxels_figure(figure = 'sphere', position = [0,0,0], size = 1):
    xmin, xmax = position[0]-size,  position[0]+size
    ymin, ymax = position[1]-size,  position[1]+size
    zmin, zmax = position[2]-size,  position[2]+size

    voxels = []

    if figure == 'cube':
        for local_z, world_z in zip(range(zmax-zmin), range(zmin, zmax)):
            for local_y, world_y in zip(range(ymax-ymin), range(ymin, ymax)):
                for local_x, world_x in zip(range(xmax-xmin), range(xmin, xmax)):
                    voxels.append([world_x,world_y,world_z])

    elif figure == 'sphere':
        for local_z, world_z in zip(range(zmax-zmin), range(zmin, zmax)):
            for local_y, world_y in zip(range(ymax-ymin), range(ymin, ymax)):
                for local_x, world_x in zip(range(xmax-xmin), range(xmin, xmax)):
                    radius = distance_dimension(xyz0 = [world_x, world_y,world_z], xyz1 = position)
                    if  radius < size:
                        voxels.append([world_x,world_y,world_z])

    return voxels

voxels = voxels_figure(figure = 'sphere', position = [0,0,0], size = 3)

After you will get voxels indexes, you can apply ~ones for cube matrix.

Answer 5

Instead of using loops, I propose to use a meshgrid + sphere equation + np.where

import numpy as np 
def generate_sphere(volumeSize):
    x_ = np.linspace(0,volumeSize, volumeSize)
    y_ = np.linspace(0,volumeSize, volumeSize)
    z_ = np.linspace(0,volumeSize, volumeSize)
    r = int(volumeSize/2) # radius can be changed by changing r value
    center = int(volumeSize/2) # center can be changed here
    u,v,w = np.meshgrid(x_, y_, z_, indexing='ij')
    a = np.power(u-center, 2)+np.power(v-center, 2)+np.power(w-center, 2)
    b = np.where(a<=r*r,1,0)
    return b