Drawing/sampling a sphere in a 3D numpy grid

Question

I want to do a voxel-based measurement of spherical objects, represented in a numpy array. Because of the sampling, these spheres are represented as a group of cubes (because they are sampled in the array). I want to do a simulation of the introduced error by this grid-restriction. Is there any way to paint a 3D sphere in a numpy grid to run my simulations on? (So basically, a sphere of unit length one, would be one point in the array)

Or is there another way of calculating the error introduced by sampling?

In 2-D that seems to be easy...

Answer 1

The most direct approach is to create a bounding box array, holding at each point the distance to the center of the sphere:

>>> radius = 3
>>> r2 = np.arange(-radius, radius+1)**2
>>> dist2 = r2[:, None, None] + r2[:, None] + r2
>>> volume = np.sum(dist2 <= radius**2)
>>> volume
123

The 2D case is easier to visualize:

>>> dist2 = r2[:, None] + r2
>>> (dist2 <= radius**2).astype(np.int)
array([[0, 0, 0, 1, 0, 0, 0],
       [0, 1, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 0],
       [1, 1, 1, 1, 1, 1, 1],
       [0, 1, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 0],
       [0, 0, 0, 1, 0, 0, 0]])
>>> np.sum(dist2 <= radius**2)
29