Make rotations in degrees exact

[plut]

This PR adds a second type parameter to angular rotations type (Angle2d and all RotX etc.). This prevents automatic conversion of the angle to floating-point and the resulting inexactitude when angles are in degrees.

[codecov[bot]]

Codecov Report

Merging #197 (980cd77) into master (15d1c32) will decrease coverage by 0.28%. The diff coverage is 73.33%.

@@            Coverage Diff             @@
##           master     #197      +/-   ##
==========================================
- Coverage   86.12%   85.84%   -0.29%     
==========================================
  Files          14       14              
  Lines        1341     1349       +8     
==========================================
+ Hits         1155     1158       +3     
- Misses        186      191       +5     

Impacted Files	Coverage Δ
src/euler_types.jl	91.56% <66.66%> (-0.86%)	:arrow_down:
src/core_types.jl	89.70% <100.00%> (-0.60%)	:arrow_down:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 15d1c32...980cd77. Read the comment docs.

[hyrodium]

Great! I didn't know that the Unitful works exactly.

julia> import Unitful.°

julia> sin(30°)
0.5

julia> sin(π/6)
0.49999999999999994

[hyrodium]

Could you change AngleAxis too?

julia> using Rotations; import Unitful.°

julia> RotX(30°)
3×3 RotX{Float64, Quantity{Int64, NoDims, FreeUnits{(°,), NoDims, nothing}}} with indices SOneTo(3)×SOneTo(3)(30°):
 1.0  0.0        0.0
 0.0  0.866025  -0.5
 0.0  0.5        0.866025

julia> RotX(30°)[2,3]
-0.5

julia> AngleAxis(30°,1,0,0)
3×3 AngleAxis{Float64} with indices SOneTo(3)×SOneTo(3)(0.523599, 1.0, 0.0, 0.0):
 1.0  0.0        0.0
 0.0  0.866025  -0.5
 0.0  0.5        0.866025

julia> AngleAxis(30°,1,0,0)[2,3]
-0.49999999999999994

If AngleAxis is updated, #181 will be fixed.

[plut]

AngleAxis modified. Also, while I was looking in there, I was a bit surprised to see the details of this struct: wouldn't it have been a bit more natural to group x, y and z in a SVector{3,T}? (this doesn't change the memory layout, and of course it is possible to preserve all existing interfaces). This package already depends on StaticVectors anyway.

I'm in the process of writing extra tests.

[hyrodium]

Also, while I was looking in there, I was a bit surprised to see the details of this struct: wouldn't it have been a bit more natural to group x, y and z in a SVector{3,T}? (this doesn't change the memory layout, and of course it is possible to preserve all existing interfaces).

I agree with that. That was already implemented before I joined this package.

@andyferris Do you have any thoughts on this?

[plut]

Exactitude of coefficients is lost for AngleAxis, since the matrix is obtained through quaternions; this involves both angle-halving and L^2 normalization. In the general case, since the axis itself is floating-point, there are very few cases with exact coordinates: apart from RotX etc., only a few finite groups are involved (e.g. the cube isometry AngleAxis(60°,1,1,1)). (The existence of such cases however mean that it is sill an interesting goal to have exact formulas if possible!)

[hyrodium]

Exactitude of coefficients is lost for AngleAxis, since the matrix is obtained through quaternions; this involves both angle-halving and L^2 normalization.

I think we are not using quaternion here, just Rodrigues' rotation formula, (code in this repo).

In the general case, since the axis itself is floating-point, there are very few cases with exact coordinates: apart from RotX etc., only a few finite groups are involved

Exactly. At first, I thought it would be better to have AngleAxis{T,A} like RotX{T,A}, and AngleAxis(30°,1,0,0) can be computed exactly. But now, I'm suspicious to change AngleAxis... because there are only a few cases where the error is reduced, as you said.

We can choose from the following:

• Do not change AngleAxis
  • We can convert RotX to AngleAxis, but that will lose exactness.
• Change AngleAxis
  • I think most users don't care about the exactness, but it will be great if there's no performance degression.

(e.g. the cube isometry AngleAxis(60°,1,1,1)).

