Rotations in 3D

Question

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?

There is a question about a "Good geometry library in python?" that is more general in scope:
http://stackoverflow.com/questions/1076778/good-geometry-library-in-python — Gabriel Schreiber, Oct 30 '12 at 10:50
Instead of writing another SO-question as an answer it would be better to mark this question as duplicate of the other question. ;) — erikbstack, Oct 30 '12 at 11:03
@erikb85 The other question is much more general in scope; my question may have an answer that doesn't apply to Stefano Borinis question. — Gabriel Schreiber, Oct 30 '12 at 11:15
Then it should be comment and you could edit your question to embrace the difference more clearly. Yet, I think a scientific/engineering approach to solving problems is to transform a problem to a more general one that already has a solution, then apply the general solution. — erikbstack, Nov 01 '12 at 12:00

score 1 · Answer 1 · answered Jan 26 '13 at 01:11

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

score 0 · Answer 2 · answered Nov 06 '12 at 00:14

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.

Rotations in 3D

2 Answers2