Grouping Rotation objects in a vector

[baptnz]

I’m working with 3D rotations, and my current code implements rotation matrices from Euler angles.

I would like to switch to a more general approach that leverages Rotations.jl (or ReferenceFrameRotations.jl, but that’s a separate discussion). I like the idea of abstracting away the specific convention that the user may prefer (parametrising with Euler angles, or axis-angle, or quaternions, etc.) but I haven’t been able to adapt my code to work with the supertype common to them all. Here’s what I mean:

using StaticArrays

# current definition: each object is given 3 Euler angles to describe the rotation
struct current_clust{T}

    angles::Vector{SVector{3,T}}

end

current_clust(SVector.(1:4,2,3))
# current_clust{Int64}(SVector{3, Int64}[[1, 2, 3], [2, 2, 3], [3, 2, 3], [4, 2, 3]])

vs my attempt at a new approach:

using Rotations

# new implementation: each object would be assigned a generic "Rotation" object
struct new_clust{T}

    angles::Vector{Rotation{3,T}}

end

r = rand(RotMatrix{3})
q = UnitQuaternion(r)
e = RotXYZ(0,1,2)
aa = AngleAxis(r)

new_clust([r, q, e, aa])
# ERROR: DimensionMismatch("No precise constructor for Rotation{3, Float64} found. Length of input was 9.")

Is my idea flawed (e.g. would the new_clust objects be inefficient because of arbitrary types)? Or do I just have the syntax wrong? (e.g I need a different container)?

I don’t really know how to probe the different subtypes – from what is printed they all seem “compatible”, as in

supertype(typeof(r))
Rotation{3, Float64}

etc.

giving me the impression that I could just group them together in a vector. But at the same time I expect that internally they’re stored quite differently (e.g. 3, 4, or maybe 9 elements, etc.). How can I check this kind of thing?

Of course I could coerce everything into 3x3 matrices, but that would seem a bit wasteful.

[rdeits]

@baptnz

e.g. would the new_clust objects be inefficient because of arbitrary types

Your instinct here is correct–Rotation is an abstract type, so you’d be running into this performance tip: Performance Tips · The Julia Language

If it were me, I would probably skip the generic container entirely and store a Vector{UnitQuaternion{T}}. Dealing with collections of differently-typed objects can be tricky, and it’s not actually obvious that you need that here. Since a quaternion is already a nice representation of a rotation, it might make sense to just pick that fundamental type to store.

Note that this doesn’t restrict you at all from populating that vector with whatever type is convenient. You can push! to that vector and any element you add will automatically try to convert into a quaternion:

julia> using Rotations

julia> vec = Vector{UnitQuaternion{Float64}}()
UnitQuaternion{Float64}[]

julia> push!(vec, AngleAxis(0, 1, 0, 0))
1-element Vector{UnitQuaternion{Float64}}:
 [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]

julia> push!(vec, RotXYZ(0, 1, 2))
2-element Vector{UnitQuaternion{Float64}}:
 [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
 [-0.22484509536615288 -0.49129549643388193 0.8414709848078966; 0.9092974268256818 -0.41614683654714246 -5.551115123125783e-17; 0.3501754883740148 0.7651474012342927 0.5403023058681398]

julia> vec
2-element Vector{UnitQuaternion{Float64}}:
 [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
 [-0.22484509536615288 -0.49129549643388193 0.8414709848078966; 0.9092974268256818 -0.41614683654714246 -5.551115123125783e-17; 0.3501754883740148 0.7651474012342927 0.5403023058681398]

The only reason to actually store the separate rotations as different types is if you actually care about accessing them in their original types later, rather than just as some totally equivalent representation.

[baptnz]

Thanks, I think I agree with your points but just two quick clarifications:

Since a quaternion is already a nice representation of a rotation, it might make sense to just pick that fundamental type to store.

Picking just one type makes sense to me, but I’m not sure this is the best for my use-case. I worry quaternions might get in the way of ForwardDiff (I couldn’t quite make sense of some thread discussing the matter), and they’re also the least-intuitive (to me at least). 3x3 matrices feel like one obvious alternative, but they’re quite wasteful. Euler angles have their own problems (gimbal lock etc.). Maybe I should go with axis-angle; they’re not friendly to combine, but as you say I could always use different formulations beforehand, and coerce into this specific form for convenient storage.

The only reason to actually store the separate rotations as different types is if you actually care about accessing them in their original types later, rather than just as some totally equivalent representation.

Yes, you make a good point: I was thinking of convenience to generate those rotations (using whichever convention is easier for a given problem, possibly combining various rotations into one, etc.) but those steps don’t need to be tied to the final storage mode: I can just coerce into one format (provided it’s a good one…).

[rdeits]

FWIW I’ve done lots of autodiff computation through quaternions, and it certainly can work just fine. You may run into edge cases, but that’s true with any representation, so you might as well pick the one that isn’t a minefield of singularities.

@baptnz

Yes, you make a good point: I was thinking of convenience to generate those rotations (using whichever convention is easier for a given problem, possibly combining various rotations into one, etc.) but those steps don’t need to be tied to the final storage mode

Right–that all sounds great, and would make lots of sense as one or more helpful methods for your new_clust{T} type.

[baptnz]

@rdeits

FWIW I’ve done lots of autodiff computation through quaternions, and it certainly can work just fine.

That’s reassuring, thanks! I’ll give them a try then.

[Orbots]

@baptnz

they’re also the least-intuitive (to me at least)

Quaternions are actually quite intuitive if you first learn what a bivector is. No hyperspheres needed. Then you also get a very nice connection between angle-axis and quaternions via exponential/log maps as a bonus. exp(angle_axis*i) = quaternion.

Any intro to geometric algebra will have a section on quaternions ( they are called rotors in GA ).

This here is pretty easy to follow. Note that quaternions are really the same thing as rotors ( dual ).

https://marctenbosch.com/quaternions/

[baptnz]

Thanks – I’ve had geometric algebra books on my reading list for years.

[baptnz]

For info, I’ve now refactored my code to use Rotations.jl and store the parameters as quaternions. Pretty happy with the change so far, it’s quite convenient to be able to choose the easiest parametrisation for a given geometry. And quaternions do come in handy for composition. I have yet to try ForwardDiff on the code though.

Amusingly, three.js quaternions are stored with a different convention, which tripped me at first.

[Orbots]

@baptnz

Amusingly, three.js quaternions are stored with a different convention, which tripped me at first.

Yeah, you need to read the code when looking at quaternion packages, since conventions differ. There has been some interest in unifying quaternions here. Taking Quaternions Seriously