Edit on GitHub

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}.\]
Warning

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
Note

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
Info

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.

Danger

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.

Note

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.5887392979926158, [0.9851461509508204, 0.5871472754686167, 0.18988173568255662])
 PointMass{Float64}(0.16447626657374115, [0.15917411711754115, 0.21781301182062407, 0.6154574686908942])
 PointMass{Float64}(0.6294684337802217, [0.6007476400909275, 0.8297934043579711, 0.011699945759655606])
 PointMass{Float64}(0.3674998635545208, [0.24634732086407585, 0.2279264662827194, 0.7397194326716405])
 PointMass{Float64}(0.14298635160866868, [0.6206875249356907, 0.26670397043925376, 0.6435613516386134])

julia> center_of_mass(body)
3-element Array{Float64,1}:
 0.6146350376396904
 0.54180467556897
 0.30860986988838707

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.5887392979926158, [0.2866060103346143, 0.6945457950059541, 0.4586029163376565])
 PointMass{Float64}(0.16447626657374115, [1.1000983352765157, 0.550726349005324, -0.104689643277118])
 PointMass{Float64}(0.6294684337802217, [0.4445578554829961, 0.23861021818764017, 0.5331822276913696])
 PointMass{Float64}(0.3674998635545208, [1.1344567669596526, 0.6883119760028129, -0.049624422666217916])
 PointMass{Float64}(0.14298635160866868, [0.8195510387085607, 0.8608408899763599, 0.0985269050993488])

julia> center_of_mass(rotated)
3-element Array{Float64,1}:
 0.6146350376396904
 0.5418046755689702
 0.308609869888387