CODEWITHSUNDEEP
pythonnumpy3dgeometrymasking
Soltius
Soltius

Reputation: 2261

Masking a 3D numpy array with a tilted disc

I'm working in Python2.7 with 3D numpy arrays, and trying to retrieve only pixels who fall on a 2D tilted disc.

Here is my code to plot the border of the disc (= a circle) I am interested in

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


#creating a 3d numpy array (empty in this example, but will represent a binary 3D image in my application)
space=np.zeros((40,40,20))

r = 8 #radius of the circle
theta = np.pi / 4 # "tilt" of the circle
phirange = np.linspace(0, 2 * np.pi) #to make a full circle

#center of the circle
center=[20,20,10]

#computing the values of the circle in spherical coordinates and converting them
#back to cartesian
for phi in phirange:
    x = r * np.cos(theta) * np.cos(phi) + center[0]
    y=  r*np.sin(phi)  + center[1]
    z=  r*np.sin(theta)* np.cos(phi) + center[2]
    space[int(round(x)),int(round(y)),int(round(z))]=1


x,y,z = space.nonzero()

#plotting
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, zdir='z', c= 'red')
plt.show()

The plot gives the following figure :

circle

which is a good start, but now I want a way to retrieve only the values of the pixels of space which are located in the disc defined by the circle : the ones in the pink zone in the following image (in my application, space will be a 3D binary image, here it is numpy.zeros() just to be able to plot and show you the disc I want):

disc

How should I procede ? I guess there is some numpy masking involved, an I understand how you would do it in 2D (like this question) but I'm having trouble applying this to 3D.

Upvotes: 0

Views: 1214

Answers (2)

Paul Panzer
Paul Panzer

Reputation: 53029

One easy way would be to calculate the normal vector to your disc plane. You can use your spherical coordinates for that. Be sure not to add the centre, set phi at zero and swap cos and sin theta, also stick a minus sign to the sin.

lets call that vector v. The plane is given by v0*x0 + v1*x1 + v2*x2 == c you can calculate c by inserting a point from your circle for x.

Next you can make a 2d grid for x0 and x1 and solve for x2. this gives you the height x2 as a function of the x0, x1 mesh. for these points you can calculate the distance from your disc centre and discard the points that are too far off. This you would indeed do using a mask.

Finally, depending on how precisely you want to plot you could round the x2 values to grid units, but for example for a surface plot I wouldn't do that.

To get a 3d mask as you describe you would round x2 and then starting from an all zero space set the disc pixels using space[x0, x1, x2] = True. This assumes that you have masked x0, x1, x2 as described earlier.

Upvotes: 1

Chazeon
Chazeon

Reputation: 546

Well that is a math problem, you should ask it in the Mathematics Stack Exchange site.

From my perspective, you should first find the surface your disc is in, and do the area calculation within that surface, by, for example, the method you mentioned in the linked question.

numpy or matplotlib here definitely do not responsible for the projection, you do.

Without clearly point out which (or which kind of) surface they are in, and the equation does not guarantee it is a plane, the area does not mean anything.

Upvotes: 1

Related Questions