Finding quaternion representing the rotation from one vector to another

Question

I have two vectors u and v. Is there a way of finding a quaternion representing the rotation from u to v?

Uniparous

See <laustep.github.io/stlahblog/posts/ReorientTransformation2.html>. – Karlow

Answer 1

Quaternion q;
vector a = crossproduct(v1, v2);
q.xyz = a;
q.w = sqrt((v1.Length ^ 2) * (v2.Length ^ 2)) + dotproduct(v1, v2);

Don't forget to normalize q.

Richard is right about there not being a unique rotation, but the above should give the "shortest arc," which is probably what you need.

Arraign

Be aware that this does not handle the case of parallel vectors (both in the same direction or pointing in opposite directions). crossproduct will not be valid in these cases, so you first need to check dot(v1, v2) > 0.999999 and dot(v1, v2) < -0.999999, respectively, and either return an identity quat for parallel vectors, or return a 180 degree rotation (about any axis) for opposite vectors. – Lydia

@Lydia Actually, if v1 = v2, crossproduct would be (0,0,0) and w would be positive, which normalizes to identity. According to gamedev.net/topic/… it should work just fine also for v1 = -v2 and in their close vicinity. – Contraceptive

How has anyone got this technique to work? For one, sqrt((v1.Length ^ 2) * (v2.Length ^ 2)) simplifies to v1.Length * v2.Length. I couldn't get any variation of this to produce sensible results. – Vaish

Yes, this works. See source code. L61 handles if the vectors face opposite directions (return PI, otherwise it'd return identity per @jpa's remark). L67 handles parallel vectors: mathematically unnecessary, but faster. L72 is Polaris878's answer, assuming both vectors are unit length (avoids a sqrt). See also unit tests. – Lydia

@Lydia You're right, it does work, although the power of two and the sqrt are superfluous. And I've managed to figure out what this equation is doing - it's taking the average of a zero-rotation quaternion and a quaternion for rotating twice the required angle, which, when normalized, produces the desired rotation. In essence, it finds the half-way quaternion, instead of the half-way vector like my method. This method may actually be marginally more efficient, as it seems to require fewer floating point operations. – Vaish

Here's some coffeescript via lolengine.net/blog/2013/09/18/…: THREE.Quaternion.prototype.setFromVectors = (a, b)-> axis = (new THREE.Vector3).crossVectors(a, b) @set(axis.x, axis.y, axis.z, 1 + a.dot(b)) @normalize() @ – Candlestand

In trigonometric terms, the vector part is vec=|v1|*|v2|*sin(arc)*axis, where axis is the rotation axis obtained as normalized cross product. The scalar part is w=|v1|*|v2|*(1+cos(arc)). Now since sin(arc)=2*sin(arc/2)*cos(arc/2) and 1+cos(arc)=2*cos(arc/2)^2, the computed vector is the required unit quaternion with half the angle, scaled by the factor 2*|v1|*|v2|*cos(arc/2). Which means that the procedure of Peter Ehrlich only works if a and b are already normalized. Which is nicely explained in his blog post. – Fullbodied

Can this be simplified if v2 is 0,0,0 and/or v1 is a normalized vector? – Miscegenation

@Miscegenation The question makes not sense if v2 is 0,0,0. If v1 and v2 are normalized, then the last line reduces to q.w = 1 + dotproduct(v1, v2); – Unveiling

@JosephThomson It is usually faster to get v.Length^2 than v.Length itself, that is why I believe it was worded like this. – Torruella

Approach in the answer DOES NOT work! – Clippers

What is v1.Length? – Corked

I believe v1.Length is supposed to be the magnitude of v1 @JuneWang. It took me several minutes of looking at it before the word choice made sense. – Disenable

Answer 2

Half-Way Vector Solution

I came up with the solution that I believe Imbrondir was trying to present (albeit with a minor mistake, which was probably why sinisterchipmunk had trouble verifying it).

Given that we can construct a quaternion representing a rotation around an axis like so:

q.w == cos(angle / 2)
q.x == sin(angle / 2) * axis.x
q.y == sin(angle / 2) * axis.y
q.z == sin(angle / 2) * axis.z

