Search code examples
python3dcoordinate-transformationrotational-matrices

Coordinate transformations from a randomly generated normal vector


I'm trying to randomly generate coordinate transformations for a fitting routine I'm writing in python. I want to rotate my data (a bunch of [x,y,z] coordinates) about the origin, ideally using a bunch of randomly generated normal vectors I've already created to define planes -- I just want to shift each plane I've defined so that it lies in the z=0 plane.

Here's a snippet of my code that should take care of things once I have my transformation matrix. I'm just not sure how to get my transformation matrix from my normal vector and if I need something more complicated than numpy for this.

import matplotlib as plt
import numpy as np
import math

origin = np.array([35,35,35])
normal = np.array([np.random.uniform(-1,1),np.random.uniform(-1,1),np.random.uniform(0,1)])
mag = np.sum(np.multiply(normal,normal))
normal = normal/mag

a = normal[0]
b = normal[1]
c = normal[2]

#I know this is not the right transformation matrix but I'm not sure what is...
#Looking for the steps that will take me from the normal vector to this transformation matrix
rotation = np.array([[a, 0, 0], [0, b, 0], [0, 0, c]])

#Here v would be a datapoint I'm trying to shift?
v=(test_x,test_y,test_z)
s = np.subtract(v,origin) #shift points in the plane so that the center of rotation is at the origin
so = np.multiply(rotation,s) #apply the rotation about the origin
vo = np.add(so,origin) #shift again so the origin goes back to the desired center of rotation

x_new = vo[0]
y_new = vo[1]
z_new = vo[2]

fig = plt.figure(figsize=(9,9))
plt3d = fig.gca(projection='3d')
plt3d.scatter(x_new, y_new, z_new, s=50, c='g', edgecolor='none')

Solution

  • Thanks to the people over in the math stack exchange, I have an answer that works. But note that it would not work if you also needed to perform a translation, which I didn't because I'm defining my planes by a normal vector and a point, and the normal vector changes but the point does not. Here's what worked for me.

    import matplotlib as plt
    import numpy as np
    import math
    
    def unit_vector(vector):
        """ Returns the unit vector of the vector.  """
        return vector / np.linalg.norm(vector)
    
    cen_x, cen_y, cen_z = 35.112, 35.112, 35.112
    origin = np.array([[cen_x,cen_y,cen_z]])
    
    z_plane_norm = np.array([1,1,0])
    z_plane_norm = unit_vector(z_plane_norm)
    
    normal = np.array([np.random.uniform(-1,1),np.random.uniform(-1,1),np.random.uniform(0,1)])
    normal = unit_vector(normal)
    
    a1 = normal[0]
    b1 = normal[1]
    c1 = normal[2]
    
    rot = np.matrix([[b1/math.sqrt(a1**2+b1**2), -1*a1/math.sqrt(a1**2+b1**2), 0], [a1*c1/math.sqrt(a1**2+b1**2), b1*c1/math.sqrt(a1**2+b1**2), -1*math.sqrt(a1**2+b1**2)], [a1, b1, c1]])
    
    init = np.matrix(normal)
    
    fin = rot*init.T
    fin = np.array(fin)
    
    # equation for a plane is a*x+b*y+c*z+d=0 where [a,b,c] is the normal
    # so calculate d from the normal
    d1 = -origin.dot(normal)
    
    # create x,y
    xx, yy = np.meshgrid(np.arange(cen_x-0.5,cen_x+0.5,0.05),np.arange(cen_y-0.5,cen_y+0.5,0.05))
    
    # calculate corresponding z
    z1 = (-a1 * xx - b1 * yy - d1) * 1./c1
    
    #-------------
    
    a2 = fin[0][0]
    b2 = fin[1][0]
    c2 = fin[2][0]
    
    d2 = -origin.dot(fin)
    d2 = np.array(d2)
    d2 = d2[0][0]
    
    z2 = (-a2 * xx - b2 * yy - d2) * 1./c2
    
    #-------------
    
    # plot the surface
    fig = plt.figure(figsize=(9,9))
    plt3d = fig.gca(projection='3d')
    
    plt3d.plot_surface(xx, yy, z1, color='r', alpha=0.5, label = "original")   
    plt3d.plot_surface(xx, yy, z2, color='b', alpha=0.5, label = "rotated")  
    
    
    plt3d.set_xlabel('X (Mpc)')
    plt3d.set_ylabel('Y (Mpc)')
    plt3d.set_zlabel('Z (Mpc)')
    
    plt.show()
    

    If you do need to perform a translation as well, see the full answer I worked off of here.