# Representing rotations by quaternions

The most usual ways to represent a rotation in three-dimensional space are rotations matrices or Euler angles. In this package, nevertheless, they are represented using unit quaternions.

Quaternions are a four-dimensional non-commutative algebra with the property that any rotation may be represented by a quaternion of norm 1. Some of their advantages consist in that we only need 4 parameters to represent a quaternion (as opposed to the 9 needed by rotation matrices) and that they do not suffer from the phenomenon of gimbal lock (as opposed to Euler angles).

If a quaternion $q$ represents a rotation matrix $R$, their action on a vector is defined as

\[R v = q v q^{-1}.\]

This package assumes that your coordinate system follows the right-hand rule. If, for whatever reason, you want to use a left-handed system, $q$ and $q^{-1}$ must be transposed on the above formula.

## The `rotate`

function

All rotations on `DeformableBodies.jl`

are done by the function `rotate`

. Its signature consists of

`rotate(v::Vector, q::Quaternion, center::Vector)`

This function rotates the vector `v`

by the rotation represented by `q`

in the frame of reference centered on `center`

.

If the third argument is left blank, it defaults to the origin.

```
julia> v = [3,0,0]
3-element Array{Int64,1}:
3
0
0
julia> q = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> rotate(v, q, [0,0,0])
3-element Array{Float64,1}:
-2.0
2.0
1.0
julia> rotate(v, q)
3-element Array{Float64,1}:
-2.0
2.0
1.0
```

## Representing the identity rotation

The identity rotation is the matrix $I$ satisfying $I v = v$ for all $v$. This is represented by the quaternion $1$, which can be gotten using Julia's multiple dispatch via the expression `one(Quaternion)`

.

```
julia> q = one(Quaternion)
1 + 0i + 0j + 0k
julia> rotate([1,2,3], q)
3-element Array{Float64,1}:
1.0
2.0
3.0
```

## Composition of rotations

The composition of rotations translates to the quaternion world as the product of quaternions. Thus, it is the same to rotate a vector by $q_1$ and then by $q_2$ and to rotate by $q_2\,q_1$.

```
julia> v = [9, 0, 0]
3-element Array{Int64,1}:
9
0
0
julia> q1 = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> q2 = Quaternion(4,3,2,1)
4 + 3i + 2j + 1k
julia> rotate(v, q2*q1)
3-element Array{Float64,1}:
-1.0
-4.0
8.0
julia> rotate(rotate(v, q1), q2)
3-element Array{Float64,1}:
-1.0
-4.0
8.0
```

Rotations are non-commutative. Applying $q_2$ after $q_1$ is the same as multiplying $q_2$ to the **left** of $q_1$.

## Axis-angle representation

An important property of three-dimensional space it that every rotation fixes a line. Therefore, they accept an axis-angle representation. That is, a rotation $R$ is defined by a unit vector $\hat{n}$ and an angle $\theta$ such that $R$ rotates a vector counterclockwise by an amount of $\theta$ around the line defined by $\hat{n}$.

To get the quaternion representing an rotation of `θ`

around `n`

, use the method `axistoquaternion`

. You may pass an unormalized vector and the method will normalize it.

```
julia> q = axistoquaternion(axis = [0,0,2], angle = π/2)
0.7071067811865476 + 0.0i + 0.0j + 0.7071067811865475k
julia> rotate([1,0,0], q)
3-element Array{Float64,1}:
2.220446049250313e-16
1.0
0.0
```

In fact, this combination is so useful that the function `rotate`

is overloaded to directly convert from an axis-angle pair.

```
julia> rotate([1,0,0], axis=[0,0,1], angle=π/2)
3-element Array{Float64,1}:
2.220446049250313e-16
1.0
0.0
```

This version of `rotate`

also accepts the optional argument `center`

, which defaults to the origin.

`rotate([1,0,0], axis=[0,0,1], angle=π/2, center=[0,1,0])`

## Converting between matrices and quaternions

Sometimes a rotation may already come represented as a matrix or the transition rotations from `solve!`

may be needed in matrix form some application. To deal with these cases, the package provides two auxiliary functions `matrixtoquaternion`

and `quaterniontomatrix`

.

To convert a rotation matrix to a unit quaternion, do

