Gabriel Schreiber
Gabriel Schreiber

Reputation: 2226

Rotations in 3D

I have three vectors in 3D a,b,c. Now I want to calculate a rotation r that when applied to a yields a result parallel to b. Then the rotation r needs to be applied to c.

How do I do this in python? Is it possible to do this with numpy/scipy?

Upvotes: 1

Views: 1820

Answers (2)

ionelberdin
ionelberdin

Reputation: 55

Using numpy:

import numpy as np
from numpy import linalg as LA
from math import pi,cos,sin,acos 

def rotate(v,angle=0,ref=np.array([0,0,1]),deg=False):
    '''Rotates a vector a given angle respect the 
        given reference. Option: deg=False (default)'''
    if(abs(angle) < 1e-5):
        return v
    if(deg):
        angle = angle*pi/180
    # Define rotation reference system
    ref = versor(ref) # rotation axis
    # n1 & n2 are perpendicular to ref, and to themselves 
    n1 = versor(np.cross(ref,np.array([-ref[1],ref[2],ref[0]])))
    n2 = np.cross(ref,n1)
    vp = np.inner(v,ref)*ref # parallel to ref vector
    vn = v-vp # perpendicular to ref vector
    vn_abs = LA.norm(vn)
    if(vn_abs < 1e-5):
        return v
    alp = acos(np.inner(vn,n1)/vn_abs) # angle between vn & n1
    if(triprod(ref,n1,vn) < 0):
        alp = -alp # correct if necesary
    return vp+vn_abs*(n1*cos(alp+angle)+n2*sin(alp+angle))

def triprod(a,b,c):
    '''Triple product of vectors: a·(b x c)'''
    return np.inner(a,np.cross(b,c))

def versor(v):
    '''Unitary vector in the direction of the one given as input'''
    v = np.array(v)
    return v/LA.norm(v)

###### Test ################################################
a = np.array([3,4,1])
b = np.array([0,-1,2])
c = np.array([1,1,5])
r = acos(np.inner(a,b)/(LA.norm(a)*LA.norm(b)))
ref = versor(np.cross(a,b))
print rotate(c,angle=r,ref=ref)
print r
print ref

Upvotes: 1

eh9
eh9

Reputation: 7430

I'll assume the "geometry library for python" already answered in the comments on the question. So once you have a transformation that takes 'a' parallel to 'b', you'll just apply it to 'c'

The vectors 'a' and 'b' uniquely define a plane. Each vector has a canonical representation as a point difference from the origin, so you have three points: the head of 'a', the head of 'b', and the origin. First compute this plane. It will have an equation in the form Ax + By + Cz = 0.

A normal vector to this plane defines both the axis of rotation and the sign convention for the direction of rotation. All you need is one normal vector to the plane, since they're all collinear. You can solve for such a vector by picking at two non-collinear vectors in the plane and taking the dot product with the normal vector. This gives a pair of linear equations in two variables that you can solve with standard methods such as Cramer's rule. In all of these manipulations, if any of A, B, or C are zero, you have a special case to handle.

The angle of the rotation is given by the cosine relation for the dot product of 'a' and 'b' and their lengths. The sign of the angle is determined by the triple product of 'a', 'b', and the normal vector. Now you've got all the data to construct a rotation matrix in one of the many canonical forms you can look up.

Upvotes: 0

Related Questions