Fix `slerp` and deprecate `linpol`

[hyrodium]

This PR fixes #51.

Deprecate linpol

As discussed in https://github.com/JuliaGeometry/Quaternions.jl/issues/51#issuecomment-1046887775, slerp seems a more common name.

slerp now normalizes input quaternions

We had the following choices, and I choose the first one:

• Normalize all input quaternions, and the return value is always unit quaternion. (like linpol)
  • linpol is deprecated for its name, and the original slerp for non-unit quaternions seems buggy, so I thought it's better to have the same behavior as linpol.
• Throws DomainError for non-unit quaternion.
  • We will drop norm field with #75, and checking
• Return exp(log(q1)*(1-t) + log(q2)*t) for any quaternions.
  • Slerp is shorthand for spherical linear interpolation, so it seems better to avoid non-unit quaternions output.

The returned Quaternion has norm flag true.

The flag will be removed with #75, but seems better to fix the flag for now.

TODO:

• [x] Add more tests
• [x] Check benchmark

[codecov[bot]]

Codecov Report

Merging #78 (80887d7) into master (1f25259) will decrease coverage by 0.06%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #78      +/-   ##
==========================================
- Coverage   97.27%   97.21%   -0.07%     
==========================================
  Files           3        3              
  Lines         404      395       -9     
==========================================
- Hits          393      384       -9     
  Misses         11       11              

Impacted Files	Coverage Δ
src/Quaternion.jl	98.50% <100.00%> (-0.07%)	:arrow_down:

:mega: Codecov can now indicate which changes are the most critical in Pull Requests. Learn more

[sethaxen]

Return exp(log(q1)*(1-t) + log(q2)*t) for any quaternions.

* Slerp is shorthand for spherical linear interpolation, so it seems better to avoid non-unit quaternions output.

I actually quite like this option. Although, I'll note that this expression is not quite correct. If q1 and q2 are on the same half-sphere, then slerp(q1, q2, t) is equivalent to (q2 / q1)^t * q1. If q1 = a1 u1 and q2 = a2 u2, where u1 and u2 are unit quaternions and a1>0 and a2>0, then (q2 / q1)^t * q1 = a * slerp(u1, u2, t), where log(a) = (1-t)*log(a1) + t * log(a2). i.e. this extension of slerp still spherically interpolates between the direction of the quaternions while it logarithmically interpolates between magnitudes, which seems sensible. I suspect it's more efficient still to normalize first and use slerp as implemented here, but the extension could be preserved by a single scalar multiplication at the end.

[hyrodium]

Thanks for the correction! The exp(log(q1)*(1-t) + log(q2)*t) was my mistake. However,

(q2 / q1)^t * q1 = a * slerp(u1, u2, t), where log(a) = (1-t)*log(a1) + t * log(a2)

is also not quite correct with current slerp (and linpol) implementation.

julia> using Quaternions

julia> q1 = Quaternion(1.,0,0,0)
QuaternionF64(1.0, 0.0, 0.0, 0.0, false)

julia> q2 = Quaternion(-1/√2,-1/√2,0,0)
QuaternionF64(-0.7071067811865475, -0.7071067811865475, 0.0, 0.0, false)

julia> t = 1.0
1.0

julia> slerp(q1, q2, t)
QuaternionF64(0.7071067811865475, 0.7071067811865475, 0.0, 0.0, false)

julia> linpol(q1, q2, t)
QuaternionF64(0.7071067811865476, 0.7071067811865476, 0.0, 0.0, true)

julia> q1*(q2/q1)^t
QuaternionF64(-0.7071067811865474, -0.7071067811865475, 0.0, 0.0, false)

The current slerp and linpol equate antipodal points, but I think this sometimes does not seem to be sensible. Can I add a flag shortest_path=false argument in slerp?

• http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/slerp/

  • Since (w,x,y,z) and (-w,-x,-y,-z) represent the same rotation we should make sure that are not sensitive to whether the positive of inverted form of the quaternion are used.

  • It might be best if this modification is optional?

• https://en.wikipedia.org/wiki/Slerp

  • However, because the covering is double (q and −q map to the same rotation), the rotation path may turn either the "short way" (less than 180°) or the "long way" (more than 180°). Long paths can be prevented by negating one end if the dot product, cos Ω, is negative, thus ensuring that −90° ≤ Ω ≤ 90°.

[hyrodium]

If we have the antipodal identification, then the magnitude interpolation might be useless, I guess. If someone needs the interpolation, one can just evaluate q1*(q2/q1)^t instead of slerp.

[sethaxen]

