Slerp demystified

This method accurately replicates Vector3.SlerpUnclamed

Code (csharp):

1. // vectors need not be unit
2. static public Vector3 SLerp(Vector3 a, Vector3 b, float t) {
3.   float magOf(Vector3 v) => MathF.Sqrt(Vector3.Dot(v, v));
4.   float theta(Vector3 v0, Vector3 v1) => MathF.Atan2(magOf(Vector3.Cross(v0, v1)), Vector3.Dot(v0, v1));
5.   var it = 1f - t;
6.   var m0 = magOf(a);
7.   var m1 = magOf(b);
8.   var v0 = a / m0;
9.   var v1 = b / m1;
10.   var th = theta(a, b);
11.   var rs = 1f / MathF.Sin(th);
12.   var j = MathF.Sin(it * th);
13.   var k = MathF.Sin(t * th);
14.   return (m0 * it + m1 * t) * (rs * j * v0 + rs * k * v1);
15. }

Of course, Slerp on its own isn't actually a mystery, but the trick was to make it work the same as Unity's. Because it works with non-unit vectors I had to experiment. In the end, the solution was to slerp the units, and then apply lerp of the original magnitudes. However, what's funky is that the theta is computed using the originals.

Here, I've deliberately introduced

magOf

and

theta

methods to illustrate the amount of instructions needed. While

Atan2

itself is very fast, we also need 3

Sqrt

, 3

Sin

, 7 dots (cross counts as 3), and 7 divisions. So I've tried to simplify this a bit. This is what I came up with.

Code (csharp):

1. // vectors need not be unit
2. static public Vector3 SLerp(Vector3 a, Vector3 b, float t) {
3.   float magOf(Vector3 v) => MathF.Sqrt(Vector3.Dot(v, v));
4.   float theta(Vector3 v0, Vector3 v1) => MathF.Atan2(magOf(Vector3.Cross(v0, v1)), Vector3.Dot(v0, v1));
5.   var it = 1f - t;
6.   var m0 = magOf(a);
7.   var m1 = magOf(b);
8.   var th = theta(a, b);
9.   var s = MathF.Sin(th);
10.   var j = MathF.Sin(it * th);
11.   var k = MathF.Sin(t * th);
12.   return (m0 * it + m1 * t) * (j * m1 * a + k * m0 * b) / (s * m0 * m1);
13. }

3

Sqrt

, 7 dots, 1

Atan2

, 3

Sin

, 1 3 divisions (edit: oops it's still a vector scalar division, can't avoid that). I've managed to avoid pre-emptive normalization (because it got baked in the final division).

However, if we work with unit vectors instead, things simplify A LOT. We no longer need to find magnitudes, and we also don't need a magnitude of a cross, for example, because it's now 1. Now we have

Atan2(1, dot)

which simplifies to

Acos(dot)

. Finally, there's no need to lerp between the magnitudes anymore.

Code (csharp):

1. // vectors must be unit
2. static public Vector3 SLerpUnit(Vector3 a, Vector3 b, float t) {
3.   var th = MathF.Acos(Vector3.Dot(a, b));
4.   var s = MathF.Sin(th);
5.   var j = MathF.Sin((1f - t) * th);
6.   var k = MathF.Sin(t * th);
7.   return (j * a + k * b) / s;
8. }

This is significantly faster than the other one (1

Acos

, 3

Sin

, 1 dot, 1 3 divisions), however, each has its own use.

An incredible thing with this unit Slerp is that it always runs on the sphere's great circle, guaranteeing the shortest path in the sphere's surface subspace. Unit vectors in that case, simply define two end points on such an arc.

What if we could autodetect whether the two vectors were unit and pick the algorithm accordingly?
This is what I did -- but it does introduce a slight overhead to both methods. Obviously not the brightest idea, but I left it in for posterity.

Code (csharp):

1. // autodetects unit vectors
2. static public Vector3 SLerpAuto(Vector3 a, Vector3 b, float t) {
3.   bool isUnit(Vector3 v) => Math.Abs(Vector3.Dot(v, v) - 1f) <= 1E-4f;
4.   if(isUnit(a) && isUnit(b)) return SLerpUnit(a, b, t);
5.   return SLerp(a, b, t);
6. }

Here's the full component I used for testing.

Spoiler

Code (csharp):

