Usage and Testing of rotation.jl

By: Yi-Xin Liu

Re-posted from: http://www.yxliu.group/2020/04/rotation-test

Numerical experiments are carried out to demonstrate the usage as well as serving as a test for the rotation.jl submodule. This is an appendix to the previous post. The format of this post mimics the Jupyter notebook. The output of a code block after run is captured by Literate.jl and rendered here as a markdown code block in the format of the text language.

First, we load some necessary packages

1
2
3
4





using Revise
using Scattering
using LinearAlgebra
using Test

Some Experiments on Rotation Operations

$\vu, \vv, \vw$ are the basis vectors of the UVW internal coordinate system with their components expressed in the reference XYZ coordinate system. Let $\vw$ point to the principle direction of a scatterer.

1
2
3
4
5
6
7
8
9
10
11





w = [1.0, 2.0, 3.0]
normalize!(w)
# pick a vector that is not parallel to w
a = [0, 1, 2]
# find u by cross product
u = cross(w, a)
normalize!(u)
# find v perpendicular to both w and u
v = cross(w, u)
normalize!(v)
w









1
2
3
4





3-element Array{Float64,1}:
 0.2672612419124244
 0.5345224838248488
 0.8017837257372732

$\mP$ is the rotation matrix which rotates XYZ coordinate system (basis vectors) to UVW coordinate system: $[\vu\;\vv\;\vw] = [\ve_x\;\ve_y\;\ve_z]\mP$.

1





P = hcat(u, v, w)









1
2
3
4





3×3 Array{Float64,2}:
  0.408248   0.872872  0.267261
 -0.816497   0.218218  0.534522
  0.408248  -0.436436  0.801784

$\mR$ is the rotation matrix which rotates a vector in XYZ system to UVW coordinate system: $\mR = \mP^{-1} = \mP^T$ and $\mR\vu = [1\;0\;0]^T, \mR\vv = [0\;1\;0]^T, \mR\vw = [0\;0\;1]^T$

1
2





@test transpose(P) ≈ inv(P)
R = transpose(P)









1
2
3
4





3×3 LinearAlgebra.Transpose{Float64,Array{Float64,2}}:
 0.408248  -0.816497   0.408248
 0.872872   0.218218  -0.436436
 0.267261   0.534522   0.801784









1





@test det(R) ≈ 1









1





Test Passed

Verify that $\mR$ indeed transforms basis vectors of the XYZ coordinate system to $\vu, \vv, \vw$, whose components are expressed in the XYZ coordinate system.

1
2
3





@test R * u ≈ [1, 0, 0]
@test R * v ≈ [0, 1, 0]
@test R * w ≈ [0, 0, 1]









1





Test Passed









1





tr(R)









1





1.4282499064371286

Compute the rotation angle:

1
2





θ = acos((tr(R)-1)/2)
rad2deg(θ)









1





77.63580469304388

Compute the rotation axis expressed in the XYZ system:

1





uu = [R[3,2]-R[2,3], R[1,3]-R[3,1], R[2,1]-R[1,2]] / (2sin(θ))









1
2
3
4





3-element Array{Float64,1}:
 0.49700656431759865
 0.07216735383016143
 0.8647406247345901

Alternative way to compute the rotation axis. Note that the result may be different from the one computed above with only a sign.

1
2
3
4





vals, vecs = eigen(R)
idx = findfirst(x -> x ≈ one(x), vals)
uu2 = real(vecs[:, idx])
@test uu2 ≈ uu || uu2 ≈ -uu









1





Test Passed

Alternative way to compute the rotation angle

1
2





vv = [R[3,2]-R[2,3], R[1,3]-R[3,1], R[2,1]-R[1,2]]
@test asin(norm(vv)/2) ≈ θ









1





Test Passed

A proper rotation matrix should have $\det(\mR) = +1$

1





det(R)









1





1.0000000000000002









1





transpose(R)









1
2
3
4





3×3 Array{Float64,2}:
  0.408248   0.872872  0.267261
 -0.816497   0.218218  0.534522
  0.408248  -0.436436  0.801784