The current slerp and linpol equate antipodal points, but I think this sometimes does not seem to be sensible.

When would this not be sensible? The unit quaternions doubly cover the rotations, so that q and -q correspond to the same rotation. While the unit quaternions form a 3-sphere, the rotations form a 3-hemisphere. The point of slerp is to smoothly interpolate between 2 rotations/orientations. Specifically, it computes the shortest path (geodesic) between the two rotations and then traverses the specified distance on that path. One could compute a geodesic on the unit quaternions, but this wouldn't be slerp, and it wouldn't be all that useful.

If we have the antipodal identification, then the magnitude interpolation might be useless, I guess.

I don't think so. The non-zero quaternions can be interpreted as the direct product of two manifolds: the (right-)quaternionic projective space (the hemisphere of unit quaternions corresponding to rotations) and the scale group (strictly positive reals). This is the manifold on which the slerp I'm proposing would be computing shortest paths.

If someone needs the interpolation, one can just evaluate q1*(q2/q1)^t instead of slerp.

When I benchmarked this, I found it to be quite inefficient compared to slerp. Probably because it currently computes log and exp without any simplification.

I suppose I haven't argued for the benefits of doing it this way. Here are a few:

1. It's more general. I'd like to minimize how much rotation-specific code we have here, since quaternions are useful and interesting in their own right, not just as representations of rotations. If code has generalizations to the non-unit quaternions without much decrease in efficiency, it would be nice to implement the generalization.
2. It reduces sensitivity to normalization. A user doesn't need to worry about whether their quaternion is normalized before or after computing slerp. They can just normalize it whenever they need it to be normalized (which is probably right before computing a rotation matrix or rotating a pure quaternion). And, there might be a way to implement slerp this way that never requires computing the norm, and is therefore more efficient.
3. A person might want to use a quaternion to store both rotation and scale information. In this case, the generalization of slerp would just do the right thing.

I don't think this is super important, and if you feel very strongly that slerp should only emit unit quaternions, then we can go that way.

[hyrodium]

When would this not be sensible?

One could compute a geodesic on the unit quaternions, but this wouldn't be slerp, and it wouldn't be all that useful.

The sphere S³ ≃ SU(2) ≃ (set of unit quaternions) is not just a double cover of SO(3) but also a universal cover. This means each point on S³ can be regarded as a path (or equivalence class of paths) in SO(3).

These videos might be helpful to understand the coverings: https://vimeo.com/62228139 http://www.gregegan.net/APPLETS/21/21.html https://youtu.be/8mhvnXWzlHw?t=146

Therefore, one quaternion can be interpreted as a path of rotation, so I think it would be useful to distinguish q and -q (with an optional argument like shortest_path).

The non-zero quaternions can be interpreted as the direct product of two manifolds

A person might want to use a quaternion to store both rotation and scale information.

I was thinking that unit quaternions is just a subset of ℍ, but I now agree with the interpolation of magnitude! (ℍ-{0} ≃ ℝ₊ × S³) However, that will be a breaking change with linpol, so I'd like to keep the normalization in this PR.

[sethaxen]

The sphere S³ ≃ SU(2) ≃ (set of unit quaternions) is not just a double cover of SO(3) but also a universal cover. This means each point on S³ can be regarded as a path (or equivalence class of paths) in SO(3).

I understand S³ universally covers SO(3). I do not understand the assertion that each point on S³ is a path on SO(3), and I did not see this assertion made in any of the provided links.

I also don't understand how this relates to the points I made above about what slerp is useful for. Are you aware of any descriptions of slerp that explicitly draw paths between quaternions from different hemispheres (as opposed to omitting the identification of antipodal quaternions for simplicity of exposition)? I understand when considering the unit quaternions as a manifold (or considering any manifold for that matter), this is something one might want to do, but slerp is concerned with 3D rotations, hence why it interpolates on a single hemisphere.

I was thinking that unit quaternions is just a subset of ℍ, but I now agree with the interpolation of magnitude! (ℍ-{0} ≃ ℝ₊ × S³) However, that will be a breaking change with linpol, so I'd like to keep the normalization in this PR.

Ah, yes, then this is the better way to go for now.

[hyrodium]

Sorry for the late reply.

I do not understand the assertion that each point on S³ is a path on SO(3)

A universal covering manifold M̃ of a manifold M can be constructed by the following quotient set.

See section 1.5 (d) in Shigeyuki Morita, Geometry of Differential Forms for more infomation.

I did not see this assertion made in any of the provided links.

Sorry for the insufficient explanation. I was trying to explain that the 0° rotation and 360° rotation are different, but 0° and 720° are equivalent because a path from 0° to 720° is contractible.

