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

Docs tutorial, "How to work with Rotations" #355

Closed
wants to merge 12 commits into from
5 changes: 4 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,10 @@ makedocs(
pages=[
"Home" => "index.md",
"ManifoldsBase.jl" => "interface.md",
"Examples" => ["How to implement a Manifold" => "examples/manifold.md"],
"Examples" => [
"How to implement a Manifold" => "examples/manifold.md",
"How to work with Rotations" => "examples/dimetourrotations.md",
],
"Manifolds" => [
"Basic manifolds" => [
"Centered matrices" => "manifolds/centeredmatrices.md",
Expand Down
161 changes: 161 additions & 0 deletions docs/src/examples/dimetourrotations.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
# [How to work with Rotations](@id how_to_work_with_rotations)

This tutorial is meant to give the briefest of overviews on how to use Manifolds.jl in a manner familiar to those needing rigid body transforms. This tutorial will introduce some of the common function calls needed to convert between the various data types and hopefully show the user something more about to combine some of the function calls listed elsewhere in the documentation.

## Rotations with SO(2)

Consider rotations on an xy-plane, commonly known (among others) as rotation matrices `R`, Direction Cosine Matrices `DCM`, `SpecialOrthogonal(2)` Lie Groups and associated Lie Algebra. Let's load some packages first:
```julia
using Manifolds
using LinearAlgebra
using StaticArrays
```

### Manifolds and Defaults

The associated manifolds and groups are defined by:
```julia
G = SpecialOrthogonal(2)
M = base_manifold(G)
@assert M == Rotations(2)
```

Pretty soon we will require some default definitions:
```julia
# default basis
B = DefaultOrthogonalBasis()
# default data type
p0 = @SMatrix [1.0 0; 0 1]

# Group identity element with zero aligned to the x-axis
xR0 = identity(G, p0)
```

!!! note
Throughout the Manifolds code you will likely find the point `p` on the manifold with `X` a tangent vector at `p` (for example on the sphere `p=[1.0,0.0,0.0]` and `X=[0.0,0.5,0.5]`, and notice that points on the sphere are represented as unit vectors and all tangent space vectors at `p` are orthogonal to `p`).

Let's say we want to define a manifold point `p_i` some rotation θ from the [`identity`](@ref) reference rotation `xR0` (another point on the manifold that we will use as reference)
```julia
# + radians rotation from x-axis on plane to point i
xθi = π/6
```

!!! note
The angle θ is easier to envision with `Rotations(2)`. Depending on the manifold, more generalized notions of distance between points exist.

### From Coordinates

To get our first Lie algebra element using the text book [`hat`](@ref), or equivaliently a more generalized [`get_vector`](@ref), function:
```julia
X_ = hat(G, xR0, xθi) # specific definition to Lie groups
xXi = get_vector(G, xR0, xθi, B) # generalized definition beyond Lie groups
# 2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
# 0.0 -0.523599
# 0.523599 0.0
@assert isapprox( X_, xXi )
```

Note that `hat` here assumes a default basis for the more general `get_vector`.

!!! note
In this case, the same would work given the base manifold [`Rotations(2)`](@ref):
```julia
_X_ = hat(M, xR0, xθi) # Lie groups definition
_X = get_vector(M, xR0, xθi, B) # generalized definition
@assert _X_ == xXi; @assert _X == xXi
```
One more caveat here is that for the Rotation matrices, the tangent vectors are always stored as elements from the Lie algebra.

Now, let's place this algebra element on the manifold using the exponential map [`exp`](@ref):
```julia
xRi = exp(G, xR0, xXi)
# similarly for known underlying manifold
xRi_ = exp(M, xR0, xXi)

@assert isapprox( xRi, xRi_ )
```

### To Coordinates

The logarithm transform from the group back to algebra (or coordinates) is:
```julia
xXi_ = log(G, xR0, xRi)
xXi__ = log(M, xR0, xRi)
@assert xXi == xXi__
```

Similarly, the coordinate value can be extracted from the algebra using [`vee`](@ref), or directly from the group using the more generalized [`get_coordinates`](@ref):
```julia
# extracting coordinates using vee
xθi__ = vee(G, xR0, xXi_)[1]
_xθi__ = vee(M, xR0, xXi_)[1]

# OR, the preferred generalized get_coordinate function
xθi_ = get_coordinates(G, xR0, xXi_, B)[1]
_xθi_ = get_coordinates(M, xR0, xXi_, B)[1]

# confirm all versions are correct
@assert isapprox( xθi, xθi_ ); @assert isapprox( xθi, _xθi_ )
@assert isapprox( xθi, xθi__ ); @assert isapprox( xθi, _xθi__ )
```

!!! note
The disadvantage might be that the representation of `X` is not nice, i.e. it uses too much space or doing vector-calculations is not so easy. E.g. fixed rank matrices are overloaded for all vector operations, but maybe that is “not enough” for a general user application that really wants vectors. But: Given a basis `B` one can look at the coefficients of the tangent vector `X` with respect to basis `B`. From the Sphere example note above the basis would be `Y1=[0.0,1.0,0.0]` and `Y2=[0.0,0.0,1.0]`, the so to get the coordinates would be `c = get_coordinates(Sphere(2), p, X, B)`. Visa versa, if you have a coordinate vector with respect to a basis `B` of the tangent space at `p` and want the vector back, then you do `X2 = get_vector(M, p, c, B)` (and you should have `X2==X`). The coordinate vector `c` might also have the advantage of saving memory. E.g. SPD matrix tangent vectors take n^2 entries to save, i.e. storing the full matrix, but the coordinate vectors only take n(n+1)/2.
dehann marked this conversation as resolved.
Show resolved Hide resolved

### Actions and Operations

With the basics in hand on how to move between the coordinate, algebra, and group representations, let's briefly look at composition and application of points on the manifold. For example, a `Rotation(n)` manifold is the mathematical representation, but the points have an application purpose in retaining information regarding a specific rotation.
dehann marked this conversation as resolved.
Show resolved Hide resolved

Points from a Lie group may have an associated action (i.e. a rotation) which we [`apply`](@ref). Consider rotating through `θ = π/6` three vectors `V` from their native domain `Euclidean(2)`, from the reference point `a` to a new point `b`. Engineering disciplines sometimes refer to the action of a manifold point `a` or `b` as reference frames. More generally, by taking the tangent space at point `p`, we are defining a local coordinate frame with basis `B`, and should not be confused with "reference frame" `a` or `b`.
Copy link
Member

Choose a reason for hiding this comment

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

In the language of Manifolds.jl, what the code below does is rotating three points aV1, aV2 and aV3 from Euclidean(2) around the origin (point [0, 0]). These points can, of course, be identified with vectors (specifically, free vectors that have two coordinates. You can also keep separately their origins if they are meant to represent bound vectors. FYI, I'm not particularly attached to this bound/free vector terminology so feel free to use different names, I just added these links to make the point more clear.

Manifolds.jl also offers a more differential-geometric way for representing bound vectors. Vectors aV1, aV2 and aV3 that are bound at a can be thought of as belonging to the tangent space at a, that is TangentSpaceAtPoint(Euclidean(2), a). If you then wanted to rotate aV1, aV2 and aV3 through a rotation that sends a to b, you'd actually apply differential of the group action on them, more or less like that:

using Manifolds, StaticArrays

a = [2, 3]
aV1 = [1;0]
aV2 = [0;1]
aV3 = [10;10]
G = SpecialOrthogonal(2)
A_left = RotationAction(Euclidean(2), G)
bθa = π/6

p0 = @SMatrix [1.0 0; 0 1]
xR0 = identity(G, p0)

B = DefaultOrthogonalBasis()
bXa = get_vector(base_manifold(G), xR0, bθa, B)
bRa = exp(G, xR0, bXa)

b = apply(A_left, bRa, a)
bV1 = apply_diff(A_left, bRa, a, aV1)

Of course, this is an extremely convoluted way of rotating a vector but that's what actually can be extended to arbitrary manifolds and actions. I don't know how useful would that perspective be, so don't take it as a recommendation for this tutorial, just a point that might be worth keeping in mind.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

offers a more differential-geometric way for representing bound vectors

I agree with that, let me try work this in. Will take a bit before I can though, busy week coming up this side.

I don't know how useful would that perspective be, so don't take it as a recommendation for this tutorial, just a point that might be worth keeping in mind.

I think it is important to get the language right, otherwise this will always be something that comes up with more people wanting to add contributions. Better to get it right now :-) thanks!


Keeping with our two-dimensional example above:
```julia
aV1 = [1;0]
aV2 = [0;1]
aV3 = [10;10]

A_left = RotationAction(Euclidean(2), G)

bθa = π/6
bXa = get_vector(base_manifold(G), xR0, bθa, B)

bRa = exp(G, R0, bXa)
Copy link
Member

Choose a reason for hiding this comment

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

Where is R0 defined? Did you mean xR0?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Oops, yes I think so. I also want to review all these variables again, especially in the context of the next docs in PR #366 which needs a lot of variable name updates.

What I'd like to get right here is that the language from rigid transforms (aka coordinate transformations) through to the differential geometry (and also differential forms) is as consistent as possible. The most important thing to get right now, in my view, is how folks understand all the conventions with a clear path to the future. I think manifold and differential form/geometry language is the right way to go. Getting consistency across terminology and convention is hard, but should eventually be the stronger outcome.

Copy link
Member

Choose a reason for hiding this comment

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

Yes, naming is hard, especially when we need to accommodate for people with different backgrounds. We have a long (and growing) table with notation used in Manifolds.jl: https://juliamanifolds.github.io/Manifolds.jl/stable/misc/notation.html but it still assumes basic knowledge of differential geometry. An introduction aimed at people who work with rigid transforms would definitely be valuable.

On the topic of differential forms, so far I haven't found any use for them outside of integration on manifolds so they are not available in Manifolds.jl yet but it's definitely something I'd like to see here 🙂 .

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sounds good, my hope is that these two tutorial examples "how to rotations" and the next "how to rigid transforms" give enough of runway to learn the terminology and pick up the JuliaManifolds packages.

Aside, over at IncrementalInference.jl we are working to adopt the differential forms / manifolds language as default (at least internally), while most of our audience is familiar with vector calculus and linear algebra lingo. All the issue comments and documentation you guys have added the past couple of weeks have been a great help, thanks!

Copy link
Member

Choose a reason for hiding this comment

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

And the other way around: Thanks for the questions and hence improving the documentation.


for aV in [aV1; aV2; aV3]
bV = apply(A_left, bRa, aV)
# test we are getting the rotated vectors in Euclidean(2) as expected
@assert isapprox( bV[1], norm(aV)*cos(bθa) )
@assert isapprox( bV[2], norm(aV)*sin(bθa) )
end
```

!!! note
In general, actions are usually non-commutative and the user must therefore be weary of [`LeftAction`](@ref) or [`RightAction`](@ref) needs. As in this case, the default `LeftAction()` is used.

Finally, the actions (i.e. points from a manifold) can be [`compose`](@ref)d together. Consider putting together two rotations `aRb` and `bRc` such that a single composite rotation `aRc` is found. The next bit of code composes five rotations of `π/4` increments:
```julia
A_left = RotationAction(M, G)
aRi = deepcopy(xR0)

iθi_ = π/4
x_θ = get_vector(M, xR0, iθi_, B) #hat(Rn, R0, θ)
iRi_ = exp(M, xR0, x_θ)

# do 5 times over:
# Ri_ = Ri*iRi_
for i in 1:5
aRi = compose(A_left, aRi, iRi_)
end

# drop back to a algebra, then coordinate
aXi = log(G, xR0, aRi)
aθi = get_coordinates(G, xR0, aXi, B)

# should wrap around to 3rd quadrant of xy-plane
@assert isapprox( -3π/4, aθi[1])
```

!!! warning
`compose` or `apply` must be done with group (not algebra) elements. This example shows how these two element types can easily be confused, since both the manifold group and algebra elements can have exactly the same data storage type -- i.e. a 2x2 matrix.

As a last note, other rotation representations, including quaternions, Pauli matrices, etc., have similar features. A contrasting example in rotations, however, are Euler angles which can also store rotation information but quickly becomes problematic with familiar problems such as ["gimbal-lock"](https://en.wikipedia.org/wiki/Gimbal_lock).