I have a ndarray of size 32x32x32. I want to create a sphere inside the array with the center at (x,y) and a radius of 4 pixels. The value of the sphere is 1 while value of the array is 0. How can this be done in python?
This is the code to generate the array:
import numpy as np
A = np.zeros((32,32,32))
print (A)
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.
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:
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.
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.
Nothing above worked for me so there is my attempt:
def create_bin_sphere(arr_size, center, r):
coords = np.ogrid[:arr_size[0], :arr_size[1], :arr_size[2]]
distance = np.sqrt((coords[0] - center[0])**2 + (coords[1]-center[1])**2 + (coords[2]-center[2])**2)
return 1*(distance <= r)
where:
Example:
arr_size = (30,30,30)
sphere_center = (15,15,15)
r=10
sphere = create_bin_sphere(arr_size,sphere_center, r)
#Plot the result
fig =plt.figure(figsize=(6,6))
ax = fig.gca(projection='3d')
ax.voxels(sphere, facecolors=colors, edgecolor='k')
plt.show()
Eventhough it's a little late - I recently faced the same problem and solved it somewhat like the solution supposed by Mstaino. Also, Mstainos solution doesn't work for asymetricly sized array's, because shapes will not match when calculating the distance.
So here's my approach in 3D which produces a sphere in the center of the array. :
# define array size and sphere radius
size = [size_x, size_y, size_z]
radius = sphere_radius
# compute center index of array
center = [int(size[0]/2), int(size[1]/2), int(size[2]/2)]
# create index grid for array
ind_0, ind_1, ind_2 = np.indices((size[0], size[1], size[2]))
# calculate "distance" of indices to center index
distance = ((ind_0 - center[0])**2 + (ind_1 - center[1])**2 + (ind_2 - center[2])**2)**.5
# create output
output = np.ones(shape = (size[0], size[1], size[2])) * (distance <= radius)
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.]]]