У меня есть строка и некоторые точки, которые находятся на этой линии в 3D-пространстве. Я знаю, что в этой точке есть определенная ошибка, но ошибка распространяется только перпендикулярно линии. Для этого я хотел бы иметь диски с радиусом ошибки, которые центрированы по прямой и ортогональны направлению линии. Я нашел это solution, но я не могу заставить его работать.Matplotlib Вращающийся 3d-диск
Если я запустил код и состояние, я хочу, чтобы выход нормалился к оси «z», я получаю waht, который я ожидаю. Диски с заданным радиусом на линии и ориентированы по оси 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()
