Beautiful maths simplification: quaternion from two vectors

In this article I would like to explain the thought process behind the derivation of a widely used formula: finding a quaternion representing the rotation between two 3D vectors. Nothing really new, but hopefully a few ideas could be reused at other times.

Note: the routine presented here is incomplete on purpose. For a version that can be used in production code, see the next article instead, Quaternion from two vectors: the final version.

Naive method

A rotation is best visualised using a rotation axis and an angle. Except in degenerate cases, the rotation axis can be obtained by computing the cross product of the two original vectors:

Then the angle can be obtained using the properties of the cross product and/or the dot product:

\begin{align}
\boldsymbol{u} . \boldsymbol{v} &= |\boldsymbol{u}| . |\boldsymbol{v}| . \cos\theta \\
||\boldsymbol{u} \times \boldsymbol{v}|| &= |\boldsymbol{u}| . |\boldsymbol{v}| . |\sin\theta|
\end{align}

Since θ is always between 0 and π, we only care about the dot product. This gives us some obvious code to create the quaternion (omitting corner cases such as θ = 0 for clarity):

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float cos_theta = dot(normalize(u), normalize(v));
    float angle = acos(cos_theta);
    vec3 w = normalize(cross(u, v));
    return quat::fromaxisangle(angle, w);
}

This is naive but it works. Googling for “quaternion from two vectors” shows this forum post and this SO question where this method or a variation thereof is suggested.

Looking under the hood

Let’s have a look at what happens when building the quaternion from an axis and an angle:

$$ q = \cos\frac{\theta}{2} + \boldsymbol{i}\sin\frac{\theta}{2}v_x + \boldsymbol{j}\sin\frac{\theta}{2}v_y + \boldsymbol{k}\sin\frac{\theta}{2}v_z $$

This means the code for quat::fromaxisangle would look somewhat like this:

quat quat::fromaxisangle(float angle, vec3 axis)
{
    float half_sin = sin(0.5f * angle);
    float half_cos = cos(0.5f * angle);
    return quat(half_cos,
                half_sin * axis.x,
                half_sin * axis.y,
                half_sin * axis.z);
}

Avoiding trigonometry

If you read Iñigo Quilez’s recent article about avoiding trigonometry you’ll have probably frowned at the fact that we computed θ from cos(θ), then computed sin(θ/2) and cos(θ/2).

Indeed, it happens that there is a much simpler way to do it; the half-angle formulas from precalculus tell us the following:

\begin{align}
\sin\frac{\theta}{2} &= \sqrt{\frac{1 - \cos\theta}{2}} \\
\cos\frac{\theta}{2} &= \sqrt{\frac{1 + \cos\theta}{2}}
\end{align}

This allows us to simplify our quaternion creation code:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float cos_theta = dot(normalize(u), normalize(v));
    float half_cos = sqrt(0.5f * (1.f + cos_theta));
    float half_sin = sqrt(0.5f * (1.f - cos_theta));
    vec3 w = normalize(cross(u, v));
    return quat(half_cos,
                half_sin * w.x,
                half_sin * w.y,
                half_sin * w.z);
}

This is pretty nice. By using well known trigonometry formulas, we got rid of all trigonometry function calls!

Avoiding square roots

It happens that we can do slightly better. Note that we normalize three vectors: u, v and cross(u, v). That’s three square roots. The thing is, we already know the norm of w through this formula:

$$||\boldsymbol{u} \times \boldsymbol{v}|| = |\boldsymbol{u}| . |\boldsymbol{v}| . |\sin\theta|$$

And we know sin(θ) from precalculus again:

$$\sin\theta = 2 \sin\frac{\theta}{2} \cos\frac{\theta}{2}$$

Also, using the fact that sqrt(a)sqrt(b) = sqrt(ab) lets us perform one less square root.

