Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix slerp and deprecate linpol #78

Merged
merged 18 commits into from
Apr 9, 2022
Merged

Conversation

hyrodium
Copy link
Collaborator

@hyrodium hyrodium commented Mar 4, 2022

This PR fixes #51.

Deprecate linpol

As discussed in #51 (comment), 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.
  • 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:

  • Add more tests
  • Check benchmark

@codecov
Copy link

codecov bot commented Mar 4, 2022

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%) ⬇️

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

@sethaxen
Copy link
Collaborator

sethaxen commented Mar 5, 2022

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
Copy link
Collaborator Author

hyrodium commented Mar 6, 2022

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
Copy link
Collaborator Author

hyrodium commented Mar 6, 2022

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
Copy link
Collaborator

sethaxen commented Mar 6, 2022

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
Copy link
Collaborator Author

hyrodium commented Mar 8, 2022

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
Copy link
Collaborator

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
Copy link
Collaborator Author

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.

image

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
Copy link
Collaborator Author

hyrodium commented Mar 21, 2022

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 hyrodium marked this pull request as ready for review March 31, 2022 12:34
@hyrodium
Copy link
Collaborator Author

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 hyrodium requested a review from sethaxen March 31, 2022 13:29
@hyrodium
Copy link
Collaborator Author

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

Copy link
Collaborator

@sethaxen sethaxen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Overall code looks good, and I agree on the scope of this PR. I made a few suggestions. Also, we need tests for the three added features:

  • the bug fix
  • norm=true
  • linpol is deprecated

@hyrodium
Copy link
Collaborator Author

hyrodium commented Apr 6, 2022

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
Copy link
Collaborator Author

hyrodium commented Apr 6, 2022

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
Copy link
Collaborator

sethaxen commented Apr 6, 2022

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
Copy link
Collaborator Author

hyrodium commented Apr 8, 2022

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
Copy link
Collaborator

sethaxen commented Apr 8, 2022

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

Did you push this commit?

@hyrodium
Copy link
Collaborator Author

hyrodium commented Apr 8, 2022

Did you push this commit?

Oh, I was forgetting. 😂

Co-authored-by: Seth Axen <seth.axen@gmail.com>
Copy link
Collaborator

@sethaxen sethaxen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM! I'll approve, but before merging, can you merge master into this PR and increment the patch number? Then after merging you can register a new release.

@hyrodium hyrodium merged commit 4cf7f2d into JuliaGeometry:master Apr 9, 2022
@hyrodium hyrodium mentioned this pull request Apr 29, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

slerp breaks flag for quaternion normalization
2 participants