【发布时间】：2016-09-16 05:06:50
【问题描述】：
我有一条线，在 3D 空间中，这条线上的一些点位于该线上。我知道该点存在一定的误差，但误差仅垂直于线延伸。要查看这一点，我希望磁盘的误差半径以直线为中心并与直线的方向正交。我找到了这个solution，但我无法让它完全正常工作。
如果我运行代码并声明我希望输出垂直于“z”轴，我会得到我期望的结果。在线上具有给定半径并沿 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()
【问题讨论】：
标签： python matplotlib 3d