We can therefore come up with the following performance improvement:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float norm_u_norm_v = sqrt(sqlength(u) * sqlength(v));
    float cos_theta = dot(u, v) / norm_u_norm_v;
    float half_cos = sqrt(0.5f * (1.f + cos_theta));
    float half_sin = sqrt(0.5f * (1.f - cos_theta));
    vec3 w = cross(u, v) / (norm_u_norm_v * 2.f * half_sin * half_cos);
    return quat(half_cos,
                half_sin * w.x,
                half_sin * w.y,
                half_sin * w.z);
}

Oh wait! We divide by sin(θ/2) to compute w, then we multiply by sin(θ/2) again. This means we don’t even need that variable, and we can simplify even further:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float norm_u_norm_v = sqrt(sqlength(u) * sqlength(v));
    float cos_theta = dot(u, v) / norm_u_norm_v;
    float half_cos = sqrt(0.5f * (1.f + cos_theta));
    vec3 w = cross(u, v) / (norm_u_norm_v * 2.f * half_cos);
    return quat(half_cos, w.x, w.y, w.z);
}

This is more or less the code used by the Ogre3D engine in OgreVector3.h, except they perform an additional normalisation step on the final result. This is mathematically useless, but due to numerical stability issues, it is probably safe to do so nonetheless.

This final normalisation step is actually an opportunity to simplify the code even further.

Improving on Ogre3D

We are down to two square roots and four divisions, plus quite a few mul/adds. Depending on the platform that we are running on, it is possible to simplify even further and improve performance. For instance, on many SIMD architectures, normalising a quaternion can be very fast.

This is the code we get if we multiply every component of the quaternion by 2.f * half_cos and let normalize() do the rest of the job:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float norm_u_norm_v = sqrt(sqlength(u) * sqlength(v));
    float cos_theta = dot(u, v) / norm_u_norm_v;
    float half_cos = sqrt(0.5f * (1.f + cos_theta));
    vec3 w = cross(u, v) / norm_u_norm_v;
    return normalize(quat(2.f * half_cos * half_cos, w.x, w.y, w.z));
}

Now half_cos only appears in its squared form, and since it comes from a square root, we can simply omit that square root:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float norm_u_norm_v = sqrt(sqlength(u) * sqlength(v));
    float cos_theta = dot(u, v) / norm_u_norm_v;
    vec3 w = cross(u, v) / norm_u_norm_v;
    return normalize(quat(1.f + cos_theta, w.x, w.y, w.z));
}

And using the same reasoning we can multiply every quaternion component by norm_u_norm_v:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float norm_u_norm_v = sqrt(sqlength(u) * sqlength(v));
    vec3 w = cross(u, v);
    quat q = quat(norm_u_norm_v + dot(u, v), w.x, w.y, w.z);
    return normalize(q);
}

We are still doing two square roots and four divisions, some of which are hidden in normalize(), but the code is considerably shorter now.

Final form

If u and v can be enforced to be unit vectors, norm_u_norm_v can be omitted and simply replaced with 1.0f:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    vec3 w = cross(u, v);
    quat q = quat(1.f + dot(u, v), w.x, w.y, w.z);
    return normalize(q);
}

Isn’t it beautiful, considering the sin(), cos() and acos() ridden mess we started with?

This algorithm can be found all over the Internet, but I do not know who first came up with it. Also, a lot of 3D engines (both publicly available and slightly more private) could benefit from it.

Update (06/01/2014)

In the comments below, Michael Norel provides the following improvement to the non-unit version. Since the values d = dot(u, v) and w = cross(u, v) are computed no matter what, the value sqlength(u) * sqlength(v) could be computed in a different way, *i.e.* d * d + sqlength(w). The following code does at least three multiplications less:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    vec3 w = cross(u, v);
    quat q = quat(dot(u, v), w.x, w.y, w.z);
    q.w += length(q);
    return normalize(q);
}

