2016-09-15 7 views
0

Ich habe eine Linie und einige Punkte, die auf dieser Linie im 3D-Raum sind. Ich weiß, dass es einen gewissen Fehler im Punkt gibt, aber der Fehler erstreckt sich nur senkrecht zur Linie. Um dies zu sehen, möchte ich Platten mit einem Radius des Fehlers, die auf der Linie und orthogonal zur Richtung der Linie zentriert sind. Ich habe dieses solution gefunden, aber ich kann es nicht zur Arbeit bringen.Matplotlib Rotierende 3D-Disk

Wenn ich den Code und den Zustand ich möchte die Ausgabe normal auf der 'Z' Achse Ich bekomme, was ich erwarten würde. Platten mit einem bestimmten Radius auf der Linie und ausgerichtet auf der Z-Achse.

pathpatch_2d_to_3d(p, z=z,normal='z') 

Image with normal='z'

Ich möchte die gedreht Scheiben. Um den Well-Vektor an diesem Punkt zu finden, benutze ich einen Punkt, indem ich diesen Vektor verwende. Dies ist der Vektor, den ich als normal=(vx,vy,vz) setze, aber wenn ich das tue, sind die Festplatten nicht einmal auf dem Diagramm. Ich bin mir nicht sicher, wo ich falsch liege. Hat jemand einen Rat?

Hier ist mein Code.

import matplotlib.pyplot as plt 
from matplotlib.patches import Circle, PathPatch 
from mpl_toolkits.mplot3d import Axes3D 
import mpl_toolkits.mplot3d.art3d as art3d 
import numpy as np 
from scipy.interpolate import interp1d 

md,wellz,wellx,welly=np.genfromtxt("./well.csv",delimiter=",",unpack=True) 

# Building interpolation function that map a measured depth to its repsective x,y,z coordinates 
fz = interp1d(md,wellz) 
fx = interp1d(md,wellx) 
fy = interp1d(md,welly) 

pointDepth = np.array([15790,15554,15215,14911,14274,13927,13625,13284,12983,12640,12345,12004,11704,11361,11061,10717,10418,10080,9771]) 


def rotation_matrix(d): 
    """ 
Calculates a rotation matrix given a vector d. The direction of d 
corresponds to the rotation axis. The length of d corresponds to 
the sin of the angle of rotation. 

Variant of: http://mail.scipy.org/pipermail/numpy-discussion/2009-March/040806.html 
    """ 
    sin_angle = np.linalg.norm(d) 

    if sin_angle == 0: 
    return np.identity(3) 

    d = d/sin_angle 

    eye = np.eye(3) 
    ddt = np.outer(d, d) 
    skew = np.array([[ 0, d[2], -d[1]], 
        [-d[2],  0, d[0]], 
        [d[1], -d[0], 0]], dtype=np.float64) 

    M = ddt + np.sqrt(1 - sin_angle**2) * (eye - ddt) + sin_angle * skew 
    return M 

def pathpatch_2d_to_3d(pathpatch, z = 0, normal = 'z'): 
    """ 
    Transforms a 2D Patch to a 3D patch using the given normal vector. 

    The patch is projected into they XY plane, rotated about the origin 
    and finally translated by z. 
    """ 
    if type(normal) is str: #Translate strings to normal vectors 
     index = "xyz".index(normal) 
     normal = np.roll((1,0,0), index) 

    normal = normal/np.linalg.norm(normal) #Make sure the vector is normalised 

    path = pathpatch.get_path() #Get the path and the associated transform 
    trans = pathpatch.get_patch_transform() 

    path = trans.transform_path(path) #Apply the transform 
    pathpatch.__class__ = art3d.PathPatch3D #Change the class 
    pathpatch._code3d = path.codes #Copy the codes 
    pathpatch._facecolor3d = pathpatch.get_facecolor #Get the face color  

    verts = path.vertices #Get the vertices in 2D 

    d = np.cross(normal, (0, 0, 1)) #Obtain the rotation vector 
    M = rotation_matrix(d) #Get the rotation matrix 
    pathpatch._segment3d = np.array([np.dot(M, (x, y, 0)) + (0, 0, z) for x, y in verts]) 

def pathpatch_translate(pathpatch, delta): 
    """ 
    Translates the 3D pathpatch by the amount delta. 
    """ 
    pathpatch._segment3d += delta 


fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d') 
ax.plot(wellx,welly,wellz,c='k') 

for n,pd in enumerate(pointDepth): 
    x,y,z = fx(pd),fy(pd),fz(pd) 

    # figure out a vector from the point 
    vx,vy,vz = (fx(pd-10)-x),(fy(pd-10)-y),(fz(pd-10)-z) 

    #Draw Circle 
    p = Circle((x,y), 100) 
    ax.add_patch(p) 
    pathpatch_2d_to_3d(p, z=z,normal=(vx,vy,vz)) 
    pathpatch_translate(p,(0,0,0)) 

