2017-02-17 1 views
0

Je travaille en Python2.7 avec des tableaux 3D numpy, et j'essaie de ne récupérer que les pixels qui tombent sur un disque 2D incliné.Masquage d'un tableau numpy 3D avec un disque incliné

Voici mon code pour tracer la frontière du disque (= un cercle) Je suis intéressé par

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() 

L'intrigue donne la figure suivante:

circle

qui est un bon début, mais maintenant je veux un moyen de récupérer uniquement les valeurs des pixels de space qui se trouvent dans le disque défini par le cercle: ceux de la zone rose dans l'image suivante (dans mon application, space sera être une image binaire 3D, ici il est numpy.zeros() juste pour être en mesure de tracer et vous montrer le disque que je veux):

disc

Comment dois-je procede? Je suppose qu'il y a un certain nombre de masques impliqués, je comprends comment vous le feriez en 2D (comme this question) mais j'ai du mal à l'appliquer en 3D.

+0

Que voulez-vous dire par "tomber"? Si un pixel est un point, il n'y a aucune garantie qu'il coupe l'avion du tout, encore moins à l'intérieur d'un cercle. Voulez-vous plutôt regarder à l'intérieur d'un cylindre de hauteur 1? – Eric

+0

Oui, je suis conscient de cela, travailler avec l'approximation est bon pour cette application, d'où le round() dans le code que je montre: les points sur la figure sont approximativement sur un cercle. Je suppose que vous pourriez dire que cela équivaut à chercher des points dans un cylindre haute de 1pixel, oui. – Soltius

Répondre

1

Un moyen simple serait de calculer le vecteur normal à votre plan de disque. Vous pouvez utiliser vos coordonnées sphériques pour cela. Soyez sûr que pas pour ajouter le centre, mettre phi à zéro et échanger cos et sin theta, également coller un signe moins au péché.

Appelons ce vecteur v. Le plan est donné par v0 * x0 + v1 * x1 + v2 * x2 == c vous pouvez calculer c en insérant un point de votre cercle pour x.

Ensuite, vous pouvez faire une grille 2d pour x0 et x1 et résoudre pour x2. cela vous donne la hauteur x2 en fonction du maillage x0, x1. pour ces points, vous pouvez calculer la distance à partir du centre du disque et rejeter les points qui sont trop éloignés. C'est ce que vous feriez en utilisant un masque. Enfin, en fonction de la précision avec laquelle vous voulez tracer, vous pouvez arrondir les valeurs x2 aux unités de la grille, mais par exemple pour un tracé de surface, je ne ferais pas cela. Pour obtenir un masque 3D tel que vous le décrivez, vous devez arrondir x2 puis, à partir d'un espace entièrement nul, définir les pixels du disque en utilisant l'espace [x0, x1, x2] = True. Cela suppose que vous avez masqué x0, x1, x2 comme décrit précédemment.

+0

Merci pour votre contribution, qui m'a conduit dans la bonne direction. Je n'ai pas fait exactement ce que vous avez proposé, mais mon code est basé sur le calcul de vecteurs normaux et autres. Plus exactement, j'ai utilisé la technique de http://demonstrations.wolfram.com/ParametricEquationOfACircleIn3D/ pour obtenir le paramétrage du cercle extérieur. Puis je suis passé par un rayon de 0 à mon maxradius pour obtenir un disque! – Soltius

1

Eh bien c'est un problème math, vous devriez le demander dans le site Mathematics Stack Exchange. De mon point de vue, vous devriez d'abord trouver la surface dans laquelle se trouve votre disque, et faire le calcul de surface à l'intérieur de cette surface, par exemple, la méthode que vous avez mentionnée dans la question liée.

numpy ou matplotlib ici certainement pas responsable de la projection, vous faites. Sans indiquer clairement quelle (ou quelle sorte) de surface ils sont, et l'équation ne garantit pas qu'il s'agit d'un plan, la zone ne veut rien dire.