Also, Marc B. Reynolds notes that, in the unit version, the final normalisation factor is sqrt((1 + dot(u, v))² + sqlength(cross(u, v))) which reduces to sqrt(2 + 2 dot(u, v)) thanks to the sin² + cos² identity. It leads to the following possibly improved version:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    float m = sqrt(2.f + 2.f * dot(u, v));
    vec3 w = (1.f / m) * cross(u, v);
    return quat(0.5f * m, w.x, w.y, w.z);
}
  • Posted: 2013-09-18 18:13 (Updated: 2014-06-01 11:32)
  • Author: sam
  • Categories: maths optim

Attachments (2)

Download all attachments as: .zip

Comments

1. anonymous -- 2013-11-14 10:12

Thank you for sharing this!! It is both beautiful and amazing. Well done.

2. ponce -- 2014-01-23 11:18

Toujours intéressant ces articles :)

3. minorlogic@yahoo.com -- 2014-05-09 13:02

The complete version from nonunit vectors looks like

quat quat::fromtwovectors(vec3 u, vec3 v) {

vec3 w = cross(u, v); quat q = quat(dot(u, v), w.x, w.y, w.z); q.w += q.magnitude(); return normalize(q);

}

  1. Faster. 4 muls for quaternion magnitude instead of 6 muls for vectors magnitude.
  2. Numerical more stable (produce less roundoff error) than use of vectors length.

Derivation of original method can be found http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/minorlogic.htm

4. marc.b.reynolds@gmail.com -- 2014-05-28 11:25

This reduces even more. Assuming unit vector inputs and ignoring the degenerate case:

vec4 getRot(vec3 a, vec3 b)
{
  vec4 r
  float m = sqrt(1+dot(a,b));
  r.xyz = (1/m)*cross(a,b);
  r.w   = 0.5*m;
}
5. anonymous -- 2014-05-28 11:34

Opps: don't see how to edit. should be m = sqrt(2*(1+dot(a,b)))

7. sam -- 2014-05-28 17:16

@marc.b.reynolds — sorry, this blog comment system is quite mediocre. But I think you need m = sqrt(2 + dot(a,b)); because that’s the final normalisation value, and then r.w = sqrt(1 + dot(a,b)) / m; will mean one square root too much.

8. marc.b.reynolds@gmail.com -- 2014-05-29 11:54

Sorry for the formatting. Where am I wrong?

r = ab*

r = cross(a,b) + dot(a,b)

r = sin(t)U + cos(t)

q = sqrt(r)

q = sin(t/2)U + cos(t/2)

w = cos(t) = dot(a,b)

s = sin(t/2)/sin(t)

s = sqrt[(1-w))/2] / sqrt(1-ww)

s = 1/sqrt(2+2w)

w' = cos(t/2)

w' = sqrt[(1+cos(t))/2]

w' = sqrt[(1+w)/2]

w' = sqrt[2+2w]/2

m = sqrt(2+2w) = sqrt(2+2dot(a,b))

q = cross(a,b)/m + m/2

9. sam -- 2014-05-29 12:23

@marc.b.reynolds The problem is that this final quaternion is not a unit quaternion. Its squared length is 1/m² + m²/4 and that’s 1 only when m = sqrt(2)

10. marc.b.reynolds@gmail.com -- 2014-05-30 09:48

As a simple example:

a = (1,0,0)

b = (x,y,0)

r = cross(a,b)+dot(a,b) = (0,0,y) + x

t = sqrt(2+2x) = sqrt(2)sqrt(1+x)

q = (0,0,y/t) + t/2