And that the dot and cross product of two normalized vectors are:

dot     == cos(theta)
cross.x == sin(theta) * perpendicular.x
cross.y == sin(theta) * perpendicular.y
cross.z == sin(theta) * perpendicular.z

Seeing as a rotation from u to v can be achieved by rotating by theta (the angle between the vectors) around the perpendicular vector, it looks as though we can directly construct a quaternion representing such a rotation from the results of the dot and cross products; however, as it stands, theta = angle / 2, which means that doing so would result in twice the desired rotation.

One solution is to compute a vector half-way between u and v, and use the dot and cross product of u and the half-way vector to construct a quaternion representing a rotation of twice the angle between u and the half-way vector, which takes us all the way to v!

There is a special case, where u == -v and a unique half-way vector becomes impossible to calculate. This is expected, given the infinitely many "shortest arc" rotations which can take us from u to v, and we must simply rotate by 180 degrees around any vector orthogonal to u (or v) as our special-case solution. This is done by taking the normalized cross product of u with any other vector not parallel to u.

Pseudo code follows (obviously, in reality the special case would have to account for floating point inaccuracies -- probably by checking the dot products against some threshold rather than an absolute value).

Also note that there is no special case when u == v (the identity quaternion is produced -- check and see for yourself).

// N.B. the arguments are _not_ axis and angle, but rather the
// raw scalar-vector components.
Quaternion(float w, Vector3 xyz);

Quaternion get_rotation_between(Vector3 u, Vector3 v)
{
  // It is important that the inputs are of equal length when
  // calculating the half-way vector.
  u = normalized(u);
  v = normalized(v);

  // Unfortunately, we have to check for when u == -v, as u + v
  // in this case will be (0, 0, 0), which cannot be normalized.
  if (u == -v)
  {
    // 180 degree rotation around any orthogonal vector
    return Quaternion(0, normalized(orthogonal(u)));
  }

  Vector3 half = normalized(u + v);
  return Quaternion(dot(u, half), cross(u, half));
}

The orthogonal function returns any vector orthogonal to the given vector. This implementation uses the cross product with the most orthogonal basis vector.

Vector3 orthogonal(Vector3 v)
{
    float x = abs(v.x);
    float y = abs(v.y);
    float z = abs(v.z);

    Vector3 other = x < y ? (x < z ? X_AXIS : Z_AXIS) : (y < z ? Y_AXIS : Z_AXIS);
    return cross(v, other);
}

Half-Way Quaternion Solution

This is actually the solution presented in the accepted answer, and it seems to be marginally faster than the half-way vector solution (~20% faster by my measurements, though don't take my word for it). I'm adding it here in case others like myself are interested in an explanation.

Essentially, instead of calculating a quaternion using a half-way vector, you can calculate the quaternion which results in twice the required rotation (as detailed in the other solution), and find the quaternion half-way between that and zero degrees.

As I explained before, the quaternion for double the required rotation is:

q.w   == dot(u, v)
q.xyz == cross(u, v)

And the quaternion for zero rotation is:

q.w   == 1
q.xyz == (0, 0, 0)

Calculating the half-way quaternion is simply a matter of summing the quaternions and normalizing the result, just like with vectors. However, as is also the case with vectors, the quaternions must have the same magnitude, otherwise the result will be skewed towards the quaternion with the larger magnitude.

A quaternion constructed from the dot and cross product of two vectors will have the same magnitude as those products: length(u) * length(v). Rather than dividing all four components by this factor, we can instead scale up the identity quaternion. And if you were wondering why the accepted answer seemingly complicates matters by using sqrt(length(u) ^ 2 * length(v) ^ 2), it's because the squared length of a vector is quicker to calculate than the length, so we can save one sqrt calculation. The result is:

q.w   = dot(u, v) + sqrt(length_2(u) * length_2(v))
q.xyz = cross(u, v)

And then normalize the result. Pseudo code follows:

Quaternion get_rotation_between(Vector3 u, Vector3 v)
{
  float k_cos_theta = dot(u, v);
  float k = sqrt(length_2(u) * length_2(v));

  if (k_cos_theta / k == -1)
  {
    // 180 degree rotation around any orthogonal vector
    return Quaternion(0, normalized(orthogonal(u)));
  }

  return normalized(Quaternion(k_cos_theta + k, cross(u, v)));
}

Vaish

+1: Great! This worked as a charm. Should be the accepted answer. – Secretive

Quaternion syntaxis is switched on some examples (Quaternion(xyz, w) and Quaternion(w, xyz)). Also seems that in last code block radians and degrees are mixed to express angles (180 vs. k_cos_theta + k). – Necessity

Quaternion(float, Vector3) is construction from scalar-vector, whereas Quaternion(Vector3, float) is construction from axis-angle. Perhaps potentially confusing, but I think it is correct. Correct me if you still think it is wrong! – Vaish

It worked! Thanks! However, I found another similar and well explained link to perform above operation. Thought I should share for the record ;) – Fruitarian

