Wie kann man die Orientierung von 3D-Normalen in großen Punktwolken effizient berechnen?

Question

Ich arbeite an einer Python-Bibliothek für die Verwaltung von Punktwolken in Python (mehr wie "learning by doing").

Ich habe eine Funktion geschrieben, um die Ausrichtung jeder Normalen in einer Punktwolke zu berechnen, die als ein numpy strukturiertes Array gespeichert ist, aber ich bin nicht glücklich genug mit der endgültigen Funktion (dachte, es funktioniert und ziemlich schnell genug) und ich war frage mich, ob es noch einen effizienteren/pythonschen Ansatz zur Berechnung der Orientierung in großen Punktwolken gibt.

Dies ist, wie die Punktwolke strukturiert:

esfera = PyntCloud.from_ply('Sphere.ply') 

esfera.vertex 
Out[3]: 
array([ (0.2515081465244293, 0.05602749437093735, 1.9830318689346313, 0.12660565972328186, 0.02801010198891163, 0.9915575981140137, 7.450349807739258, 77.52488708496094), 
     (0.09723527729511261, 0.02066999115049839, 1.9934484958648682, 0.048643846064805984, 0.011384730227291584, 0.9987513422966003, 2.863548517227173, 76.82744598388672), 
     (0.17640848457813263, 0.028193067759275436, 1.9881943464279175, 0.08916780352592468, 0.01611466333270073, 0.9958862066268921, 5.198856830596924, 79.75591278076172), 
     ..., 
     (0.17817874252796173, -0.046098098158836365, -1.9879237413406372, 0.08992616087198257, -0.02275240235030651, -0.9956884980201721, 5.322407245635986, 284.19854736328125), 
     (0.2002459168434143, -0.002330917865037918, -1.986855149269104, 0.09960971027612686, -0.0010710721835494041, -0.9950260519981384, 5.717002868652344, 270.6160583496094), 
     (0.12885123491287231, -0.03245270624756813, -1.9912745952606201, 0.06637085974216461, -0.01580258458852768, -0.9976698756217957, 3.912114381790161, 283.3924865722656)], 
     dtype=[('x', '<f4'), ('y', '<f4'), ('z', '<f4'), ('nx', '<f4'), ('ny', '<f4'), ('nz', '<f4'), ('scalar_Dip_(degrees)', '<f4'), ('scalar_Dip_direction_(degrees)', '<f4')]) 

esfera.vertex['nx'] 
Out[4]: 
array([ 0.12660566, 0.04864385, 0.0891678 , ..., 0.08992616, 
     0.09960971, 0.06637086], dtype=float32) 

esfera.vertex[-1]['nx'] 
Out[5]: 0.06637086 

Und dies ist die Orientierungsfunktion:

def add_orientation(self, degrees=True): 

    """ Adds orientation (with respect to y-axis) values to PyntCloud.vertex 

    This function expects the PyntCloud to have a numpy structured array 
    with normals x,y,z values (correctly named) as the corresponding vertex 
    atribute. 

    Args: 
     degrees (Optional[bool]): Set the oputput orientation units. 
      If True(Default) set units to degrees. 
      If False set units to radians. 
    """ 

    #: set copy to False for efficience in large pointclouds 
    nx = self.vertex['nx'].astype(np.float64, copy=False) 
    ny = self.vertex['ny'].astype(np.float64, copy=False) 

    #: get orientations 
    angle = np.arctan(np.absolute(nx/ny)) 

    #: mask for every quadrant 
    q2 = np.logical_and((self.vertex['nx']>0),(self.vertex['ny']<0)) 
    q3 = np.logical_and((self.vertex['nx']<0),(self.vertex['ny']<0)) 
    q4 = np.logical_and((self.vertex['nx']<0),(self.vertex['ny']>0)) 

    #: apply modification for every quadrant 
    angle[q2] = np.pi - angle[q2] 
    angle[q3] = np.pi + angle[q3] 
    angle[q4] = (2*np.pi) - angle[q4] 

    if degrees == False: 
     orientation = np.array(angle, dtype=[('orir', 'f4')]) 
    else: 
     orientation = np.array((180 * angle/np.pi), dtype=[('orid', 'f4')]) 

    #: merge the structured arrays and replace the old vertex attribute 
    self.vertex = join_struct_arrays([self.vertex, orientation]) 

Und die in CloudCompare visualisierte Ergebnisse (nicht enoght rep hat für die Zeit nach Bildern.):

https://raw.githubusercontent.com/daavoo/sa/master/Captura%20de%20pantalla%20de%202016-03-21%2013%3A28%3A39.png

Vielen Dank für Ihre hel p.

Answer 1

Nun, ich schäme mich für mich. xD

Diese numpy eingebaute Funktion war genau das, was ich suchte.

Dank @Dan.

Hier ist die neue Funktion:

def add_orientation(self, degrees=True): 

     """ Adds orientation (with respect to y-axis) values to PyntCloud.vertex 

     This function expects the PyntCloud to have a numpy structured array 
     with normals x,y,z values (correctly named) as the corresponding vertex 
     atribute. 

     Args: 
      degrees (Optional[bool]): Set the oputput orientation units. 
       If True(Default) set units to degrees. 
       If False set units to radians. 
     """ 

     #: set copy to False for efficience in large pointclouds 
     nx = self.vertex['nx'].astype(np.float64, copy=False) 
     ny = self.vertex['ny'].astype(np.float64, copy=False) 

     #: get orientations 
     angle = np.arctan2(nx,ny) 

     #: convert (-180 , 180) to (0 , 360) 
     angle[(np.where(angle < 0))] = (2*np.pi) + angle[(np.where(angle < 0))] 

     if degrees: 
      orientation = np.array(np.rad2deg(angle), dtype=[("orid2",'f4')]) 
     else: 
      orientation = np.array(angle, dtype=[("orir2",'f4')]) 

     self.vertex = join_struct_arrays([self.vertex, orientation]) 

Wich ist einfacher und schneller.

t0 = t.time() 
esfera.add_orientation() 
t1 = t.time() 
dif = t1-t0 
dif 
Out[18]: 0.34514379501342773 

t0 = t.time() 
esfera.add_orientation2() 
t1 = t.time() 
dif = t1-t0 
dif 
Out[20]: 0.291456937789917 

Jetzt bin ich so glücklich wie beschämt.

Nächstes Mal werde ich einen tieferen Blick auf die numpy Docs werfen, bevor Sie eine Frage stellen.

Danke.

comp = esfera.vertex['orid'] == esfera.vertex['orid2'] 

np.all(comp) 
Out[15]: True