```
julia> R = [cos(π/4) -sin(π/4) 0;
sin(π/4) cos(π/4) 0;
0 0 1]
3×3 Array{Float64,2}:
0.707107 -0.707107 0.0
0.707107 0.707107 0.0
0.0 0.0 1.0
julia> R * [sqrt(2), 0, 0]
3-element Array{Float64,1}:
1.0000000000000002
1.0
0.0
julia> q = matrixtoquaternion(R)
0.9238795325112867 + 0.0i + 0.0j + 0.3826834323650898k
julia> rotate([sqrt(2), 0, 0], q)
3-element Array{Float64,1}:
0.9999999999999999
1.0
0.0
```

Some minor differences may occur due to floating-point rounding errors.

The function `matrixtoquaternion`

assumes that the input is a *rotation matrix* but, for efficency reasons, no check is done in this regard. If you do not make sure beforehand that the matrix is orthogonal, bad things may happen.

To convert a quaternion to a matrix simply do

```
julia> q = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> R = quaterniontomatrix(q)
3×3 Array{Float64,2}:
-0.666667 0.133333 0.733333
0.666667 -0.333333 0.666667
0.333333 0.933333 0.133333
julia> rotate([3,0,0], q)
3-element Array{Float64,1}:
-2.0
2.0
1.0
julia> R * [3,0,0]
3-element Array{Float64,1}:
-2.0
2.0
1.0
```

The function `quaterniontomatrix`

works for every quaternion, and does not require the input to be a unit quaternion.

There is a **unique** rotation matrix representing a given quaternion but there are **two** unit quaternions representing the same matrix.

This means that `quaterniontomatrix ∘ matrixtoquaternion`

equals the identity (disconsidering floating-point rouding errors) but the opposite is not in general true. For a simple example.

```
q = Quaternion(-1.0)
(matrixtoquaternion ∘ quaterniontomatrix)(q)
```

Nevertheless, both these quaternions produce the same rotations.

## Rotating a PointMass

All the previous functionalities only require the submodule `Quaternions`

and work directly with vectors. Nevertheless, the models on `DeformableBodies.jl`

are constructed with respect to the type `PointMass`

. To help with that, the function rotate is overloaded to directly rotate the position of a `PointMass`

without interfering with its mass.

```
rotate(p::PointMass, q::Quaternion, center::Vector)
rotate(p::PointMass; axis::Vector, angle::Real, center::Vector)
```

The usage is identical to the version for vectors including the fact that the argument `center`

must be a Vector and not another PointMass.

An usual application consists in rotating a body around its center of mass.

```
julia> body = [ PointMass(rand(), rand(3)) for i in 1:5 ]
5-element Array{PointMass{Float64},1}:
PointMass{Float64}(0.027623475555394572, [0.622502144842952, 0.11921732534587282, 0.8843304426741481])
PointMass{Float64}(0.08235034045511802, [0.9932241501452683, 0.334779563612172, 0.8007546565223924])
PointMass{Float64}(0.16330720003612642, [0.386699947607279, 0.8099171674657109, 0.5908987643400341])
PointMass{Float64}(0.8618919657812625, [0.000659859384851913, 0.7755364055477827, 0.37741989797604925])
PointMass{Float64}(0.5806780634969708, [0.922564768890463, 0.04356399040074832, 0.18765714901619934])
julia> center_of_mass(body)
3-element Array{Float64,1}:
0.40704043526795597
0.4993749112573767
0.3619966261272671
julia> q = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> rotated = [ rotate(p, q, center_of_mass(body)) for p in body ]
5-element Array{PointMass{Float64},1}:
PointMass{Float64}(0.027623475555394572, [0.5957564162308041, 1.1179577906424625, 0.14864795800777955])
PointMass{Float64}(0.08235034045511802, [0.31606113462014573, 1.2375345239874034, 0.46226994400353016])
PointMass{Float64}(0.16330720003612642, [0.629867962558881, 0.5349019262226586, 0.6755761877965223])
PointMass{Float64}(0.8618919657812625, [0.7260927511178531, 0.14668287713769346, 0.48634359841711555])
PointMass{Float64}(0.5806780634969708, [-0.1252661931427156, 0.8787684558838791, 0.08516928092044052])
julia> center_of_mass(rotated)
3-element Array{Float64,1}:
0.40704043526795597
0.49937491125737676
0.361996626127267
```