Rotation Matrices 2
WCR 2017-04-21
Outline
Markers and Measured Orientation Vectors
Rotation-of-Points Matrix
Rotation-of-Coordinate-System Matrix
Summary
Helical Axis and Angles
Simulating Movement with Helical Angles
Euler Angles
Simulating Movement with Euler Angles
Shoulder Joint
Hip Joint
Markers and Measured Orientation Vectors
Marker locations are determined in the global (lab) reference system (GRS), using motion capture techniques. These marker locations are used to estimate orthogonal unit vectors iprox,jprox,kprox for the proximal segment and idistal,jdistal,kdistalfor the distal segment, at each time point. Since the markers are in the GRS, the unit vectors are also expressed in the GRS. We assume vectors are expressed as column vectors. We also assume that i,j,k of the proximal and distal segments are defined so that they are aligned in the neutral position. Another way of saying this is that the joint angles (flexion/extension, abduction/adduction, internal/external rotation) are assumed to be zero when i,j,k point in the same directions for both segments. (The ankle joint may be an exception to this general rule.)
Our goal is to use our knowledge of the unit vector orientations at each time to determine the angles of flex/ext, abd/add, and IR/ER, of the distal segment relative to the proximal, at each time.
At each time, we arrange the orthogonal unit vectors of the proximal and distal segments into two 3x3 matrices. Each column is a unit vector expressed in the GRS:
(1)
(2)
Rotation-of-Points Matrix
The rotation matrix that will move the vectors in the proximal segment into alignment with the distal segment isRp2d :
(3)
(4)
(5)
(6)
In going from eq. 5 to eq. 6, we used the fact that R-1=RT when R is a rotation matrix. Note that the rotation matrix Rp2d rotates the vectors, or (equivalently) the points on the segment. Rp2d isnot a rotation of the coordinate system. The rotation matrices in Winter, 4th ed., chapter 7, are rotations-of-coordinate-systems, which leave the points themselves unmoved. Therefore we calculate Rp2d according to the equation 6. Saying it again, slightly differently, because it is important to understand:Rp2d is the matrix which, if we multiply it times the unit vectors in the proximal segment (as they are at one instant), will give us a set of unit vectors that are identical to the distal segment unit vectors, as they are at that instant. At most joints (the ankle being an exception), we assume the distal and proximal segment unit vectors are aligned in the neutral position. In that case, Rp2d tells us how the distal segment is rotated relative to its neutral orientation.
Rotation-of-Coordinate-System Matrix
An alternate approach to quantifying the relative orientation of segments is to consider rotation of the coordinate reference system, instead of rotation of the distal segment relative to the proximal. Mathematically, there is no particular reason to prefer rotation-of-points or rotation-of-coordinate-systems: either approach will successfully account for rotations in three dimensions. The rotation-of-points approach is, perhaps, easier to visualize, since one rotates the segment itself. However, the rotation-of-coordinate-system approach has a practical advantage for biomechanics: when we do successive rotations of the coordinate system, the second and third rotations will be relative to the rotated(distal segment) axes, because the rotation of coordinates “takes the axes along”. Rotations-of-points do not “take the axes along”, so successive rotations are always about the original axes, and not about the moving distal segment axes.
Imagine rotating the coordinate system along with the distal segment. When we multiply the rotation-of-coordinate-system matrix by a vector, we get the expression of that vector in the rotated reference system. Therefore, the rotation-of-coordinate-system matrix representing the combined movement of the distal segment relative to the proximal will express the unit vectors of the proximal segment (represented by the identity matrix) in the distal segment coordinatesystem. Put another way, the columns of the rotation-of-coordinate-system matrix express what iprox, jprox, and kprox “look like” to an observer in the distal segment coordinate system.We know, from experimental observation of markers, Rprox and Rdistal, whose columns are the unit vectors of the proximal and distal segments, expressed in the GRS. The following equations show how we can use that information to express iprox, jprox, kprox in terms of idistal, jdistal, kdistal.
(7)
(8)
(9)
The above three equations can be written as a matrix equation:
(10)
where iprox, etc., are column vectors. Matrix Q≡(qij) is the rotation-of-coordinate-system matrix. The columns of Q are the proximal coordinate system unit vectors, expressed in terms of the distal coordinate system. We apply equations 1 and 2 to equation 10:
(11)
(12)
(13)
Compare equation 6 to equation 13. Equation 6 is the rotation-of-points matrix that moves the proximal segment unit vectors into alignment with the distal coordinate system vectors. Q in equation 13 is the rotation-of-coordinate-system matrix, that expresses the proximal segment unit vectors in terms of the distal segment coordinate system.
Summary
We have demonstrated two ways to represent the relative orientation of the distal and proximal segments: by the rotation-of-points (Rp2d) or by the rotation-of-coordinate-systems (Q). We need to estimate the proximal and distal segment unit vectors in the GRF (Rprox and Rdistal) to estimate Rp2d or Q. The next step is to use Rp2d or Q to estimate the joint angles,
Helical Axis and Angles
In the helical approach, we will do a single rotation about an axis which may not be aligned with any of the principle axes in either segment. It is always true that for any non-zero rotation of a segment, there is a single axis n=(nx;ny;nz) and an associated angle such that a rotation by angle about n will move the proximal i,j,k into the distal i,j,k. (If the distal segment is in the neutral position relative to the proximal, then the angle will be zero and the axis will be undefined.) If n is a unit vector (which we can assume without loss of generality), then the components of n along the X, Y, and Z axes tell us how to “allocate” to flexion/extension, abduction/adduction, and internal/external rotation. Now we show how to determine n=(nx;ny;nz) from Rp2d. (Remember that we have already determined Rp2d from marker locations.)
The rotation-of-points matrix arising corresponding to rotation by about n=(nx;ny;nz) is
(14)
which can also be written as
(15)
( retreived 2013-03-22) where
e0=cos(/2)(16a)
e1=nxsin(/2), e2=nysin(/2), e3=nzsin(/2)(16b)
The helical angle approach finds the one single axis n=(nx;ny;nz), about which we can rotate the segment to move it from its neutral to its observed position, and it finds the angle of rotation, . We wish to find the unit vector n and the rotation angle which will make the theoretical rotation-of-points matrix Rhelical equal to the measured rotation-of-points matrix Rp2d:
(17)
when
(18)
(19)
It can be shown by algebra ( retrieved 2013-03-24) that each column of the following matrix is a scaled version of the rotation axis, n.
(20)
where I is the identity matrix. Combining equations 17, 18, and 20, we get
(21)
(22)
To minimize round-off errors, estimate n by choosing the column from(22) above which has the largest magnitude. Normalize that column vector to have unit length:
(23)
where nmax = column of N with largest magnitude = maxnorm(n1, n2, n3).
Before proceeding further, check that the values for n and prelim generate a rotation matrix equal to Rp2d, by plugging n and prelim into equation 14 and computing Rhelix. If Rhelix does not match Rp2d to good precision, set =-prelim. This is necessary because the estimation of by equation 19 may result in a sign error for . You might think that n would then point the opposite way, in which case we would still get correct answers for angles, but this is not always the case.We do not want small round-off differences between matrices to make us erroneously conclude that the sign of ϕprelim needs changing. Therefore we use equations 14 to compute two rotation-of-points matrices: Rhelix(n,+ϕprelim) and Rhelixn,-ϕprelim). We choose the sign of ϕ that gives a matrix that is “closer to” the matrix Rp2d , the rotation of points matrix computed with equation 6. We define “closer to” by computing
(23a)
and
(23b)
where norm denotes the 2-norm, or equivalently, the Euclidean norm, of a matrix. If d- is smaller than d+, we choose ϕ=ϕprelim ; otherwise, we choose ϕ=+ϕprelim.
The projections of n onto the proximal segment axes are the respective amounts of flex/ext, abd/add, and IR/ER. If +X points right (flex/ext axis), +Y points forward (add/abd axis), and +Z points up (IR/ER axis), then
Flexion/extension:(24)
Adduction/abduction:(25)
Internal/external rotation:(26)
The algorithm just described, for finding the helical axis and angles, fails when the proximal and distal segments are exactly aligned. It is not surprising that it fails, because when the segments are aligned, the rotation angle should be zero, and in that case there is no uniquely defined axis of rotation. Mathematically, the problem occurs in equation 23: n=nmax/||nmax||. When the proximal and distal segments are aligned, Rp2d is the identity matrix, and when Rp2d=I, the matrix N in equations 20-22 is the zero matrix, and so the vector nmax has length 0, and equation 23 has an undefined (divide by zero) result. This problem is handled by checking for Rp2d=I. If Rp2d=I, the helical angles are set equal to zero, and the axis can be arbitrarily set to n=(1;0;0), or can be estimated by interpolating from nearby time points at which the axis is well defined.
Simulating Movement With Helical Angles
Simulation of movement using helical angles is fairly straightforward. Create a 3xN matrix P, whose columns Pj, j=1..N, are the N points in an object. Choose the axis of rotation n and the desired rotation angle . Use equations 15, 16, and 14 to create the rotation-of-points matrix R. Rotation about the origin moves the points P toProtated:
(27)
where P and Protated are both expressed in the same coordinate system, such as the global reference system. Rotation about a point C is achieved as follows, assuming P, C, and Protated are all expressed in the same coordinate system:
(28)
where (C…C) denotes a 3xN matrix composed of N replicates of the column vector C.
Euler Angles
One can decompose a rotation in three dimensions into three successive rotations about axes that are anatomically defined. This is in contrast to the helical approach, in which a single rotation is done about an axis that is not anatomically defined. When using Euler angles, we must define axes and choose the order of rotations. The order of rotation matters. The order dependence can be understood mathematically as a result of the fact the matrix multiplication is not commutative: if A and B are matrices, it is generally true that AB≠BA. The order dependence can be understood physically as meaning that, for example,45° flexion followed by 45° abduction (about the rotated ab/adduction axis) will result in a different segment orientation than 45° abduction followed by 45° flexion (about the rotated flexion/extension axis). The angles of rotation are called Euler angles.
We choose to use matrices that represent rotating the coordinate system while leaving the points in place: the rotation-of-coordinate-system matrices. We choose this even though it may be easier to visualize the action of rotation-of-points, which moves the segment itself. To recapitulate the reason, which was also described above: if we do successive rotations of points, each successive rotation is about the original coordinate axes, since points are always expressed in the original coordinate system. This means that the combined rotation matrix RpointsXYZ=RpointsZRpointsYRpointsX, for example, represents rotation about proximal X, then proximal Y, thenproximal Z. That’s not what we want. We want to rotate about X (which is the same for proximal and distal segments initially), then about Y’= distal segment Y, then about Z”=twice-rotated distal segment Z. For this reason, we must use rotation-of-coordinate-system matrices.
There are two kinds of rotation sequences that can be used to account for any rotation in three dimensions: Euler sequences and Cardan sequences. Euler sequences use a shared axis in the aligned proximal and distal segments for the first rotation, then an axis in the distal segment, which has been rotated by the first rotation, for the second rotation, and finally, the same distal segment axis that was used for the first rotation is used again for the final rotation – but due to the first two rotations, this axis now points in a different direction than it did initially. An example of such a sequence is X-Y’-X’, or simply X-Y-X. This sequence indicates a rotation about X, followed by a rotation about the rotated Y axis (called Y’ to indicate that it moved during the first rotation), and then about the rotated X axis (denoted X’). There are six possible Euler sequences: XYX, XZX, YXY, YZY, ZXZ, ZYZ. Euler sequence rotations can not be accomplished with rotation-of-points. If one attempts to do so, one finds that the first and final rotation axes are the same, and thus there are really only two independent axes of rotation. A rotation sequence with two pre-defined axes is not sufficiently general to account for most rotations in three dimensions.
The second type of rotation sequence is called a Cardan sequence. A Cardan sequence involves rotation about three different axes. The first one is shared between proximal and distal segment, since they are presumed to be aligned initially. The second and third rotations are about the other two orthogonal axes, which are fixed in the distal segment. Because they are fixed in the distal segment, they move with each rotation. An example of a Cardan sequence is X-Y’-Z”, or simply XYZ, denoting an initial rotation about the X axis, then about the rotated Y-axis in the distal segment (denoted Y’), and finally about the twice-rotated Z-axis (Z”) of the distal segment. There are six possible Cardan sequences: XYZ, XZY, YXZ, YZX, ZXY, ZYX.
With either kind of sequence, the rotation can be represented by multiplying three rotation matrices.
The matrices that represent rotations of the coordinate system about an axis (rather than rotations of the points themselves) are as follows:
(29)
(30)
(31)
The rotation matrices above operate on points represented as column vectors, and therefore the rotation matrix is to the left of the column vector, or the matrix of points. When several rotations are done, the matrix for the first rotation is right-most and the last rotation is left-most. This allows the first (right-most) matrix to operate first on the point or points. For example, expression of point in a coordinate system that is rotated about X, then Y’, then Z’’ is represented with matrices as follows:
(32)
(33)
where Rtotal is the rotation-of-coordinate-system matrix, and P0 and Pfinal are the same vector, expressed in the initial (proximal) and final (distal) coordinate systems respectively.Rtotal depends on the rotation sequence chosen, as shown in the following table.
Cardan Sequences / Rtotal =X then Y’ then Z” /
Y then X’ then Z” /
X then Z’ then Y” /
Z then X’ then Y” /
Y then Z’ then X” /
Z then Y’ then X” /
Euler Sequences
X then Y’ then X” /
Y then X’ then Y” /
X then Z’ then X” /
Z then X’ then Z” /
Y then Z’ then Y” /
Z then Y’ then Z” /
Table 1. Rotation-of-coordinate-system matrices for all Cardan and Euler sequences.
If P0 and Pfinal are direction vectors, as is the case when we are estimating joint angles from segment unit vectors, then it does not matter if the two systems share an origin or not, since direction vectors are independent of origin. If P0 and Pfinal are position vectors, then the origin matters, and equation 33will be correct if and only if the proximal and distal systems share an origin.
We will now complete the explanation of how to use Rtotaland the matrix Q, which we obtain from motion capture data, to estimate Euler angles. We use the rotation-of-coordinate-system matrices (equations 29-31) for RX, RY, and RZ to compute the theoretical Rtotal for our preferred rotation sequence. Then we assign values to the angles that will make the theoretical matrix Rtotal equal to the experimentally measured matrix Q.
Euler example 1
If we choose to rotate about X, then Y’, then Z”, then we see from Table 1 that
(34)
which after plugging in from equations 29-31 becomes
(35)
(36)
(37)
Choose the angles X,Y,Z to make RXYZ equal the measured Q, from equation 13:
(38)
RXYZ equals Q when
(39)
(40)
(41)
where atan2(y,x) = four quadrant arctangent function. With the above definitions, X and Z will be in the range [-,+], and Ywill be in the range [-/2,+/2].
Euler example 2
If we reverse the order of the first and second rotations, we will rotate about Y, then X’, then Z”. Using the same approach as before, we now obtain
(42)