@JosephThomson The half-way quaternion solution seems to come from here. – Anthracene

In the above implementations, the calls to axis-angle constructor, which in turn would call sin and cos can be avoided since sin (180/2) = 1 and cos (180/2) = 0. Thus the scalar-vector constructor can be called instead: return Quaternion(0.0f, normalized(cross(u, other))); – Anthracene

Just want to note that the condition abs(dot(u, X_AXIS)) < 1.0 is wrong here: unlike the rest of the code, it doesn't take into account that vector u is not normalized. E.g. for u=(0.5, 0, 0) this condition will be true, although this vector is parallel to the x axis. Also, if it was normalized, we could just take its x component instead of calculating the dot product, which gives the same because X_AXIS=(1,0,0). Same with Y_AXIS. – Oldline

@Anthracene Thanks. I have amended the answer. I hope everything is correct now. – Vaish

@Oldline Thanks. I have amended the answer. I hope everything is correct now. – Vaish

Half-way vector solution doesn't work for vectors [3, 4, 0], [-3 + 1e-15, -4, 0]. The result is [-0.8; 0. 0. -0.6]. As you can see, this quaternion fails to rotate the vector by 180 degrees. Instead It rotates the vector by about 286.26 degrees. I have used double precision, but the similar vectors of floats can be used. For example, you can try vectors [3, 4, 0] and [-3 + 1e-7, -4, 0]. – Unveiling

@user502144 Did you check u == -v using a threshold? You are going to get floating point error when this condition is nearly but not quite true. Use the second solution if it handles these situations better (though you should still check using a threshold). – Vaish

Answer 3

The problem as stated is not well-defined: there is not a unique rotation for a given pair of vectors. Consider the case, for example, where u = <1, 0, 0> and v = <0, 1, 0>. One rotation from u to v would be a pi / 2 rotation around the z-axis. Another rotation from u to v would be a pi rotation around the vector <1, 1, 0>.

Reiter

In fact isn't there an infinite number of possible answers? Because after you align the "from" vector with the "to" vector you can still freely spin the result around it's axis? Do you know what extra information can typically be used to constrain this choice and make the problem well defined? – Flotsam

Answer 4

I'm not much good on Quaternion. However I struggled for hours on this, and could not make Polaris878 solution work. I've tried pre-normalizing v1 and v2. Normalizing q. Normalizing q.xyz. Yet still I don't get it. The result still didn't give me the right result.

In the end though I found a solution that did. If it helps anyone else, here's my working (python) code:

def diffVectors(v1, v2):
    """ Get rotation Quaternion between 2 vectors """
    v1.normalize(), v2.normalize()
    v = v1+v2
    v.normalize()
    angle = v.dot(v2)
    axis = v.cross(v2)
    return Quaternion( angle, *axis )

A special case must be made if v1 and v2 are paralell like v1 == v2 or v1 == -v2 (with some tolerance), where I believe the solutions should be Quaternion(1, 0,0,0) (no rotation) or Quaternion(0, *v1) (180 degree rotation)

Pollack