Maybe AngleAxis(120°,1,1,1)? :grin:

The existence of such cases however mean that it is sill an interesting goal to have exact formulas if possible!

Yeah, but it will be hard to express a point on a sphere with a floating-point number...

[plut]

I think we are not using quaternion here, just Rodrigues' rotation formula, (code in this repo).

But just three lines above your quote:

@inline Base.getindex(aa::AngleAxis, i::Int) = UnitQuaternion(aa)[i]

Rodrigues formula would indeed be easier for exactness (no angle halving, for a start).

Maybe AngleAxis(120°,1,1,1)? grin

Rational exactness is another useful property :grinning: (i.e. AngleAxis(60°, [1,1,1]//1) could return a rational matrix).

Yeah, but it will be hard to express a point on a sphere with a floating-point number...

There is this arXiv paper claiming a rational Rodrigues formula. I'm reading it right now; at first sight, it seems legit.

Edit: after reading the paper, this comes down to a simple fact: the rotation with axis vector u and angle θ has as coefficients rational functions of (tan(θ/2)/‖u‖) u (and with a denominator at most 1+tan²(θ/2)). For example, in the case of AngleAxis(120°, 1,1,1), since tan(θ/2) = √3/3 and ‖u‖ = √3, the coefficients are rational functions of u/3, and hence rational numbers.

Unfortunately, I believe that the simplest way to cover all those cases is to write ad-hoc code for the case when the angle is a multiple of 30° and the type of the axis vector is exact (integer or rational).

[hyrodium]

But just three lines above your quote:

Ah, I missed the line. That makes the following result.

julia> aa = rand(AngleAxis)
3×3 AngleAxis{Float64} with indices SOneTo(3)×SOneTo(3)(1.78335, 0.814399, 0.0957725, -0.572347):
  0.592207   0.653918  -0.470832
 -0.465016  -0.199847  -0.862451
 -0.658066   0.729693   0.185732

julia> RotMatrix(aa)[1]
0.5922065228804426

julia> aa[1]
0.5922065228804425

These values RotMatrix(aa)[1] and aa[1] should be equal. If we relace UnitQuaternion with Rodrigues formula, we will get more performance, I guess. (This should be fixed in other PR.)

i.e. AngleAxis(60°, [1,1,1]//1) could return a rational matrix

Understood, but I think it will be hard for type stability. For example, AngleAxis(59°, [1,1,1]//1) has same argument types, but the result cannot be rational.

Unfortunately, I believe that the simplest way to cover all those cases is to write ad-hoc code for the case when the angle is a multiple of 30° and the type of the axis vector is exact (integer or rational).

I agree with that. As mentioned above, obtaining exact rational numbers seems to be incompatible with type stability here.

Do you think adding the new type parameter A to AngleAxis{T,A} will be helpful? That may produce more precise results, but I think such cases are very rare. (e.g. AngleAxis(30°,1,0,0))

[plut]

i.e. AngleAxis(60°, [1,1,1]//1) could return a rational matrix

Understood, but I think it will be hard for type stability. For example, AngleAxis(59°, [1,1,1]//1) has same argument types, but the result cannot be rational.

While searching for the exact Rodrigues formula I found (but did not read) a few papers bout rational approximations for rotation matrices. While that's an interesting option for when the user explicitly asks for a rational answer, that is work for another PR indeed.

Do you think adding the new type parameter A to AngleAxis{T,A} will be helpful? That may produce more precise results, but I think such cases are very rare. (e.g. AngleAxis(30°,1,0,0))

For now I reverted AngleAxis to the original code AngleAxis{T}, since obtaining exact results here will be quite hard. Also, here is some more bad news:

julia> RotXYZ(30°,60°,90°)[3,2]
0.7499999999999999

This is due to the fact that sqrt(3)^2 ≠ 3. It might however be possible to make this exact by linearizing products of trigonometric functions (e.g. cos(a)*cos(b) = (cos(a+b)+cos(a-b))/2), but this calls for more evaluation of trigonometric functions and is hence more computationally expensive. I'm not sure a few exact coefficients are worth this cost. Mabe some RotXYZlinearized (etc.) types could (eventually) be implemented for covering those few cases however?

[hyrodium]

While searching for the exact Rodrigues formula I found (but did not read) a few papers bout rational approximations for rotation matrices. While that's an interesting option for when the user explicitly asks for a rational answer, that is work for another PR indeed.

:+1:

For now I reverted AngleAxis to the original code AngleAxis{T}, since obtaining exact results here will be quite hard.

:+1:

Also, here is some more bad news:

julia> RotXYZ(30°,60°,90°)[3,2]
0.7499999999999999

This is due to the fact that sqrt(3)^2 ≠ 3.

Ah, bad news. Btw, on the current master branch, it was able to calculate exactly, (luckily)

julia> RotXYZ(30°,60°,90°)
3×3 RotXYZ{Float64} with indices SOneTo(3)×SOneTo(3)(0.523599, 1.0472, 1.5708):
 3.06162e-17  -0.5        0.866025
 0.866025     -0.433013  -0.25
 0.5           0.75       0.433013

julia> RotXYZ(30°,60°,90°)[3,2]
0.75

Hmm, should we drop the accurate calculation of floating-point numbers here..?

but this calls for more evaluation of trigonometric functions and is hence more computationally expensive.

Mabe some RotXYZlinearized (etc.) types could (eventually) be implemented for covering those few cases however?

I think many users understand that floating-point numbers are approximate calculations, and adding RotXYZlinearized seems confusing. And perhaps the computation of a few trigonometric functions is not very costly.

[c42f]

I didn't know that the Unitful works exactly.

In this case Unitful is not doing anything very fancy, it's just causing sin(30°) to dispatch to sind(30) rather than sin(deg2rad(30)). And btw, sind(30) immediately uses floating point internally, so it would also be fine to keep the representation as a floating point number and have sind(30.0).

For now I reverted AngleAxis to the original code AngleAxis{T}, since obtaining exact results here will be quite hard.

This seems reasonable to me. I expect that code for exact rational rotations is never going to be good as a default codepath for a library which is mainly trying to do fast numerical (floating point) rotations. However, it could be useful to have alternative utility functions which try to go to great lengths to get a good rational answer or approximation.

An alternative viewpoint here could be that we accept floating point as a practical necessity, but try to make some of the conversions to rotation matrix elements correctly rounded. For example using compensated floating point arithmetic. This is the point of view that leads to functions like sind(): rather than trying to avoid floating point, they try to produce the closest floating point number to the exact answer.

Note that writing correctly rounding code which is fast is very hard though! Practically, it may be better to just use BigFloat in the cases where you care about accuracy and round it to the nearest Float64. This will be very slow, but should work immediately without any effort.

[plut]

I didn't know that the Unitful works exactly.

In this case Unitful is not doing anything very fancy, it's just causing sin(30°) to dispatch to sind(30) rather than sin(deg2rad(30)). And btw, sind(30) immediately uses floating point internally, so it would also be fine to keep the representation as a floating point number and have sind(30.0).

The point of the second type parameter is not keeping the angle an integer, it is keeping the angle in degrees, whereas the previous code lost exactness by converting the angle to radians. I expect the standard use case to look more like Angle2d(90.0°), where the angle will be stored as a float internally, which is fine since floats are exact for (reasonably-sized) integers (and for the same reason, the intent of the user is likely to mean the exact quarter-turn rotation).

[c42f]

The point of the second type parameter is not keeping the angle an integer, it is keeping the angle in degrees, whereas the previous code lost exactness by converting the angle to radians.

Sure, this is kind of what I was trying to point out: you want the computation to eventually use sind(x) rather than sin(deg2rad(x)). It doesn't matter whether x is a float or int (as floats exactly represent integers up to very large numbers) so the internal representation could always use a float if that's convenient.

[andyferris]

Do you have any thoughts on this?

The exact layouts of the structs isn't important to me. This package predates @c42f and I contributing (as well as StaticArrays), so I'm not sure exactly how old those choices are!

[KronosTheLate]

I just stumbled upon this, and though I should bump it. My 2 uninformed cents: It seems like perfect may be the enemy of the good here.