1





inv(R)









1
2
3
4





3×3 LinearAlgebra.Transpose{Float64,Array{Float64,2}}:
  0.408248   0.872872  0.267261
 -0.816497   0.218218  0.534522
  0.408248  -0.436436  0.801784

A proper rotation matrix should be an orthogonal matrix

1





@test transpose(R) ≈ inv(R)









1





Test Passed

Convert current rotation matrix representation to Euler angles representation. Convention used: Z1Y2Z3

Procedure:

  1. Rotate about +z by eta (counter-clockwise in x-y plane),
  2. Rotate about the former y-axis (which is y’), counter clockwise in x’-z plane. Then,
  3. Rotate about the former z-axis (which is z’), counter-clockwise in x’-y’ plane.

See wiki page:

  1. https://en.wikipedia.org/wiki/Rotation_matrix
  2. https://en.wikipedia.org/wiki/Euler_angles#Intrinsic_rotations
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22





α = atan(R[2,3], R[1,3])
β = acos(R[3,3])
γ = atan(R[3,2], -R[3,1])
rad2deg.([α, β, γ])
c1 = cos(α)
c2 = cos(β)
c3 = cos(γ)
s1 = sin(α)
s2 = sin(β)
s3 = sin(γ)

Ra = zero(R)
Ra[1, 1] = c1*c2*c3 - s1*s3
Ra[1, 2] = -c3*s1 - c1*c2*s3
Ra[1, 3] = c1 * s2
Ra[2, 1] = c1*s3 + c2*c3*s1
Ra[2, 2] = c1*c3 - c2*s1*s3
Ra[2, 3] = s1*s2
Ra[3, 1] = -c3*s2
Ra[3, 2] = s2*s3
Ra[3, 3] = c2
@test Ra ≈ R









1





Test Passed









1





rad2deg.([α, β, γ])









1
2
3
4





3-element Array{Float64,1}:
 -46.91127686463718
  36.69922520048988
 116.56505117707799

Convention: Z1X2Z3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22





α = atan(R[1,3], -R[2,3])
β = acos(R[3,3])
γ = atan(R[3,1], R[3,2])
rad2deg.([α, β, γ])
c1 = cos(α)
c2 = cos(β)
c3 = cos(γ)
s1 = sin(α)
s2 = sin(β)
s3 = sin(γ)

Rb = zero(R)
Rb[1, 1] = c1*c3 - c2*s1*s3
Rb[1, 2] = -c1*s3 - c2*c3*s1
Rb[1, 3] = s1 * s2
Rb[2, 1] = c3*s1 + c1*c2*s3
Rb[2, 2] = c1*c2*c3 - s1*s3
Rb[2, 3] = -c1*s2
Rb[3, 1] = s2*s3
Rb[3, 2] = c3*s2
Rb[3, 3] = c2
@test Rb ≈ R









1





Test Passed

Convention Z1Y2X3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22





α = atan(R[2,1],R[1,1])
β = asin(-R[3,1])
γ = atan(R[3,2], R[3,3])
rad2deg.([α, β, γ])
c1 = cos(α)
c2 = cos(β)
c3 = cos(γ)
s1 = sin(α)
s2 = sin(β)
s3 = sin(γ)

Z1Y2X3 = zero(R)
Z1Y2X3[1, 1] = c1*c2
Z1Y2X3[1, 2] = c1*s2*s3 - c3*s1
Z1Y2X3[1, 3] = s1*s3 + c1*c3*s2
Z1Y2X3[2, 1] = c2*s1
Z1Y2X3[2, 2] = c1*c3 + s1*s2*s3
Z1Y2X3[2, 3] = c3*s1*s2 - c1*s3
Z1Y2X3[3, 1] = -s2
Z1Y2X3[3, 2] = c2*s3
Z1Y2X3[3, 3] = c2*c3
@test Z1Y2X3 ≈ R









1





Test Passed

Convention: Z1Y2Z3

1
2
3
4
5
6
7
8
9





Rp = hcat([-1, 0, 0], [0, 0, 1], [0, 1, 0])
α = atan(Rp[2,3], Rp[1,3])
β = acos(Rp[3,3])
γ = atan(Rp[3,2], -Rp[3,1])
rad2deg.([α, β, γ])

θp = acos((tr(Rp)-1)/2)
up =[Rp[3,2]-Rp[2,3], Rp[1,3]-Rp[3,1], Rp[2,1]-Rp[1,2]]
θpp = asin(norm(up)/2)









1





0.0









1





rad2deg(θp), rad2deg(θpp), up









1





(180.0, 0.0, [0, 0, 0])

Find rotation axis by diagonalization

1
2
3





vals, vecs = eigen(Rp)
idx = findfirst(x -> x ≈ one(x), vals)
uu = real(vecs[:, idx])









1
2
3
4





3-element Array{Float64,1}:
 0.0
 0.7071067811865475
 0.7071067811865477

Testing and Usage of Scattering/rotation.jl

Testing EulerAngle constructors

1
2





using StaticArrays: SVector, FieldVector
α, β, γ = π/4, π/3, π/2









1





(0.7853981633974483, 1.0471975511965976, 1.5707963267948966)

Constructors of EulerAngle. Initialize by a vector

1





EulerAngle([α, β, γ])









1





Scattering.Rotation.EulerAngle{Float64}(0.7853981633974483, 1.0471975511965976, 1.5707963267948966)









1





EulerAngle(1.3, 0, 0)









1





Scattering.Rotation.EulerAngle{Float64}(1.3, 0.0, 0.0)

Initialize by three angles

1





euler = EulerAngle(α, β, γ)









1





Scattering.Rotation.EulerAngle{Float64}(0.7853981633974483, 1.0471975511965976, 1.5707963267948966)

Copy constructor

1





EulerAngle(euler)









1





Scattering.Rotation.EulerAngle{Float64}(0.7853981633974483, 1.0471975511965976, 1.5707963267948966)

Initialize by a static vector

1
2





sv = SVector(1.0, 2.0, 3.0)
EulerAngle(sv)









1





Scattering.Rotation.EulerAngle{Float64}(1.0, 2.0, 3.0)

Follow the Z1Y2Z3 convention

1
2
3
4





α = atan(v[3], u[3])
β = acos(w[3])
γ = atan(w[2], -w[1])
rad2deg.([α, β, γ])









1
2
3
4





3-element Array{Float64,1}:
 -46.91127686463718
  36.69922520048988
 116.56505117707799

Testing RotationMatrix constructors

Constructor of RotationMatrix initialized with a matrix

1
2





rotmat = RotationMatrix(Vector3D(u), Vector3D(v), Vector3D(w))
@test transpose(RotMatrix(hcat(u,v,w))) ≈ rotmat.R









1





Test Passed

Convert an EulerAngle instance to a RotationMatrix instance

1
2





rotmat_euler = RotationMatrix(EulerAngle(α, β, γ))
@test rotmat_euler.R ≈ rotmat.R









1





Test Passed

Convert a RotationMatrix instance to an EulerAngle instance

1
2





euler = EulerAngle(rotmat_euler)
@test [euler.α, euler.β, euler.γ] ≈ [α, β, γ]









1





Test Passed

Testing EulerAxisAngle Constructors

Constructor of AxisAngleRepresentation initialized with a RotationMatrix instance: conversion from rotation matrix representation to axis-angle representation

1





axisangle = EulerAxisAngle(rotmat)









1





Scattering.Rotation.EulerAxisAngle{Float64}([0.49700656431759865, 0.07216735383016143, 0.8647406247345901], 1.3550004093288814, [0.0 -0.8647406247345901 0.07216735383016143; 0.8647406247345901 0.0 -0.49700656431759865; -0.07216735383016143 0.49700656431759865 0.0])

Convert Euler axis-angle representation to rotation matrix representation

1
2





rotmat_axisangle = RotationMatrix(axisangle)
@test rotmat_axisangle.R ≈ rotmat.R atol=1e-15









1





Test Passed









1
2





axisangle2 = EulerAxisAngle(axisangle.ω, axisangle.θ)
@test axisangle2.K ≈ axisangle.K









1





Test Passed

Testing conversion and promotion

1





euler









1





Scattering.Rotation.EulerAngle{Float64}(-0.8187562376025611, 0.6405223126794245, 2.0344439357957027)









1





axisangle









1





Scattering.Rotation.EulerAxisAngle{Float64}([0.49700656431759865, 0.07216735383016143, 0.8647406247345901], 1.3550004093288814, [0.0 -0.8647406247345901 0.07216735383016143; 0.8647406247345901 0.0 -0.49700656431759865; -0.07216735383016143 0.49700656431759865 0.0])









1





EulerAxisAngle(euler)









1





Scattering.Rotation.EulerAxisAngle{Float64}([0.4970065643175987, 0.07216735383016142, 0.86474062473459], 1.3550004093288814, [0.0 -0.86474062473459 0.07216735383016142; 0.86474062473459 0.0 -0.4970065643175987; -0.07216735383016142 0.4970065643175987 0.0])









1





RotationMatrix(euler) |> EulerAxisAngle









1





Scattering.Rotation.EulerAxisAngle{Float64}([0.4970065643175987, 0.07216735383016142, 0.86474062473459], 1.3550004093288814, [0.0 -0.86474062473459 0.07216735383016142; 0.86474062473459 0.0 -0.4970065643175987; -0.07216735383016142 0.4970065643175987 0.0])









1





EulerAngle(axisangle)









1





Scattering.Rotation.EulerAngle{Float64}(-0.8187562376025609, 0.6405223126794247, 2.0344439357957027)









1





convert(EulerAngle, axisangle)









1





Scattering.Rotation.EulerAngle{Float64}(-0.8187562376025609, 0.6405223126794247, 2.0344439357957027)









1





convert(EulerAngle, euler)









1





Scattering.Rotation.EulerAngle{Float64}(-0.8187562376025611, 0.6405223126794245, 2.0344439357957027)

Testing promotion, comparison, and unit rotation

1





@test euler == axisangle









1





Test Passed









1





@test euler == euler









1





Test Passed









1





@test promote(euler) == rotmat









1





Test Passed









1





@test promote(axisangle) == rotmat









1





Test Passed









1





one(rotmat)









1





Scattering.Rotation.RotationMatrix{Float64}([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])









1





one(euler)









1





Scattering.Rotation.EulerAngle{Float64}(0.0, 0.0, 0.0)









1





promote(one(euler))









1





Scattering.Rotation.RotationMatrix{Float64}([1.0 -0.0 0.0; 0.0 1.0 0.0; -0.0 0.0 1.0])









1





one(axisangle)









1





Scattering.Rotation.EulerAxisAngle{Float64}([1.0, 0.0, 0.0], 0.0, [0.0 -0.0 0.0; 0.0 0.0 -1.0; -0.0 1.0 0.0])









1





promote(one(axisangle))









1





Scattering.Rotation.RotationMatrix{Float64}([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])

Testing inv function

1





rot1 = rotmat









1





Scattering.Rotation.RotationMatrix{Float64}([0.408248290463863 -0.816496580927726 0.408248290463863; 0.8728715609439696 0.21821789023599242 -0.4364357804719848; 0.2672612419124244 0.5345224838248488 0.8017837257372732])

Test inv function for RotationMatrix

1
2





irot1 = inv(rot1)
@test transpose(rot1.R) ≈ irot1.R









1





Test Passed









1





rot2 = EulerAxisAngle(rot1)









1





Scattering.Rotation.EulerAxisAngle{Float64}([0.49700656431759865, 0.07216735383016143, 0.8647406247345901], 1.3550004093288814, [0.0 -0.8647406247345901 0.07216735383016143; 0.8647406247345901 0.0 -0.49700656431759865; -0.07216735383016143 0.49700656431759865 0.0])

Test inv function for EulerAxisAngle

1
2





irot2 = inv(rot2)
@test irot2 == irot1









1





Test Passed









1





rot3 = EulerAngle(rot1)









1





Scattering.Rotation.EulerAngle{Float64}(-0.818756237602561, 0.6405223126794245, 2.0344439357957027)

Test inv function for EulerAngle

1
2





irot3 = inv(rot3)
@test irot3 == irot1









1





Test Passed

Testing multiplication of two AbstractRotation instances

Test multiplication of two RotationMatrix instances.

1
2
3
4





R1 = rotmat
M = RotationMatrix(Rp)
R2 = R1 * M
@test R2.R ≈ R1.R*M.R









1





Test Passed

Test multiplication of a EulerAxisAngle instance and a RotationMatrix instance.

1
2
3
4
5





R1 = rotmat
R1a = EulerAxisAngle(R1)
M = RotationMatrix(Rp)
R2 = R1a * M
@test R2 == R1 * M









1





Test Passed

Test multiplication of a RotationMatrix instance and a EulerAxisAngle instance.

1
2
3
4
5





R1 = rotmat
M = RotationMatrix(Rp)
Ma = EulerAxisAngle(M)
R2 = R1 * Ma
@test R2 == R1 * M









1





Test Passed

Test multiplication of a EulerAxisAngle instance and a EulerAxisAngle instance.

1
2
3
4
5
6





R1 = rotmat
R1a = EulerAxisAngle(R1)
M = RotationMatrix(Rp)
Ma = EulerAxisAngle(M)
R2 = R1a * Ma
@test R2 == R1 * M









1





Test Passed

Testing computing power of a AbstractRotation instance

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25





function pow1(rot::AbstractRotation, n::Integer)
    rot = promote(rot)
    result = rot
    for i in 1:(n-1)
        result *= rot
    end
    result
end

function pow(rot::RotationMatrix, n::Integer)
    if n == zero(n)
        return one(rot)
    elseif n == one(n)
        return rot
    else
        if isodd(n)
            return rot * pow(rot, n-1)
        else
            rothalf = pow(rot, n÷2)
            return rothalf * rothalf
        end
    end
end

pow(rot::AbstractRotation, n::Integer) = pow(promote(rot), n)









1





pow (generic function with 2 methods)









1





pow(euler, 0)









1





Scattering.Rotation.RotationMatrix{Float64}([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])









1





pow(euler, 1)









1





Scattering.Rotation.RotationMatrix{Float64}([0.40824829046386313 -0.8164965809277259 0.4082482904638629; 0.8728715609439694 0.21821789023599242 -0.4364357804719848; 0.26726124191242434 0.5345224838248487 0.8017837257372732])









1





pow(euler, 2)









1





Scattering.Rotation.RotationMatrix{Float64}([-0.4369210333151354 -0.2932896043722907 0.8503418245705543; 0.43018214432212887 -0.8983623348886722 -0.08881687880007882; 0.7899641342195538 0.32699590703939707 0.5186813505672173])









1





pow(euler, 3)









1





Scattering.Rotation.RotationMatrix{Float64}([-0.20711300761088602 0.7472703153125533 0.6314200487243415; -0.6322711178708509 -0.5947556020638506 0.4964866637886771; 0.7465503570320742 -0.29639981387703873 0.5956590591512387])









1





pow(euler, 4)









1





Scattering.Rotation.RotationMatrix{Float64}([0.7364715816744584 0.6696830270043788 0.09557328459445946; -0.644577211376976 0.651844177986726 0.3995239494677259; 0.2052555187063243 -0.35584239624735736 0.9117271308201462])









1





pow(euler, 5).R









1
2
3
4





3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
 0.910754   -0.404104   0.0850187
 0.412606    0.882094  -0.227304
 0.0168598   0.242097   0.970106









1





pow(euler, 100).R









1
2
3
4





3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
 -0.443094   0.414668   0.794807
 -0.277188  -0.906518   0.318422
  0.852546  -0.0792197  0.516614









1





euler^1









1