I have a working implementation, but this yours is prettier, so I really wanted it to work. Unfortunately it failed all of my test cases. My tests all look something like quat = diffVectors(v1, v2); assert quat * v1 == v2. – Lydia

It’s unlikely that this will work at all since angle gets its value from a dot product. – Hinman

Where is the Quaternion() function? – Corked

I haven't tried this, but, looking at it, I think maybe you just need to remove the v.normalize(). So the scalar part of the answer will be v.dot(v2) = (v1+v2).dot(v2) = 1 + v1.dot(v2), and the vector part will be v.cross(v2) = (v1+v2).cross(v2) = v1.cross(v2). – Unpaged

Answer 5

To find the quaternion of smallest rotation which rotates u to v, use

align(quat(1, 0, 0, 0), u, v)

The Generalized Solution

function align(Q, u, v)
    U = quat(0, ux, uy, uz)
    V = quat(0, vx, vy, vz)
    return normalize(length(U*V)*Q - V*Q*U)

Why This Generalization?

R = align(Q, u, v)

Most importantly, R is the quaternion closest to Q whose local u direction points globally in the v direction. Consequently, R is the quaternion closest to Q which will rotate u to v.

This can be used to give you all possible rotations which rotate from u to v, depending on the choice of Q. If you want the minimal rotation from u to v, as the other solutions give, use Q = quat(1, 0, 0, 0).

Most commonly, I find that the real operation you want to do is a general alignment of one axis with another.

// If you find yourself often doing something like
quatFromTo(toWorldSpace(Q, localFrom), worldTo)*Q
// you should instead consider doing 
align(Q, localFrom, worldTo)

Example

Say you want the quaternion Y which only represents Q's yaw, the pure rotation about the y axis. We can compute Y with the following.

Y = align(quat(Qw, Qx, Qy, Qz), vec(0, 1, 0), vec(0, 1, 0))

// simplifies to
Y = normalize(quat(Qw, 0, Qy, 0))

Alignment as a 4x4 Projection Matrix

If you want to perform the same alignment operation repeatedly, because this operation is the same as the projection of a quaternion onto a 2D plane embedded in 4D space, we can represent this operation as the multiplication with 4x4 projection matrix, A*Q.

A = I - leftQ(V)*rightQ(U)/length(U*V)

// which expands to
A = mat4(
    1 + ux*vx + uy*vy + uz*vz, uy*vz - uz*vy, uz*vx - ux*vz, ux*vy - uy*vx,
    uy*vz - uz*vy, 1 + ux*vx - uy*vy - uz*vz, uy*vx + ux*vy, uz*vx + ux*vz,
    uz*vx - ux*vz, uy*vx + ux*vy, 1 - ux*vx + uy*vy - uz*vz, uz*vy + uy*vz,
    ux*vy - uy*vx, uz*vx + ux*vz, uz*vy + uy*vz, 1 - ux*vx - uy*vy + uz*vz)

// A can be applied to Q with the usual matrix-vector multiplication
R = normalize(A*Q)

//LeftQ is a 4x4 matrix which represents the multiplication on the left
//RightQ is a 4x4 matrix which represents the multiplication on the Right
LeftQ(w, x, y, z) = mat4(
    w, -x, -y, -z,
    x,  w, -z,  y,
    y,  z,  w, -x,
    z, -y,  x,  w)

RightQ(w, x, y, z) = mat4(
    w, -x, -y, -z,
    x,  w,  z, -y,
    y, -z,  w,  x,
    z,  y, -x,  w)

I = mat4(
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 1)

Nide

This is actually amazing! Just sad I didn't learn about this a decade or so earlier, would have helped alot with imu sensor data processing. its funny how Madgwick filter uses gradient descent optimization steps for the same problem, but now I learn that you can just align the quaternion in one shot. Really helpful! – Surplusage

Answer 6

From algorithm point of view , the fastest solution looks in pseudocode

 Quaternion shortest_arc(const vector3& v1, const vector3& v2 ) 
 {
     // input vectors NOT unit
     Quaternion q( cross(v1, v2), dot(v1, v2) );
     // reducing to half angle
     q.w += q.magnitude(); // 4 multiplication instead of 6 and more numerical stable

     // handling close to 180 degree case
     //... code skipped 

        return q.normalized(); // normalize if you need UNIT quaternion
 }