1. using System;
2. using UnityEngine;
3.  
4. public class NewSLerpTest : MonoBehaviour {
5.  
6.   [SerializeField] Vector3 point1;
7.   [SerializeField] Vector3 point2;
8.   [SerializeField] [Range(0, 10)] int innerPoints;
9.   [SerializeField] SLerpMethod mode = SLerpMethod.Unity;
10.   [SerializeField] bool forceUnit = false;
11.   [SerializeField] [Min(.001f)] float gizmoScale = 1f;
12.  
13.   private enum SLerpMethod {
14.     Unity = 1,
15.     Custom = 2,
16.     Both = 3
17.   }
18.  
19.   void OnDrawGizmos() {
20.     var p1 = point1;
21.     var p2 = point2;
22.  
23.     if(forceUnit) {
24.       p1.Normalize();
25.       p2.Normalize();
26.     }
27.  
28.     drawPoint(p1, .1f, Color.white);
29.     drawPoint(p2, .1f, Color.yellow);
30.  
31.     for(int i = 0; i < innerPoints; i++) {
32.  
33.       float t = i / (float)innerPoints;
34.  
35.       Vector3 pt;
36.       if(((int)mode & 1) > 0) {
37.         pt = Vector3.SlerpUnclamped(p1, p2, t);
38.         drawPoint(pt, .06f, Color.cyan);
39.       }
40.  
41.       if(((int)mode & 2) > 0) {
42.         pt = SLerpAuto(p1, p2, t);
43.         drawPoint(pt, .06f, Color.magenta);
44.       }
45.  
46.     }
47.   }
48.  
49.   void drawPoint(Vector3 point, float radius = 1f, Color color = default) {
50.     if(color != default) Gizmos.color = color;
51.     Gizmos.DrawSphere(point, radius * gizmoScale);
52.   }
53.  
54.   // // vectors need not be unit (original unflattened)
55.   // static public Vector3 SLerp(Vector3 a, Vector3 b, float t) {
56.   //   float magOf(Vector3 v) => MathF.Sqrt(Vector3.Dot(v, v));
57.   //   float theta(Vector3 v0, Vector3 v1) => MathF.Atan2(magOf(Vector3.Cross(v0, v1)), Vector3.Dot(v0, v1));
58.   //   var it = 1f - t;
59.   //   var m0 = magOf(a);
60.   //   var m1 = magOf(b);
61.   //   var v0 = a / m0;
62.   //   var v1 = b / m1;
63.   //   var th = theta(a, b);
64.   //   var rs = 1f / MathF.Sin(th);
65.   //   var j = MathF.Sin(it * th);
66.   //   var k = MathF.Sin(t * th);
67.   //   return (m0 * it + m1 * t) * (rs * j * v0 + rs * k * v1);
68.   // }
69.  
70.   // vectors need not be unit
71.   static public Vector3 SLerp(Vector3 a, Vector3 b, float t) {
72.     float magOf(Vector3 v) => MathF.Sqrt(Vector3.Dot(v, v));
73.     float theta(Vector3 v0, Vector3 v1) => MathF.Atan2(magOf(Vector3.Cross(v0, v1)), Vector3.Dot(v0, v1));
74.     var it = 1f - t;
75.     var m0 = magOf(a);
76.     var m1 = magOf(b);
77.     var th = theta(a, b);
78.     var s = MathF.Sin(th);
79.     var j = MathF.Sin(it * th);
80.     var k = MathF.Sin(t * th);
81.     return (m0 * it + m1 * t) * (j * m1 * a + k * m0 * b) / (s * m0 * m1);
82.   }
83.  
84.   // vectors must be unit
85.   static public Vector3 SLerpUnit(Vector3 a, Vector3 b, float t) {
86.     var th = MathF.Acos(Vector3.Dot(a, b));
87.     var s = MathF.Sin(th);
88.     var j = MathF.Sin((1f - t) * th);
89.     var k = MathF.Sin(t * th);
90.     return (j * a + k * b) / s;
91.   }
92.  
93.   // autodetects unit vectors
94.   static public Vector3 SLerpAuto(Vector3 a, Vector3 b, float t) {
95.     bool isUnit(Vector3 v) => Math.Abs(Vector3.Dot(v, v) - 1f) <= 1E-4f;
96.     if(isUnit(a) && isUnit(b)) return SLerpUnit(a, b, t);
97.     return SLerp(a, b, t);
98.   }
99.  
100. }

Haven't benchmarked any of this, performance-wise, because i'm lazy. I would be surprised if the fully-blown one matched the Unity's own SlerpUnclamped, but I believe the unit Slerp is faster (just an intuition without any evidence to back it up). Any feedback is appreciated.

 

@MelvMay Thank you!
For the record this is Quaternion.Slerp. There are some key differences, namely this is always unit Slerp.

 

"Why don't you try to replicate Slerp?" they said
"It'll be fun" they said
lol

