0

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.

Community
  • 1
  • 1
Soltius
  • 2,162
  • 1
  • 16
  • 28
  • What do you mean by "fall on"? If a pixel is a point, there's no guarantee that it intersects the plane at all, let alone inside a circle. Do you instead want to look inside a cylinder of height 1? – Eric Feb 17 '17 at 14:42
  • Yes, I'm aware of that, working with approximation is fine for this application, hence the round() in the code I show : the points on the figure are approximatively on a circle. I guess you could say this is equivalent to looking for points in a 1pixel high cylinder, yes. – Soltius Feb 20 '17 at 08:21

2 Answers2

1

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.

Community
  • 1
  • 1
Chazeon
  • 546
  • 2
  • 14
1

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.

Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
  • Thanks for your input, which lead me into the right direction. I didn't do exactly what you proposed, but my code is based on computing normal vectors and such. More exactly, I used the technique in http://demonstrations.wolfram.com/ParametricEquationOfACircleIn3D/ to get the parametrisation of the outer circle. Then I just went through a radius from 0 to my maxradius to get a disc ! – Soltius Feb 22 '17 at 15:04