I have a ndarray size of 32x32x32. I want to create a sphere inside the array with center at (x,y) and radius 4 pixels. The value of sphere is 1 while value of array is 0. How could I do it in python? Thanks so much

This is code to generate the array

import numpy as np
A = np.zeros((32,32,32))
print (A)

Since array indexes only have a certain level of specificity (i.e. you can only subdivide down to the width, in this case 32), there's no one perfect way to represent a sphere in an array. Instead, we can treat each array index as a space of cubic area, where the [x][y][z] indices of the index represent the the cubic area's center coordinates. To create the sphere, we evaluate whether the sphere's presence in that area of space meets certain criteria.

We start with the equation for a sphere. From Wikipedia:

In analytic geometry, a sphere with center (x0, y0, z0) and radius r is the locus of all points (x, y, z) such that

(x - x0)^2 + (y - y0)^2 + (z - z0)^2 <= r^2.

For an array of dimensions N, the center will have the coordinate (N - 1) / 2 for all dimensions. (because for an even-numbered dimension, the center should be between the two center points, and for an odd-numbered dimension, the center should be an integer.) The magnitude of the radius can vary depending on where you decide the boundaries of the sphere relative to our imagined cubic array representation; re-reading the question, I notice you already gave the desired radius: 4.

There are two evaluation criteria I can think of:

Simple approach

In this approach, we will simply use a test of whether the array index's cubic area's center lies within the circle equation.

You can see Siddharth Satpathy's answer for some code using this approach.

Sophisticated approach

Ideally for me, the equation would decide whether an index lies within the sphere by assessing whether the proportion of sphere for that cubic area is greater than 50%. However, this approach unfortunately goes beyond my current working mathematical knowledge.

In regards to a discussion I had in the comments, neither approach is better than the other since they represent different perspectives: I personally imagine the array as being actually representative of the cubic area for each index, while others may imagine the indexes being the center points of these cubic areas.

Very good question. 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.

Using a combination of indexing, distance calculation and masking (all with numpy):

import numpy as np
center = (31/2, 31/2, 31/2) # if it is centered
size = (32, 32, 32)
max_dist = 4
distance = np.linalg.norm(np.subtract(np.indices(size).T, np.asarray(center)), axis=2)
#print(distance)
mask = np.ones(size) * (distance < max_dist)
print(mask)

np.indices creates an index in the form [[[(i, j, k)]]], np.substract calculates the vector difference to your center, and np.linalg.norm calculates the vector norm. The rest is just using a mask operation on the distance array.

Does that work?

EDIT: an example with (3,3,3) for clarity purposes

center = (1, 1, 1)
size = (3, 3, 3)
distance = np.linalg.norm(np.subtract(np.indices(size).T,np.asarray(center)), axis=len(center))
mask = np.ones(size) * (distance<=1)
print(mask)

>>[[[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.]]]