ax.set_xlim3d(np.min(wellx),np.max(wellx)) 
ax.set_ylim3d(np.min(welly), np.max(welly)) 
ax.set_zlim3d(np.min(wellz), np.max(wellz)) 
plt.show() 

Antwort

0

Hier ist die Lösung, die ich mir ausgedacht habe. Ich entschied mich für den Unterschied zwischen dem Punkt, wo der Punkt auf die Linie fällt und dem ersten Punkt in p.._segment3d. Dies gibt mir, wie weit mein Kreis von dem Ort entfernt ist, an dem ich es haben möchte, dann übersetzte ich einfach den Ausschnitt, minus den Radius des Kreises, so dass er zentriert wäre.

Ich habe in einigen Zufallszahlen hinzugefügt, wie einige „Fehler“ und hier ist Endcode und resultierende Bild

import matplotlib.pyplot as plt 
from matplotlib.patches import Circle, PathPatch 
from mpl_toolkits.mplot3d import Axes3D 
import mpl_toolkits.mplot3d.art3d as art3d 
import numpy as np 
from scipy.interpolate import interp1d 

md,wellz,wellx,welly=np.genfromtxt("./well.csv",delimiter=",",unpack=True) 

# Building interpolation function that map a measured depth to its repsective x,y,z coordinates 
fz = interp1d(md,wellz) 
fx = interp1d(md,wellx) 
fy = interp1d(md,welly) 

pointDepth = np.array([15790,15554,15215,14911,14274,13927,13625,13284,12983,12640,12345,12004,11704,11361,11061,10717,10418,10080,9771]) 

# Some random radii 
dist = [random.random()*100 for x in pointDepth] 

def rotation_matrix(d): 
    """ 
Calculates a rotation matrix given a vector d. The direction of d 
corresponds to the rotation axis. The length of d corresponds to 
the sin of the angle of rotation. 

Variant of: http://mail.scipy.org/pipermail/numpy-discussion/2009-March/040806.html 
    """ 
    sin_angle = np.linalg.norm(d) 

    if sin_angle == 0: 
    return np.identity(3) 

    d = d/sin_angle 

    eye = np.eye(3) 
    ddt = np.outer(d, d) 
    skew = np.array([[ 0, d[2], -d[1]], 
        [-d[2],  0, d[0]], 
        [d[1], -d[0], 0]], dtype=np.float64) 

    M = ddt + np.sqrt(1 - sin_angle**2) * (eye - ddt) + sin_angle * skew 
    return M 

def pathpatch_2d_to_3d(pathpatch, z = 0, normal = 'z'): 
    """ 
    Transforms a 2D Patch to a 3D patch using the given normal vector. 

    The patch is projected into they XY plane, rotated about the origin 
    and finally translated by z. 
    """ 
    if type(normal) is str: #Translate strings to normal vectors 
     index = "xyz".index(normal) 
     normal = np.roll((1,0,0), index) 

    normal = normal/np.linalg.norm(normal) #Make sure the vector is normalised 

    path = pathpatch.get_path() #Get the path and the associated transform 
    trans = pathpatch.get_patch_transform() 

    path = trans.transform_path(path) #Apply the transform 
    pathpatch.__class__ = art3d.PathPatch3D #Change the class 
    pathpatch._code3d = path.codes #Copy the codes 
    pathpatch._facecolor3d = pathpatch.get_facecolor #Get the face color  

    verts = path.vertices #Get the vertices in 2D 

    d = np.cross(normal, (0, 0, 1)) #Obtain the rotation vector 
    M = rotation_matrix(d) #Get the rotation matrix 
    pathpatch._segment3d = np.array([np.dot(M, (x, y, 0)) + (0, 0, z) for x, y in verts]) 

def pathpatch_translate(pathpatch, delta): 
    """ 
    Translates the 3D pathpatch by the amount delta. 
    """ 
    pathpatch._segment3d += delta 


fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d') 
ax.plot(wellx,welly,wellz,c='k') 

for n,pd in enumerate(pointDepth): 
    x,y,z = fx(pd),fy(pd),fz(pd) 

    r = dist[n] 

    # figure out a vector from the point 
    vx,vy,vz = (fx(pd-10)-x),(fy(pd-10)-y),(fz(pd-10)-z) 

    #Draw Circle 
    p = Circle((x,y), r) 
    ax.add_patch(p) 
    pathpatch_2d_to_3d(p, z=z,normal=(vx,vy,vz)) 
    difs = (x,y,z)-p._segment3d[0] 
    pathpatch_translate(p,(difs[0]-r/2,difs[1]-r/2,difs[2]-r/2)) 


ax.set_xlim3d(np.min(wellx),np.max(wellx)) 
ax.set_ylim3d(np.min(welly), np.max(welly)) 
ax.set_zlim3d(np.min(wellz), np.max(wellz)) 
plt.show() 

enter image description here

dienen