Matplotlib 3Dディスクを回転させる

Question

私には線があり、その上にある点が3D空間上のその線上にあります。私はその点にある程度の誤差があることを知っていますが、誤差は線に対して垂直にしか伸びません。これを表示するには、エラーの半径が線の中央にあり、線の方向と直交するディスクが必要です。私はこれを見つけましたsolutionしかし、私はそれをかなり仕事に得ることができません。

私がコードを実行し、出力を 'z'軸に垂直にしたい場合、私はwoudlを期待しています。指定された半径の円上にあり、z軸に向けられたディスク。

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

私は、回転ディスクをしたいと思います。その点で井戸ベクトルを見つけるために、私はそのベクトルを使って点を近似しています。これは私がnormal=(vx,vy,vz)としているベクトルですが、私がそれを行うと、ディスクはグラフ上にさえありません。私はどこが間違っているのか分からない。誰にもアドバイスはありますか？

ここに私のコードです。

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

Answer 1

これは私が思いついた解決策です。私は、ポイントがライン上にある場所と最初のポイントがどこにあるかの違いを取ることにしました。p.._segment3d。これは私の円が私が望むところからどのくらい離れているかを教えてくれます。そして、その距離からその円の半径を差し引いたものを中心になるように翻訳しました。

私はここにいくつかの「エラー」として機能するためにいくつかのランダムな数字に追加されましたが、最終的なコード、そして得られた画像である

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