q = (0,0,y/(sqrt(2)sqrt(1+x)) + sqrt(1+x)/sqrt(2)

Computing 'q' by normalization yields a normalization factor of sqrt[(1+x)(1+x)+yy] = sqrt(1+2x+xx+yy) = sqrt(2+2x)

q = (0,0,y)/sqrt(2+2x) + (1+x)/sqrt(2+2x)

q = (0,0,y)/(sqrt(2)(sqrt(1+x)) + sqrt(1+x)/sqrt(2)

So the two methods are the same. Allowing 'a' & 'b' to be arbitrary unit vectors doesn't change the problem.

11. sam -- 2014-05-30 11:16

Oh, you’re absolutely right. I was making the same mistake from the beginning. Sorry for wasting your time, I’ll add your suggestion to the article!

12. marc.b.reynolds@gmail.com -- 2014-06-04 07:48

It can't be said that I was being clear and concise. The interesting part here is the square root of a unit quaterion which can be useful elsewhere.

The hateful remaining bit of the unit version is the division which can be removed if your language and hardware supports either 1/sqrt(x) or an approximation, then use x/sqrt(x) = sqrt(x) to complete.

quat quat::fromtwovectors(vec3 u, vec3 v)
{
  float d = 1+dot(u,v);
  float m = rsqrt(d+d);
  vec3  w = m * cross(u, v);
  return quat(d*m, w.x, w.y, w.z);
}
13. gfrodo -- 2014-06-04 22:30

in the non-unit version you can even reduce 3 multiplications, because sqlength(w) has to calculated twice: first time in length(q), the second time in normalize(q).

14. gfrodo -- 2014-06-04 22:37

Here is the optimized code:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    vec3 w = cross(u, v);
    quat q = quat(dot(u, v), w.x, w.y, w.z);
    float l=sqlength(w);
    q.w += sqrt(q.w*q.w+l);
    return q*(1/sqrt(q.w*q.w+l));
    //normally one '/' and 4 '*' is faster than 4 '/'
}
15. gfrodo -- 2014-06-04 22:58

additionally, if l is very small and dot(u,v) is negative, you can add a special treatment for opposite vectors (where you can drop one sqrt)

16. gfrodo -- 2014-06-04 23:22

according to the other article here the version checking for 180 degree:

quat quat::fromtwovectors(vec3 u, vec3 v)
{
    vec3 w = cross(u, v);
    float l=sqlength(w);
    float real_part=dot(u, v);
    if(l < 1.6e-12f && real_part < 0 )
    {
        w = abs(u.x) > abs(u.z) ? vec3(-u.y, u.x, 0.f) / sqrt(u.y*u.y + u.x*u.x)
                                : vec3(0.f, -u.z, u.y) / sqrt(u.y*u.y + u.z*u.z);
        return quat(0, w.x, w.y, w.z);
    }
    real_part += sqrt(real_part*real_part + l);
    return quat(real_part, w.x, w.y, w.z) * (1 / sqrt(real_part*real_part + l));
}

Maybe the very special case of opposite vectors can be even more optimized, but possible would be never called.

19. minorlogic@yahoo.com -- 2014-06-17 18:03

@gfrodo your version from "14" is identical to version from "3" , instead of implicit cache of x*x+y*y+z*z.

20. gfrodo -- 2014-06-24 15:00

@minorlogic: yes, you are right. i wanted to be sure, that the compiler doesn't compute x*x+y*y+z*z twice. i think we do this optimizations because we don't trust the compiler to do it.

21. shiqi ai -- 2016-04-17 16:22

its amazing .. thank you for your sharing.

22. Angelica Jayson -- 2017-10-06 08:37

It to a great degree is noted as a revolution. In stargazing, the dissemination of a wonderful physical make-up around its hub. unequivocal the Moon turns. it to some degree is in fundamental terms that its turn and circle around the earth are the equivalents, thus, we in essential terms see the similar edge of the Moon. the correct call is Synodic Month. The time between 2 New Moons. A sidereal month is a period it takes to make a whole 360° hover of the Earth in respect to the celebs. it to some degree is shorter than a synodic month on account of the reality the Earth itself has moved in its turn. http://www.assignmentclock.com/

23. anonymous -- 2017-10-14 08:51

run 3 game http://run23.co has the very best solution of your all problems guys .if you are feeling some boredom then here is a awesome online fun run 3 for you .just click here and enjoy

24. anonymous -- 2017-11-08 10:39

Get free amazon gift card codes from here within minute. http://freeamazongiftcardsnow.com http://onlinefreecodes.com

25. anonymous -- 2017-12-14 16:00
26. anonymous -- 2017-12-26 07:13

Don't worry about your old potty SANDAS help you to solve your former Potty. And if you have to buy a Toilet then HAMRA SANDAS show you the best option.

<a href="https://www.aapkepaikhanekihifazat.com/article">HAMRA SANDAS</a>

27. jonc113@optimum.net -- 2017-12-28 23:11

Apparently, there are TWO DIFFERENT WAYS OF NOTATING QUATERNIONS ???

This forum uses (w, x, y, z), but Unity and many others use (x, y, z, w)

I was scratching my head for over a week until I finally figured this out - and it is mentioned in VERY FEW places on the internet. (Pardon my SHOUTING)

28. Kelvim -- 2018-01-15 08:09

You have made a great tutorial, which is very informative for students. I am impressed to read the content you have posted, i will come back to read more similar articles. Purchase switch board matting online from http://switchboardmattingco.co.uk/

29. Matlab Homework Help -- 2018-01-16 07:36

This is really a great stuff for sharing. Keep it up .Thanks for sharing. https://www.matlabhelp.com

30. Eleston -- 2018-01-16 08:23

My first time to visit in the post. I see in the some math question seeing for here. I am appreciate in the work. keep it up. http://rubberflooringexpert.co.uk/

31. Programming Assignment Help -- 2018-01-16 11:29

This a good way to appreciate the teacher as they put their efforts to train students. UK dissertation Writers appreciates the teachers.<a https://programmingdoc.com/

32. Operations Management Assignment -- 2018-01-16 13:10

I genuinely appreciated understanding it. Sitting tight for some more incredible articles like this from you in the nearing days https://www.theoperationsmanagement.com/

33. Jahon -- 2018-01-25 11:28

I am very impress to visit in the post. I see in the post some mathematics question solve here. I am appreciate in the work. keep it up. http://industrialshelvinguk.co.uk/

34. anonymous -- 2018-01-26 06:56

Thank you very much sharing a very good information with us about Naive method i was looking to get information about the native method. Thanks a lot for it. To get safety flooring Click here http://safetyflooringco.co.uk/safety-flooring.html

35. anonymous -- 2018-01-26 08:45

Improving on Ogre3D is very good knowledge for me. Definitely this post enhance my knowledge. Thank you very much for sharing it and carry on it. To get playground tiles Click here http://www.rubberflooringuk.com/

36. anonymous -- 2018-02-07 02:38

Great article, very useful information, thanks for sharing

Jeff @ http://www.oakcitydrywall.com

37. Finance-Assignments.com -- 2018-02-08 10:45

https://financeassignments.xyz/ Hi buddy, your blog' s design is simple and clean and i like it. Your blog posts about Online Dissertation Help are superb. Please keep them coming. Greets!!

38. Buy Harvard Business Case Studies Solutions -- 2018-02-08 10:46

https://casehelp.xyz I appreciate your efforts in preparing this post. I really like your blog articles.

39. Kristian Zoppa -- 2018-02-14 17:31

I always had a good relationship with maths. Been a while since I've last seen vectors, but it's nevertheless a pleasure to find this. That's certainly some useful content.

40. Jack kevin -- 2018-03-14 12:15

The end-effector's revolution about its own z-pivot can be anything and does not make a difference for me. I have endeavored to ascertain the two vectors to a rotational framework from which I can without much of a stretch get the quaternions, however with each extraordinary technique for figuring the rotational network I get distinctive outcomes. If you need help in Assignment and want cheap assignment service so visit our website http://www.assignmentstar.com

41. sharonfrankklin@gmail.com -- 2018-03-17 08:09

Hi, great to see your website. I like the content and the research done behind every aspect of 70-346 questions answers

42. Kip Robinson -- 2018-03-24 16:39

I remember studying vectors in college, those were some fine days... I would like to find some window companies near me in Pensacola FL, do you know some trustworthy ones?

43. WORDAI Manual Multilevel Spinning -- 2018-03-28 11:27

only professional writers can make this kind of material, cheers https://www.contentastic.com/wordai-manual-article-spinning-12031

44. JCB Construction Equipment Made in India for the World case analysis -- 2018-03-28 11:28
45. Econometrics Assignment Help -- 2018-03-28 12:38

This a good way to appreciate the teacher as they put their efforts to train students. UK dissertation Writers appreciates the teachers. https://economicskey.com/econometrics-assignment-help-9236

46. FINANCIAL MANAGEMENT AN OVERVIEW Finance Assignment -- 2018-03-28 12:38

John arnold is an academic writer of the Dissertation-Guidance. Who writes quality academic papers for students to help them in accomplishing their goals. https://www.finance-assignments.com/financial-management-an-overview-3067

47. Nested switch Statements Java Assignment Help -- 2018-03-28 13:10

I personally like your post, you have shared good article. It will help me in great deal.https://javahelponline.com/nested-switch-statements-3080

48. Robust Control Homework Help -- 2018-03-28 13:11

This is really great work. Thank you for sharing such a useful information here in the blog. https://www.matlabhelp.com/robust-control-11131

49. Assembly Language Programming Homework Help -- 2018-03-28 13:33

Such a nice post, keep providing good resources. http://programmingdoc.com/assembly-language-10727

50. Stress Management Project Help -- 2018-03-28 13:34

Good way of telling, good post to take facts regarding my presentation subject matter, which i am going to deliver in my college http://www.medassignments.com/stress-management-9088

51. Karen -- 2018-04-08 02:06

Great post, thanks for sharing your amazing insight.

Karen at http://www.raleighawnings.com

52. Assignment Help -- 2018-05-07 07:42

We are the best in Assignment Help Online. Taking help our high expert's professionals who have extensive experience in writing assignment will not only make you avail with the quality assignment, but you will also get plagiarism free assignment within the given time. https://www.allassignmenthelp.com/

53. Olive.ethan101.co.uk@gmail.com -- 2018-05-07 12:08

Thanks for sharing it and keep on it. To get playground tiles. I am very galvanized to visit the post. Pay for Essay I see the publish some mathematics question solve right here. I admire the paintings.

54. AllAssignmentHelp reviews -- 2018-05-08 11:23

We Provide assignment help for students especially in usa getting brilliant quality reviews writing USA, essays and dissertations.We at Top Quality Assignment believe that there is no shortcut to success and to attain success, hard work, dedication, and commitment must be present.<a href="http://www.assignmentservicerating.com/a-client-review-on-allassignmenthelp-com/">AllAssignmentHelp.com reviews</a>  best in writing unique Assignment.

55. kristenwarne -- 2018-05-11 13:02

Get the best macroeconomics help services help from the experts of Students Assignment Help and at a low price. Our industry expert offers the reliable assignments services to the students. Our online expert writers are proficient in completing the assignment. For more info visit https://www.studentsassignmenthelp.com/macroeconomics-assignment-help/

56. mjjones733@gmail.com -- 2018-05-21 08:56

Wow! Happy to read such an informative writing piece. Keep doing good work. Allassignmenthelp is a site that is doing a fine job in providing online assignment help to the students. You can expect to get a fantastic assignment written by qualified writers using our service. https://www.allassignmenthelp.com/

57. anonymous -- 2018-06-07 21:55

Thank you for sharing!

http://www.arbuk.co.uk/

58. anonymous -- 2018-06-08 10:23

Fantastic blog post.

Kevin, http://kwikflow.com

59. oberoi.aishaa@gmail.com -- 2018-06-10 10:15

Commendable job with the article! It was really informative and enriching. Looking forward to more such posts in the future. Keep us updated with what’s in store! Will surely keep frequenting this website <a href=https://www.makemyassignments.com/>Make My Assignment </a>

60. anonymous -- 2018-06-15 13:09

Thank you for sharing this blog post! https://greentreesurgeonguildford.co.uk

61. anonymous -- 2018-06-20 06:38

All of this is so over my head, but I am trying to learn. Thanks for breaking it all down. http://www.brickmasonnashville.com

Add New Comment