One could compute a geodesic on the unit quaternions, but this wouldn't be slerp, and it wouldn't be all that useful.

Are you aware of any descriptions of slerp that explicitly draw paths between quaternions from different hemispheres (as opposed to omitting the identification of antipodal quaternions for simplicity of exposition)?

If a robot has a arm like human body, then the rotation 90° and -270° should be distinguishable because the path is not contractible, and one might need different slerps slerp(0°, 90°, t) and slerp(0°, -270°, t). I am not an expert in the applications and have yet to find any resources for the slerp.

[hyrodium]

To be clarified, I made two misunderstandings comments:

• exp(log(q1)*(1-t) + log(q2)*t): this should be q1*(q2/q1)^t.
• I was thinking that the slerp should distinguish antipodal points q and -q.

I now think that it would be better to have an optional argument like shortest_path or geodesic_path. But this doesn't need to be in this PR.

[hyrodium]

I merged the master branch, and I think this is ready to merge.

What I have done in this PR:

Fix bug in slerp

Current master

julia> q3 = Quaternion(2.,1.,0.,0.,true)
QuaternionF64(2.0, 1.0, 0.0, 0.0, true)

julia> q4 = Quaternion(3.,0.,0.,0.,true)
QuaternionF64(3.0, 0.0, 0.0, 0.0, true)

julia> slerp(q3, q4, 0.2)
QuaternionF64(2.0, 1.0, 0.0, 0.0, true)

With this PR

julia> q3 = Quaternion(2.,1.,0.,0.,true)
QuaternionF64(2.0, 1.0, 0.0, 0.0, true)

julia> q4 = Quaternion(3.,0.,0.,0.,true)
QuaternionF64(3.0, 0.0, 0.0, 0.0, true)

julia> slerp(q3, q4, 0.2)
QuaternionF64(0.9319949582311067, 0.3624712372475889, 0.0, 0.0, true)

Fix norm flag

Current master

julia> q1 = Quaternion(1.0, 0.0, 0.0, 0.0, true)
QuaternionF64(1.0, 0.0, 0.0, 0.0, true)

julia> q2 = Quaternion(0.0, 1.0, 0.0, 0.0, true)
QuaternionF64(0.0, 1.0, 0.0, 0.0, true)

julia> slerp(q1,q2,0.3)
QuaternionF64(0.8910065241883678, 0.45399049973954675, 0.0, 0.0, false)

With this PR

julia> q1 = Quaternion(1.0, 0.0, 0.0, 0.0, true)
QuaternionF64(1.0, 0.0, 0.0, 0.0, true)

julia> q2 = Quaternion(0.0, 1.0, 0.0, 0.0, true)
QuaternionF64(0.0, 1.0, 0.0, 0.0, true)

julia> slerp(q1,q2,0.3)
QuaternionF64(0.8910065241883678, 0.45399049973954675, 0.0, 0.0, true)

Deprecate linpol

slerp seems more common wording than linpol.

Not in this PR:

• Interpolation for quaternion magnitude
• Adding optional arguments for distinguishing antipodal points q and -q.

[hyrodium]

I've also added DomainError for zero-quaternion.

Current master

julia> slerp(quat(1),quat(0),1)
QuaternionF64(0.0, 0.0, 0.0, 0.0, false)

With this PR

julia> slerp(quat(1),quat(0),1)
ERROR: DomainError with Quaternion{Int64}(0, 0, 0, 0, false):
The input quaternion must be non-zero.
Stacktrace:
 [1] slerp(qa::Quaternion{Int64}, qb::Quaternion{Int64}, t::Int64)
   @ Quaternions ~/.julia/dev/Quaternions/src/Quaternion.jl:312
 [2] top-level scope
   @ REPL[2]:1

[hyrodium]

I've added the following tests.

• the bug fix: @testset "Normalizing input quaternions"
• norm=true: @test slerp(a, b, 0.0).norm etc.
• linpol is deprecated: @testset "Deprecated warning"

[hyrodium]

Ohh, I found some performance degressions.

current master

julia> using Quaternions, BenchmarkTools

julia> q1,q2 = normalize.(rand(QuaternionF64,2))
2-element Vector{QuaternionF64}:
 QuaternionF64(0.13376212455218736, 0.8893671325737943, 0.2937038326951311, 0.32383924436816264, true)
  QuaternionF64(0.4630216872108008, 0.4083229516075815, 0.7029702168934131, 0.35315174999391385, true)