Be sure that you need unit quaternions (usualy, it is required for interpolation).

NOTE: Nonunit quaternions can be used with some operations faster than unit.

Sava

Answer 7

Why not represent the vector using pure quaternions? It's better if you normalize them first perhaps.
q₁ = (0 u_x u_y u_z)'
q₂ = (0 v_x v_y v_z)'
q₁ q_rot = q₂
Pre-multiply with q₁^-1
q_rot = q₁^-1 q₂
where q₁^-1 = q₁^conj / q_norm
This is can be thought of as "left division". Right division, which is not what you want is:
q_rot,right = q₂^-1 q₁

Goiter

I'm lost, isn't the rotation from q1 to q2 calculated as q_2 = q_rot q_1 q_rot^-1 ? – Clarethaclaretta

You are right. I've tried this, and it is not working – Sliver

Answer 8

Some of the answers don't seem to consider possibility that cross product could be 0. Below snippet uses angle-axis representation:

//v1, v2 are assumed to be normalized
Vector3 axis = v1.cross(v2);
if (axis == Vector3::Zero())
    axis = up();
else
    axis = axis.normalized();

return toQuaternion(axis, ang);

The toQuaternion can be implemented as follows:

static Quaternion toQuaternion(const Vector3& axis, float angle)
{
    auto s = std::sin(angle / 2);
    auto u = axis.normalized();
    return Quaternion(std::cos(angle / 2), u.x() * s, u.y() * s, u.z() * s);
}

If you are using Eigen library, you can also just do:

Quaternion::FromTwoVectors(from, to)

Benedictine

toQuaternion(axis, ang) -> you forgot to specify what is ang – Schubert

2nd parameter is angle which is part of axis-angle representation of the quaternion, measured in radians. – Benedictine

You were asked to get quaternion to rotate from one vector to another. You don't have angle, you have to calculate it first. Your answer should contain the calculation of angle. Cheers! – Schubert

This is c++? what is u.x()? – Corked

Yes, this is C++. u is vector type from Eigen library (if you are using one). – Benedictine

Answer 9

One can rotate a vector u to vector v with

function fromVectors(u, v) {

  d = dot(u, v)
  w = cross(u, v)

  return Quaternion(d + sqrt(d * d + dot(w, w)), w).normalize()
}

If it is known that the vectors u to vector v are unit vectors, the function reduces to

function fromUnitVectors(u, v) {
  return Quaternion(1 + dot(u, v), cross(u, v)).normalize()
}

Depending on your use-case, handling the cases when the dot product is 1 (parallel vectors) and -1 (vectors pointing in opposite directions) may be needed.

Halyard

Honestly the best way to check for this case is just to check the magnitude before normalizing. As long as it's not 0, should be fine. The way magnitude is computed requires squaring values. So anything that requires closer than half of machine precision has a magnitude of 0. with floats, 10^-20 * 10^-20 -> 0 – Nide

Answer 10

Working just with normalized quaternions, we can express Joseph Thompson's answer in the follwing terms.

Let q_v = (0, u_x, v_y, v_z) and q_w = (0, v_x, v_y, v_z) and consider

q = q_v * q_w = (-u dot v, u x v).

So representing q as q(q_0, q_1, q_2, q_3) we have

q_r = (1 - q_0, q_1, q_2, q_3).normalize()

Indoor

Answer 11

@RicharDunlap is right: there's an infinity of solutions. And saying the "smallest rotation" is not clear.

But there is a "canonical" quaternion which maps a normalized vector u to a normalized vector v. Precisely, it maps the plane passing through the origin with normal vector u to the plane passing through the origin with normal vector v. With such a description, it is unique.

This quaternion is:

re = sqrt((1 + sum(u*v))/2)
w = crossProduct(u, v) / 2 / re
q = (re, w[1], w[2], w[3])

Afro