Anyway, I managed to fix a lot of edge-case issues, including the antiparallel case.
@MelvMay thanks a lot for the insight, that code was really useful.

This is the final version, for the moment. I haven't been able to find a discrepancy for a while now, and this time I was hunting for all kinds of contrived situations ranging from rotation ambiguity to strict orthogonality, including cases like having point 1 at (10,10,10) and point 2 at (-1,-1,-1), for example. Earlier today this gave me a beautiful symmetric pattern but alas I had to hate it because it was so wrong.

Code (csharp):

1. // vectors need not be unit
2. static public Vector3 SLerp(Vector3 a, Vector3 b, float t) {
3.   float magOf(Vector3 v) => MathF.Sqrt(Vector3.Dot(v, v));
4.   var it = 1f - t; // inverse of t
5.   var m0 = magOf(a);
6.   var m1 = magOf(b);
7.   var mm = m0 * m1; // combined magnitude
8.   if(mm == 0f) return a * it + b * t; // use lerp if one of the vectors is zero
9.   var d = Vector3.Dot(a, b) / mm; // unit dot
10.   if(1f - Math.Abs(d) <= 1E-5f) { // abs(dot) is close to 1
11.     return d > 0f? a * it + b * t // << use lerp for very small angles
12.                  : AxisAngle(FastOrthogonal(a), MathF.PI * t) * ((it + (m1/m0) * t) * a);
13.   }           // ^^ vectors are antiparallel, apply rotation on orthogonal axis, lerp mag
14.   var th = MathF.Acos(d);
15.   var s = MathF.Sin(th) * mm;
16.   var j = MathF.Sin(it * th);
17.   var k = MathF.Sin(t * th);
18.   // left-hand-side scalar part = mag lerp
19.   // right-hand-side vector part = actual slerp
20.   return (m0 * it + m1 * t) * (j * m1 * a + k * m0 * b) / s;
21. }

Anyway it turned out that I needed a much more robust orthonormal function, so I finally got to make one.

Code (csharp):

1. static public Vector3 FastOrthogonal(Vector3 v, bool normalize = true) {
2.   var sqr = v.x * v.x + v.y * v.y;
3.   if(sqr > 0f) { // (0,0,1) x (x,y,z)
4.     var im = normalize? 1f / MathF.Sqrt(sqr) : 1f;
5.     return new Vector3(-v.y * im, v.x * im, 0f);
6.   } else {       // (1,0,0) x (x,y,z)
7.     sqr = v.y * v.y + v.z * v.z;
8.     var im = normalize? 1f / MathF.Sqrt(sqr) : 1f;
9.     return new Vector3(0f, -v.z * im, v.y * im);
10.   }
11. }

This is effectively a hacked cross product, which optionally normalizes the result (aka returns orthonormal). We can also do

FastOrthogonal(a/m0, false)

but in any case

axis

has to be unit, because quaternion expects a normal.

Also the antiparallel solution is a weird one, and here I'm hijacking quaternions for the final multiplication, but I bypass normalization (and a conversion from degrees to radians). AxisAngle is thus

Code (csharp):

1. // axis must be unit, rad is angle in radians
2. static public Quaternion AxisAngle(Vector3 axis, float rad) {
3.   rad *= .5f;
4.   axis *= MathF.Sin(rad);
5.   var q = new Quaternion(axis.x, axis.y, axis.z, MathF.Cos(rad));
6.   //Debug.Assert(q == q.normalized);
7.   return q;
8. }

I couldn't notice any scenario in which this quaternion wasn't unit (if the axis was unit), so I skipped that part. I will probably stress-test this some more and analyze the math in more detail. I'm too lightheaded these days for any kind of complex algebra.

Finally, because I rotate vector a directly, mag lerp is a bit different than usual, basically I just scale it from 1 towards b/a. Hence

Code (csharp):

1. AxisAngle(FastOrthogonal(a), MathF.PI * t) * ((it + (m1/m0) * t) * a);

A much simpler SLerpUnit now looks like this

Code (csharp):

1. // vectors must be unit
2. static public Vector3 SLerpUnit(Vector3 a, Vector3 b, float t) {
3.   var d = Vector3.Dot(a, b);
4.   if(1f - Math.Abs(d) <= 1E-6f) // smaller epsilon
5.     return d > 0f? a * (1f - t) + b * t
6.                  : AxisAngle(FastOrthogonal(a, normalize: false), MathF.PI * t) * a;
7.   var th = MathF.Acos(d);
8.   var s = MathF.Sin(th);
9.   var j = MathF.Sin((1f - t) * th);
10.   var k = MathF.Sin(t * th);
11.   return (j * a + k * b) / s;
12. }