julia> @benchmark slerp($q1,$q2,0.2)
BenchmarkTools.Trial: 10000 samples with 994 evaluations.
 Range (min … max):  33.141 ns … 57.069 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     34.058 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   34.450 ns ±  1.993 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▁▂█  ▃                                                      
  ▃███▄▇█▅▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▁▁▁▁▁▂▂▂▂▁▂▂▂▁▂▁▂▂▂ ▃
  33.1 ns         Histogram: frequency by time        46.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark linpol($q1,$q2,0.2)
BenchmarkTools.Trial: 10000 samples with 989 evaluations.
 Range (min … max):  46.063 ns … 62.636 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     47.106 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   47.183 ns ±  0.779 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▁     █  ▅▂                                           
  ▃▁▂▃▁▁▄█▁▆▇▁▃█▂▂██▂▂▇▁▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  46.1 ns         Histogram: frequency by time        50.6 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

with this PR

julia> using Quaternions, BenchmarkTools

julia> q1,q2 = normalize.(rand(QuaternionF64,2))
2-element Vector{QuaternionF64}:
 QuaternionF64(0.26545001572497234, 0.31425231331646797, 0.6605279482815369, 0.6280800922381187, true)
 QuaternionF64(0.21867092776216737, 0.6423537858021358, 0.6971969708881615, 0.23125964412683261, true)

julia> @benchmark slerp($q1,$q2,0.2)
BenchmarkTools.Trial: 10000 samples with 981 evaluations.
 Range (min … max):  61.737 ns … 175.857 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     62.493 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   63.211 ns ±   2.643 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅   █ ▄                                  ▃                   ▁
  █▂▆▃█▃█▅▅▅▅▅▅▅▅▆▅▆▆▆▆▆▆▇▆▇▇▇▆▆▅▅▅▄▅▄█▆█▅▄█▅█▆▅▆▅▆▅▅▄▄▄▅▄▄▄▅▅ █
  61.7 ns       Histogram: log(frequency) by time      72.3 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark linpol($q1,$q2,0.2)

BenchmarkTools.Trial: 10000 samples with 964 evaluations.
 Range (min … max):  83.602 ns … 132.761 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     85.369 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   85.623 ns ±   1.951 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁▄▄▃▃▄▅██▂       ▁                                           ▂
  ███████████▇▆▆▆██████▇█▇████▇▅▄▄▄▄▄▅▄▄▄▅▁▄▁▄▃▁▁▄▁▄▁▁▁▄▃▁▁▆▄▇ █
  83.6 ns       Histogram: log(frequency) by time      97.1 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[sethaxen]

Can the difference be explained by the normalization? On my machine the performance regression is less extreme, but adding the sign computations on master brought them to about the same runtime.

[hyrodium]

Yeah, the problem was around normalization with sign. I replaced sign(q) with q/abs(q) because we have iszero check. Now we have the following benchmark result.

julia> using Quaternions, BenchmarkTools

julia> q1,q2 = normalize.(rand(QuaternionF64,2))
2-element Vector{QuaternionF64}:
    QuaternionF64(0.495912700698677, 0.47929964850623047, 0.5228646779504917, 0.5009540585515571, true)
 QuaternionF64(0.7828859145095284, 0.27796675430678036, 0.49666698530615594, 0.25128874640466514, true)

julia> @benchmark slerp($q1,$q2,0.2)
BenchmarkTools.Trial: 10000 samples with 991 evaluations.
 Range (min … max):  40.218 ns … 63.298 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     41.804 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   42.091 ns ±  1.616 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▃    █▂     ▂                                             
  ▂▂▅█▄▅█▁██▂▂▃▂▆█▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▁▂▁▂▁▁▂▂▁▂▂▁▂▂▂ ▃
  40.2 ns         Histogram: frequency by time        51.8 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark linpol($q1,$q2,0.2)
BenchmarkTools.Trial: 10000 samples with 974 evaluations.
 Range (min … max):  71.233 ns … 166.682 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     72.098 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   72.533 ns ±   2.557 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▂ ▇ █ ▅ ▃  ▅ ▁                                               ▂
  ███▁█▅█▅█▅▄█▇█▇█▇▇▆▆▆▆▇▆▆▆▇▇▅▆▅▅▅▄▅▄▃▅▁▄▁▁▁▄▁▄▁▁▁▁▃▃▁█▃█▇▁▅▇ █
  71.2 ns       Histogram: log(frequency) by time      82.7 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

I'm not sure why linpol takes more time than slerp, but we don't need to care about that because it is deprecated now.

[sethaxen]

I replaced sign(q) with q/abs(q)

Did you push this commit?

[hyrodium]

Did you push this commit?

Oh, I was forgetting. :joy: