Matplotlib Rotating 3d Disk
<p>I have a line and some points that are on that are on that line in 3D space. I know there is a certain amount of error in the the point, but the error only extends perpendicular to the line.    To view this I'd like to have disks with a radius of the error that are centered on the line and orthogonal to the direction of the line.  I found this <a href="http://stackoverflow.com/questions/18228966/how-can-matplotlib-2d-patches-be-transformed-to-3d-with-arbitrary-normals">solution</a> but I can't get it to quite work. </p>

<p>If I run the code and state I want the output normal to the 'z' axis I get waht I woudl expect. Disks with a given radius on the line and oriented on z axis. </p>

<pre><code>pathpatch_2d_to_3d(p, z=z,normal='z')
</code></pre>

<p><a href="http://i.stack.imgur.com/9E4OI.png" rel="nofollow"><img src="http://i.stack.imgur.com/9E4OI.png" alt="Image with normal=&#39;z&#39;"></a></p>

<p>I would like the disks rotated. In order to find the well vector at that point I'm using a point close by using that vector.  This is the vector that I am putting as <code>normal=(vx,vy,vz)</code> but when I do that, the disks are not even on the chart.  I'm not certain where I'm going wrong.  Does anyone have any advice? </p>

<p>Here is my code.</p>

<pre><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()
</code></pre>

<p>Here is the solution I've come up with.  I decided to take the difference between where the point falls on the line and where the first point in <code>p.._segment3d</code>.  This gives me how far away my circle is from where I want it to be, then I simply translated the patch that distance minus the radius of the circle so it would be centered.</p>

<p>I've added in some random numbers to serve as some "error" and here is final code and resulting image</p>

<pre><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])

# 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()
</code></pre>

<p><a href="http://i.stack.imgur.com/FxJoU.png" rel="nofollow"><img src="http://i.stack.imgur.com/FxJoU.png" alt="enter image description here"></a></p>

python matplotlib 3d