Scattering.Rotation.RotationMatrix{Float64}([0.40824829046386313 -0.8164965809277259 0.4082482904638629; 0.8728715609439694 0.21821789023599242 -0.4364357804719848; 0.26726124191242434 0.5345224838248487 0.8017837257372732])









1





euler^2









1





Scattering.Rotation.RotationMatrix{Float64}([-0.4369210333151354 -0.2932896043722907 0.8503418245705543; 0.43018214432212887 -0.8983623348886722 -0.08881687880007882; 0.7899641342195538 0.32699590703939707 0.5186813505672173])









1





euler^3









1





Scattering.Rotation.RotationMatrix{Float64}([-0.20711300761088602 0.7472703153125533 0.6314200487243415; -0.6322711178708509 -0.5947556020638506 0.4964866637886771; 0.7465503570320742 -0.29639981387703873 0.5956590591512387])









1





euler^4









1





Scattering.Rotation.RotationMatrix{Float64}([0.7364715816744584 0.6696830270043788 0.09557328459445946; -0.644577211376976 0.651844177986726 0.3995239494677259; 0.2052555187063243 -0.35584239624735736 0.9117271308201462])









1





@test pow(euler, 100).R ≈ (euler^100).R









1





Test Passed

Benchmarks

1
2
3





using BenchmarkTools

@btime pow1($euler, 1000)









1
2





142.561 μs (3004 allocations: 172.31 KiB)
Scattering.Rotation.RotationMatrix{Float64}([-0.17617351956438138 0.7712758102440354 0.6116342988556001; -0.6592241548072412 -0.553880454852226 0.50856656934094; 0.7310173764848875 -0.31360814126062625 0.606022713280731])

The implementation of pow (or Scattering.^ which is the same as pow here)
is significantly faster than the naive implementation.

1





@btime $euler^1000









1
2





2.327 μs (49 allocations: 3.02 KiB)
Scattering.Rotation.RotationMatrix{Float64}([-0.1761735195643819 0.7712758102440348 0.6116342988555798; -0.6592241548072477 -0.5538804548522176 0.5085665693409378; 0.7310173764848691 -0.3136081412606315 0.6060227132807026])

Testing applying AbstractRotation instance to a vector

Test multiplication of a RotationMatrix instance and a Vector

1
2
3





M = RotationMatrix(Rp)
v = [1.0, 2.0, 3.0]
@test M*v ≈ M.R * v









1





Test Passed

Test multiplication of a RotationMatrix instance and a RVector (SVector, QVector)

1
2
3





M = RotationMatrix(Rp)
r = RVector([1.0, 2.0, 3.0])
@test M*r ≈ M.R * r









1





Test Passed

Test multiplication of a EulerAxisAngle instance and a Vector

1
2
3
4





M = RotationMatrix(Rp)
Ma = EulerAxisAngle(M)
v = [1.0, 2.0, 3.0]
@test Ma*v ≈ M.R * v









1





Test Passed

Test multiplication of a EulerAxisAngle instance and a RVector (SVector, QVector)

1
2
3
4





M = RotationMatrix(Rp)
Ma = EulerAxisAngle(M)
v = RVector([1.0, 2.0, 3.0])
@test Ma*v ≈ M.R * v









1





Test Passed

Test multiplication of a EulerAngle instance and a Vector

1
2
3
4





M = RotationMatrix(Rp)
Ma = EulerAngle(M)
v = [1.0, 2.0, 3.0]
@test Ma*v ≈ M.R * v









1





Test Passed

Test multiplication of a EulerAngle instance and a RVector (SVector, QVector)

1
2
3
4





M = RotationMatrix(Rp)
Ma = EulerAngle(M)
v = RVector([1.0, 2.0, 3.0])
@test Ma*v ≈ M.R * v









1





Test Passed

Acknowledgements

  • This work is partially supported by the General Program of the National Natural Science Foundation of China (No. 21873021).
  • The captured output text blocks in this post is automatically generated by using Literate.jl.

Usage and Testing of rotation.jl was originally published by Yi-Xin Liu at Yi-Xin Liu on April 19, 2020.