University of Amsterdam - Amsterdam Medisch
Centrum - Biomedical Engineering and Physics
Bachelor Project Report Physics and Astronomy
Artifact in 3D-based position analysis using Euler angles
Conducted between: 09-9-2019 and 12-09-2019 Credits: 15 ec
Author: Isali Kant - Studentid. 10716505
Supervisor: G.J. Streekstra Daily Supervisor: J.G.G. Dobbe
Examiner: G.J. Strijkers
Abstract
Fracture of the distal radius is a common injury as a result of falling on outstretched arms. When the bone is healed in a malunited manner and results in serious impairment, corrective surgery is to be considered. Through preoperative planning and surgery with a patient specific casts, the precision of the surgery is increased and the invasiveness reduced. To test the precision of the method, CT scans were made and the data was evaluated. Coor-dinate systems are assigned to the marker tools on the patient-specific cast and to the distal radius. The rotation from one CS to another is expressed in Euler angles. There are various sequences of rotation possible to analyse the precision of the method. The Euler angles from the rotation sequence used in the research by G. Caiti have been reproduced, and the sets of Euler angles for the all six rotation sequences have been calculated as well. The discrep-ancy found in the research has been checked for Gimbal Lock. In that research, the rotation sequence YXZ was used with a total standard deviation of 3,503, while the best rotation se-quence would have been sese-quence ZYX with total standard deviation of 3,441. It is therefore important that all rotation sequences are analysed when using an Euler angle representation. Suggestions for other representations such as axis angle and quaternions are presented.
Samenvatting voor scholieren
In het plaatje hieronder is de anatomie te zien van de onderarm en pols. Een breuk van de radius, het rode bot in het plaatje, komt veel voor als gevolg van vallen op uitgestrekte armen. Dit wordt behandeld met gips, of het bot wordt gefixeerd met pinnen. Het kan zijn dat het bot toch niet op de goede manier vastgroeit, waardoor je last van je pols kan krijgen. In dit geval kunnen artsen er voor kiezen je te opereren, en het bovenste deel van het bot in een an-dere positie te brengen. Voor de operatie worden er foto’s gemaakt van je onderarm. Via de computer wordt er een planning gemaakt hoe het bot in positie gebracht moet worden. Met behulp van een mal om de arm, wordt deze planning overgebracht naar je arm. Onderzoekers zijn bezig technieken te ontwikkelen waarbij de operatie zo min van je arm open legt. Voor een nieuw ontwikkelde mal, zijn tests gedaan hoe precies de planning overgebracht kan wor-den op een arm. Hierbij wordt er berekend welke rotatie het bovenste deel van het bot moet maken, ten opzichte van de mal. Deze rotatie wordt uitgedrukt in drie hoeken, die een ro-tatie voorstellen rond de x, y en z as. In dit onderzoek is uitgerekend wat het verschil is in verschillende rotatievolgordes van die drie hoeken.
Contents
1 Introduction 4
1.1 The clinical case . . . 4
1.1.1 Purpose . . . 5 1.1.2 Methods . . . 5 1.1.3 Results . . . 7 1.2 Motivation . . . 9 2 Theoretical framework 10 2.1 Rotation matrices . . . 10
2.2 Computing Euler angles from a rotation matrix . . . 11
2.3 Challenges while using the Euler angle representation . . . 13
2.3.1 Tait-Bryan rotations . . . 13
2.3.2 Intrinsic and extrinsic rotations . . . 13
2.3.3 Rotation axes . . . 13 2.3.4 Gimbal Lock . . . 14 2.3.5 Rotation Order . . . 15 3 Experiments 16 4 Discussion 23 5 Conclusion 24 6 References 25 7 Acknowledgements 26 8 Appendices 27 8.1 Appendix A . . . 27
1
Introduction
In the field of biomedical engineering and physics, the orientation of a 3D object with respect to a reference coordinate system (CS) is usually represented by so called Euler angles. The precision of the method is often expressed in terms of the standard deviation of these angles. Unfortunately, because of intrinsic properties of Euler angle calculations, there are some diffi-culties in using this representation.
In this thesis we will study a case where the rotation matrices yielded Euler angles that did not seem to be reliable. The precision of the method used in the research by G. Caiti [2] is expressed in terms of the standard deviation on the calculated Euler angles. Paragraph (1.1) will give an introduction into this research where this discrepancy was noticed. In chapter (2) we will take a closer look into the underlying mathematics, and in chapters (3) and (4) we will make a proposition on how to identify whether or not a researcher may have encountered the same problem, and how to select a new representation to avoid the discrepancy.
1.1
The clinical case
A human forearm contains two large bones, the radius and the ulna, see figure 1. The end part of the radius closest to the hand, the distal radius, plays an important role in the movement of the hand relative to the forearm. Fracture of the distal radius is a common injury as a result of falling on outstretched hands. Of all forearm injuries, up to 72% consists of fracture of the distal radius [1]. Treatment of this injury is divided into two main categories; casting the forearm or fixating the bones internally through surgery. However, up to 31% of these cases result in complications [3], [6]. The most common complication is malunion, where the bone segments did not heal in an optimal position. This results in divergent kinematics of the wrist joint, and therefor may result into serious impairment and pain. In this case, corrective surgery is to be considered.
Figure 1: Wrist anatomy, showing the two forearm bones; radius and ulna and two carpal bones; lunate and scaphoid. Adapted from G. Caiti [2].
To correct the malunion, the radius is cut at the height of the old fracture and the distal segment is repositioned, as presented in the schematic representation in figure 2.
Figure 2: Schematic representation of the radius. a) The malunited bone is cut at the height of the old fracture. b) The distal radius is repositioned to the right orientation. Adapted from G. Caiti [2].
1.1.1 Purpose
Researchers are constantly developing new methods to increase the accuracy and reduce the invasiveness of surgeries. In corrective osteotomy of the distal radius, patient specific guides are used to accurately transfer the virtual planning for the osteotomy to the patient. In most cases however, a large incision has to be made to optimally fit the guide to the patient. G. Caiti [2] designed and evaluated a patient-specific guiding system to transfer the virtually planned osteotomy of the distal radius in a minimally invasive fashion.
1.1.2 Methods
3D images of the malunited bone and the contralateral healthy radius are created. Using CT scans of both the patients forearms, a 3D image of each radius is created. These 3D images then are used to create polygonal meshes of the bones by segmentation by a region growing algorithm, as described by Dobbe et al. [5]. The distal and proximal ends of the effected ra-dius are then registered onto the mirrored image of the healthy rara-dius. The guiding system designed by G. Caiti [2], consists of a cast with two drilling slots and an external cutting guide. The patient’s arm is placed in the patient-specific cast and the self-tapping pins are be inserted in the two holes. The guide is slid over the pins after the cast is removed. See figure (3).
Reproducibility of this method was tested on five cadaver arm specimens. The cast con-tains a marker tool, representing the distal coordinate system, see figure (3). The marker tool consist of three aluminum spheres in orthogonal configuration. This coordinate system serves to quantify the positioning of the arm in the cast, to test the reproducibility of the arm place-ment. To quantify the position of the radius relative to the cast, each cadaver specimen was placed in, and removed from the cast 10 times. After each placement, a CT scan was acquired and a 3D polygonal model of the radius was registered to the CT scan. This yields 10 positions
of the arm with respect to the cast’s coordinate system. These 10 positions are referred to later in this study as 10 frames.
Figure 3: a) The arm is placed in the patient-specific cast. The cast features marker tools on the distal and proximal end (red circle), and two drilling slots (yellow circle); b) After removing the cast and inserting the percutaneous pins, the guide is slid over the pins. The cutting guide features a slit (yellow circle) to guide the cutting blade. Adapted from G. Caiti [2].
A Radius right-handed coordinate system (RCS) was defined for the radius in each frame, as seen in figure (4). The z-axis was chosen to match the central axis of the radius, the x-axis was directed toward the styloid process, perpendicular to the z-axis, and the y-axis was ori-ented to form an orthogonal right-handed triplet of vectors. Using the marker tools on the cast, the coordinate systems of the Distal radius (DCS) is defined. The relative displacement of the radius with respect to the cast was represented by the RCS-to-DCS transformation. The precision of this method is measured in the variability of the transformation from the RCS-to-DCS. Each transformation matrix yields information about six degrees of freedom (DOF), three translations and three rotations. The errors are reported as the standard deviation on the two sets of Euler angles calculated from each of the 10 rotation matrices.
Figure 4: After segmentation, three right-handed coordinate systems are defined. A coordinate system for the radius bone (RCS), and coordinate systems for the distal (DCS) marker tool. Adapted from G. Caiti [2].
1.1.3 Results
In the study, the standard deviations on the rotations from the rotation matrix of RCS-to-DCS were reported. First, however, the rotation matrices from DCS-to-RCS were examined. During calculations to retrieve the two sets of Euler angles, the researchers noted that the standard deviation on one of the rotations was disproportionately large. Figure (5) shows the coordinate systems of 2 out of the 10 frames, with respect to the global coordinate system. They are both DCS-to-RCS transformations. It is evident from image (5) that CS1 (the smaller multi color CS) and CS2 (the bigger multi color CS) are almost aligned. Still, the rotation angles that are calculated from their rotation matrices, lie far apart, see table (1). These Euler angles are calculated from the rotation matrix from rotation order ZXY. Note that, in the order of rotation ZXY, the first rotation about the Y axis, the second rotation about the X axis, and the final rotation about the Z axis. By convention, this order of rotation is written as ZXY. This will be thoroughly discussed in chapter (2).
Table 1: Both sets of Euler angles for the two coordinate systems CS1 and CS2, with respect to the global CS, for rotation sequence ZXY. Note that the values of φyand φzlie far apart.
ZXY CS1 to CS CS2 to CS
φx φy φz φx φy φz
Set 1 -79,080 158,066 121,957 -78,370 -152,958 172,193 Set 2 259,080 -21,934 -58,043 258,370 27,042 7,807
Figure 5: Three coordinate systems. The global CS represented in orange, the smaller multi color CS represents CS1, and the bigger multi color CS represents CS2. Note that CS1 and CS2 are almost aligned. CS1 and CS2 represent 2 frames from the series of 10 frames.
As shown in table (1), the ZXY rotation order from the DCS-to-RCS rotation matrix yields very different angles for φyand φz. Since the precision of the method is reported in terms of
this error, another rotation order was used to retrieve other Euler angles. Instead of acquiring the DCS-to-RCS rotation matrix, the inverse direction was was used (RCS-to-DCS) by taking the inverse of the rotation matrix before retrieving the Euler angles. The RCS-to-DCS rotation matrix yielded the following sets of Euler angles, with the corresponding errors are presented in table (2).
Table 2: Mean and Standard Deviation of both sets of Euler angles φx, φy, φz for rotation
sequence ZXY for the rotation matrix from the RCS-to-DCS transformation.
ZXY φx φy φz
set 1 Mean 54,200311 -108,845740 105,676688 Stdev 1,102709 1,798452 2,797610 set 2 Mean 125,799689 71,154260 -74,323312
Stdev 1,102709 1,798451 2,797610
It is evident from table (2) that the standard deviation on these Euler angles is low. Still, the reason for the disproportionally high standard deviation in the first original DCS-to-RCS transformation remained unknown.
1.2
Motivation
The aim of this study is to reproduce the approach used by G. Caiti [2] to calculate the Euler angles, and to clarify what caused the spread in the Euler angles calculated from the rotation sequence ZXY. Further, it is important to gain a good understanding of this phenomenon, to avoid the discrepancy in future experiments. To achieve this understanding, the relation between Euler angles and the chosen rotation sequence will be thoroughly elaborated on in chapter (2). Alternative rotation sequences within this methodology will be investigated and presented in chapters (2) and (3). A brief overview of the results will be provided in chapter (4), recommendations for further research will be suggested. Finally, in chapter (5), a summary will be presented.
2
Theoretical framework
There are different ways to describe the orientation of an object in three dimensions. In the field of biomedical engineering and physics, the orientation of an object is often represented by a matrix. The values in the matrix represent in which direction the three basis vectors point. These basis vectors are described using a reference CS, called the global CS. This yields a relative orientation of the object CS with respect to the global CS. In the following paragraphs it will be explained how to extract Euler angles from a rotation matrix and what challenges could be encountered during this process. [11]
2.1
Rotation matrices
A general rotation matrix could have the form,
M = M11 M12 M13 M21 M22 M23 M31 M32 M33 , (1)
where the columns of the matrix represent three orthogonal vectors. Matrices of this kind are part of the Special Orthogonal Group SO(3), where the Det(M) = +1 [13]. This rotation matrix is constructed from three rotations about the three cardinal axes. These rotations are called φxfor a rotation about the x axis, φy for a rotation about the y axis and finally φzfor a
rotation about the z axis. The matrix transformation about each axis is given by
Mx(φ) = 1 0 0 0 cos(φx) −sin(φx) 0 sin(φx) cos(φx) , My(φ) = cos(φy) 0 sin(φy) 0 1 0) −sin(φy) 0 cos(φy) , Mz(φ) = cos(φz) −sin(φz) 0 sin(φz) cos(φz) 0 0 0 1 . (2)
In the field of biomedical engineering and physics, The Visualization Toolkit (The VTK 7.0.0 in the study of G. Caiti [2]) is a commonly used toolkit to visualise scientific data [12]. In this program, the order of rotation is first about the y axis, then the x axis and finally about the z axis. For convenience, this order of rotation ZXY is used in the following example. Note that, in the rotation sequence ZXY, the first rotation about the Y axis is executed by the Y rotation matrix. In Linear Algebra, a matrix multiplication works on the matrix on it’s right. Hence, when the matrix ZXY works on an object in 3D space, the object is first rotated about the y axis, then about the x axis and finally about the z axis. As a convention, all sequences in this thesis are presented in this manner. See appendix A for the other rotation sequences. Combining these three rotation matrices yields the rotation matrix
By substituting equation (2) in equation (3), we find M = c(φy)c(φz) − s(φy)s(φx)s(φz) −s(φz)c(φx) s(φy)c(φz) + c(φy)s(φx)s(φz) c(φy)s(φz) + s(φy)s(φz)c(φx) c(φz)c(φx) s(φz)s(φy) − c(φy)s(φx)c(φz) −s(φy)c(φx) s(φx) c(φy)c(φx) , (4) where ”s” and ”c” are introduced as abbreviated notations for sine and cosine. With this matrix and object can be rotated in three dimensions.
2.2
Computing Euler angles from a rotation matrix
The angles mentioned in the previous paragraph, are called Euler angles. They represent a sequence of three rotations about the three axis to describe the orientation of an object. It is possible to subtract the Euler angles from a matrix with help of some mathematics [4]. From equation (4) it is evident that
M32= sin(φx)
φx= sin
−1
(M32).
(5) The sine function is symmetric around 12π so the other result for φxis given by
φx1= sin
−1
(M32)
φx2= π − φx1
(6) Computing φy is less straightforward; from (4) it is evident that M31 and M33 are both
dependent on φyand φx. By division, we will find the following relation;
tan(φy) =
−M31
M33
. (7)
Solving equation (7) for φy, there will be two options available depending on the sign from
cos(φx) that is hidden in M31and M33;
φy= atan2(−M31, M33) if cos(φx) > 0
φy= atan2(M31, −M33) if cos(φx) < 0
(8) Atan2 is an arctan function, used in coding languages like Python. The values of the atan2 function are described in equation (11). Making up for this sign in equation (8), it is easier to write these two options for φyas one,
φy= atan2( −M31 cos(φx1) , M33 cos(φx1) ) (9)
where division by cos(φx) solves the problem with the changing sign. Depending on the
two solutions for φx, the following two solutions for φyare yielded;
φy1= atan2( −M31 cos(φx1) , M33 cos(φx1) ) φy2= atan2( −M31 cos(φx2) , M33 cos(φx2) ). (10)
The atan2 function in equations (8-10) is defined as atan2(y, x) = arctan(yx) if x > 0, arctan(yx) + π if x < 0 and y ≥ 0, arctan(yx) − π if x < 0 and y < 0, −π 2 if x = 0 and y > 0, +π2 if x = 0 and y < 0, undefined if x = 0 and y = 0 (11)
and visualised in figure (6).
Figure 6: Plot of the atan2 function. Arctan(y/x) is plotted against y/x.
The calculations of φz follow the same steps as above. From (4) it is visible that M12 and
M22are both dependent on φzand φx. Division yields
tan(φz) =
−M12
M22
. (12)
Again, there is a cos(φx) hidden in M12and M22that gives two options
φz= atan2(−M12, M22) if cos(φx) > 0
φz= atan2(M12, −M22) if cos(φx) < 0.
Dividing by cos(φx) and filling in both solutions for φx, there are two solutions for φy; φz1= atan2( M12 cos(φx1), − M22 cos(φx1) ) φz2= atan2( M12 cos(φx2) , − M22 cos(φx2) ). (14)
For cos(φx) , 0 we found two solution sets, (φx1, φy1, φz1) and (φx2, φy2, φz2).
In the case φx= ±90°, and thus cos(φx) = 0, we encounter a problem. This would mean a
division by zero in equations 10 and 14. This occurrence is known as Gimbal Lock ([4]), where two angles effectively rotate about the same principal axis. In this case one degree of freedom is lost. In paragraph (2.3.4), this will be explained more in depth.
2.3
Challenges while using the Euler angle representation
Euler angles are a common representation to describe an orientation in the biomedical field, but are also widely used in the gaming and graphical designing industry. There are a couple of difficulties and different conventions that could be encountered while working with an Euler angle representation.
2.3.1 Tait-Bryan rotations
Often, when reading about Euler angle representations, it is not specified what type of Euler angle rotation is used. There are different types of Euler angle rotations. One type of Euler rotations, rotate about just two different axes. For example, a rotation about the x axis, followed by a rotation about the y axis, and finished by a rotation about the x axis again. These Euler rotations are calledClassic Euler rotations or Proper Euler rotations [11]. Rotations about three
different axis, are called Tait-Bryan Euler rotations. That is the kind of Euler rotations that are used in this thesis.
2.3.2 Intrinsic and extrinsic rotations
There is a difference in intrinsic and extrinsic rotations [11]. When using extrinsic rotations, an object in 3D is rotated about the axis of a global CS. This global CS does not move. When using intrinsic rotations, an object CS is used. Each rotation is performed on the CS as rotated by the former rotation(s). Euler rotations are intrinsic rotations.
2.3.3 Rotation axes
The difficulty in understanding the importance of the order or rotation, is the result of a con-ceptual problem. The axes are defined on the object, and as the rotations are executed, each following rotation also affects the axis of the former rotation(s) since intrinsic rotations are used. This is best visualized using an object with three coaxial rings as seen in figure (7), to represent the plane in which the three rotation work. The axis of rotation runs through the middle of the rings. The first rotation axis is represented by the innermost (blue) ring. The second rotation axis is oriented perpendicular to the first rotation axis, represented by the midsize (red) ring and the final rotation axis is perpendicular to the first two, represented by the outermost (green) ring. The first rotation rotates the object about the first axis. The sec-ond rotation rotates the object about the secsec-ond axis, and also rotates the first axis. The third
rotation rotates the object about the third axis and also rotates the other two axes. This is no problem when you rotate about each axis once, but might give rise to problems when multiple rotations about the same axis are performed. Note that the first rotation stays aligned with the object’s local axis, and that the third axis stays aligned with the global CS axis.
a) b) c)
Figure 7: a) The object in an orientation. The axes of rotation are visualised by three coaxial rings. b) Rotation about the first axis. Only the object is rotated. c) Rotation about the second axis. Note that the first axis is rotated together with the object. Images adapted from [8]
2.3.4 Gimbal Lock
Next to the problems that may ocur because of the rotation about the multiple axes themselves, there is a phenomenon called Gimbal Lock. This so-called Gimbal Lock occurs if the second axis of rotation is rotated over 90°[7]. As seen in figure (8), the coaxial rings representing the first and third rotation, align. This means that a rotation about the first and third axis has the same effect. In other words; one degree of freedom is lost.
a) b) c)
d) e) f)
Figure 8: a) The object in an orientation. The axis of rotation are visualised by three coaxial rings. b) Rotation about the first axis. c) Rotation about the second axis. d) When rotating the second axis 90°in total, Gimbal Lock occurs. e) and f): Rotating the first axis and the third axis have the same effect. One degree of freedom is lost. Images adapted from [8]
2.3.5 Rotation Order
To avoid running into a Gimbal Lock, another rotation order may be considered. Gimbal Lock occurs when the second axis rotates over 90 °. Also, other problems could occur, for example the discrepancy found by G. Caiti [2]. By selecting another sequence order, this problem may be avoided. There are six orders of rotation, ZYX, YZX, ZXY, XZY, YXZ, XYZ. The rotation matrices for those orders of rotation can be found in Appendix A. Since the Euler angles are calculated with different entries from the matrix, other Euler angles will be found for the different rotation sequences.
3
Experiments
The analysis of the data of G. Caiti [2] yield that different rotation sequences yield different values for the two sets of Euler angles, and their standard deviations. For the sequence ZXY, the following values for the mean and standard deviation for the two sets of Euler Angles are calculated, see table (3). The matrix used to obtain these Euler angles, is from the DCS-to-RCS transformation.
Table 3: Mean and Standard Deviation of both sets of Euler angles φz, φx, φy for rotation
sequence ZXY for the rotation matrix from the DCS-to-RCS transformation.
φz φx φy
set 1 Mean 147,002803 -78,717827 1,301234 Stdev 16,064603 0,660473 176,091266 set 2 Mean -32,997197 258,717827 1,301234
Stdev 16,064602 0,660473 16,386973
To visualise the individual measurements for the Euler angles in set 1 and set 2, plots of these values were created. Figure (9) shows the two sets of Euler angles for each of the 10 frames. For set 1, it is evident that for one of the rotation axes, the Euler angles flip it’s sign from positive to negative or vice versa at 180°. This results in a very high standard deviation for φy.
Figure 9: Both of Euler angles for rotation sequence ZXY for the DCS-to-RCS transformation. Each point represents one of the 10 frames. The average of the angles is given by the dashed line.
To provide indepth understanding why this angle changes from 180°to -180°, a plot of the atan2 function was created. In this plot, the arctan function is plotted versus it’s input, arctan(y/x) versus (y/x). The matrix entries that are used to calculate the specific Euler angle are represented by y and x respectively. The two Euler angles that are calculated with help of this function, namely φyand φz for rotation sequence ZXY, are plotted in figure (10). It is
evident that the values of φy1 change from one quadrant on the arctan curve to another. As
seen in figure (10), the values for the other angles are widely spread on the arctan curve. This results in a high standard deviation, and thus a low precision for this method. To resolve this
Figure 10: Atan2 function. Arctan(y/x) in radians against y/x, where y and x represent entries from the matrix. The values for φy1change from one quadrant on the arctan curve, to another.
The other angles are spread widely on the curve. For rotation sequence ZXY, in DCS-to-RCS transformation order.
The next step is evaluating what happened by changing the order of reference CS from DCS-to-RCS to RCS-to-DCS. Effectively, this change of the reference CS, RCS-to-DCS, will give the same standard deviations as the original DCS-to-RCS with inverse rotation sequence. To demonstrate this, the two orders DCS-to-RCS in ZXY and RCS-to-DCS in YXZ will be given next to each other in figures (11)-(13). As evident in figures (11) and (12), the plots are mirrored in the horizontal axis of the graph. The same entries from the original matrix, albeit inverted, are used to calculate these angles. To further prove this point, both arctan curves are plotted in figure (13). As seen in these plots, the angles have the same spread.
ZXY DCS-to-RCS YXZ RCS-to-DCS
Figure 11: First set of Euler angles for the two orders DCS-to-RCS in ZXY and RCS-to-DCS in YXZ. Each point represents one of the 10 frames. The average of the angles is represented by the dashed line.
ZXY DCS-to-RCS YXZ RCS-to-DCS
Figure 12: Second set of Euler angles for the two orders DCS-to-RCS in ZXY and RCS-to-DCS (called inverse matrix) in YXZ. Each point represents one of the 10 frames. The average of the angles is represented by the dashed line.
Table 4: Mean and Standard Deviation of both sets of Euler angles φz, φx, φy for rotation
sequence ZXY for the rotation matrix from the DCS-to-RCS transformation and for rotation sequence YXZ for the rotation matrix from the RCS-to-DCS transformation. Note that the standard deviations for the sets of the Euler angles are identical when comparing ZXY in DCS-to-RCS and YXZ in RCS-to-DCS
ZXY DCS-RCS φz φx φy YXZ RCS-DCS φy φx φz
set 1 Mean 147,003 -78,718 1,301 -1,301 78,718 -147,003
Stdev 16,065 0,660 176,091 176,091 0,660 16,065
set 2 Mean -32,997 258,718 1,301 -1,301 101,282 32,997
YXZ RCS-to-DCS ZXY DCS-to-RCS
Figure 13: Two atan2 functions for the two orders DCS-to-RCS in rotation sequence ZXY and RCS-to-DCS in rotation sequence YXZ. Arctan(y/x) in radians agains y/x, where y and x rep-resent entries from the matrix.
trix. Note, that the standard deviations for for the sets of the Euler angles are identical. This proves that for inverse matrix and reverse sequence, the same spread in the data point are found.
The next step is, to calculate all Euler angles for the other 5 rotation sequences. It is not necessary to take a look at the inverse order RCD-to-DCS matrices, as these give the same spread as the normal matrices in another rotation sequence. Tables (5)-(9) give an overview of all sets of Euler angles for the matrix from the DCS-to-RCS transformation.
Table 5: Mean and Standard Deviation of both sets of Euler angles φz, φy, φx for rotation
sequence ZYX for the rotation matrix from the DCS-to-RCS transformation.
ZYX φz φy φx
set 1 Mean -34,289084 -0,185308 -100,877176 Stdev 1,166839 3,152495 0,738839 set 2 Mean 145,710916 180,185308 79,122824
Stdev 1,166839 3,152495 0,738839
Table 6: Mean and Standard Deviation of both sets of Euler angles φy, φz, φx for rotation
sequence YZX for the rotation matrix from the DCS-to-RCS transformation. Note that the standard deviation for φyis extremely high in the second set of Euler angles.
YZX φy φz φx
set 1 Mean -0,226498 -34,235141 -101,008435 Stdev 3,826420 1,151790 1,962928 set 2 Mean -0,226498 214,235141 78,991565
Stdev 186,523610 1,151790 1,962929
Table 7: Mean and Standard Deviation of both sets of Euler angles φx, φz, φy for rotation
sequence XZY for the rotation matrix from the DCS-to-RCS transformation.
XZY φx φz φy
set 1 Mean -99,138556 5,946294 33,852873 Stdev 1,612507 2,741465 1,338404 set 2 Mean 80,861444 174,053706 -146,147127
Table 8: Mean and Standard Deviation of both sets of Euler angles φy, φx, φz for rotation
sequence YXZ for the rotation matrix from the DCS-to-RCS transformation.
YXZ φy φx φz
set 1 Mean 108,845722 -54,200302 -105,676695 Stdev 1,798464 1,102717 2,797613 set 2 Mean -71,154278 234,200302 74,323305
Stdev 1,798464 1,102717 2,797613
Table 9: Mean and Standard Deviation of both sets of Euler angles φx, φy, φz for rotation
sequence XYZ for the rotation matrix from the DCS-to-RCS transformation.
XYZ φx φy φz
set 1 Mean -103,072530 33,610403 7,125818 Stdev 0,705018 1,408177 3,252640 set 2 Mean 76,927470 146,389597 -172,874182
As shown in tables (5)-(9), one other rotation sequence has a high standard deviation. Ro-tation order YZX from table (6) has a high standard deviation in the second set of Euler angles, for angle φy. As appears from tables (4) and (6) neither of the second rotations are close to 90°,
and therefore the phenomenon Gimbal Lock is ruled out. The discrepancy is a consequence of the intrinsic behavior of Euler angles, to change quadrants in the arctan funcion when ex-ceding 180°. Therefore it is very important to evaluate all rotation orders, and select the best option.
The absolute length of vector spanned by the three standard deviations are determined per rotation order, see table (10).
Table 10: Table with absolute values of the vectors created by the standard deviations in each of the three directions x, y, and z. Only the rotation orders with the smaller standard deviations are accounted for.
ZYX XZY YXZ XYZ
Length of Stdev vector 3,441 3,450 3,503 3,613
G. Caiti selected the equivalent of rotation order YXZ from table (8) [2]. Table (10) shows an overview of the total standard deviation of the four options without the discrepancy. Evi-dently, the rotation order ZYX would give the lowest standard deviation in this case, followed by rotation order XYZ. This means that a suboptimal rotation order has been used in the study by G. Caiti [2].
4
Discussion
In this thesis, the data sets from G. Caiti [2] have been evaluated and the calculations for the Euler angles have been reproduced. An overview of all rotation sequences has been given. With this data set, the representation with the lowest total standard deviation would be se-quence ZYX (with a total standard deviation of 3,441), instead of the chosen equivalent of the YXZ rotation sequence (with a total standard deviation of 3,503). It is recommended to evalu-ate both sets of Euler angles for all rotation sequences, when using this representation for data sets.
There are a few remarks to be made.
1: The researcher, G. Caiti [2] found something peculiar, and instead of searching for the un-derlying mathematical cause, another method was tried without understanding the root cause of the problem. It is important to always try to understand the mathematical foundation of what has been calculated. In this case, the reported outcome was not the most precise option within the applied method. The reported standard deviation of the method was, in hindsight and with more attention for the mathematical background, smaller than reported in her study. 2: The order of the axes of the CS determine the order of the Euler angles. If the axes of the CS were named in a different order, the Euler angles that were found would have been assigned to a different order of rotation too. The chance that the researcher found Euler angle sets where this discrepancy was found, were 1 out of 3 because 2 out of 6 rotation orders showed a dis-crepancy.
3: The reporting of the precision in this method, namely based on the error on the calculated Euler angles, is not the good representation of the precision. The calculated angles depend on several aspects. The way the CS of the reference CS and the object CS are chosen, determine the outcome of the Euler angles. If one of these CSs had been chosen under a slight different angle, the discrepancy might even not have occurred at all. In other words: The angles that are found, and their errors, are highly dependent on the orientation of the two CSs. It would be more reliable to look for a method of reporting the precision that is only dependent on the method itself and not also on the orientation of the chosen CSs.
4: It is important to deliberately share knowledge between domains. Gimbal Lock is already widely known and described in other domains, for example the gaming and animation indus-try. It is also known that Euler angles tend to flip from 180°to -180°. This is also a recognised problem in physics, other representations are suggested such as axis angle or quaternions [9], [10].
5
Conclusion
The precision of the method used in the research by G. Caiti is expressed in terms of the standard deviation on the calculated Euler angles. These angles are retrieved from the rotation matrix from a DCS-to-RCS transformation, from rotation sequence ZXY. The wide spread in these values for the Euler angles are the cause of a high standard deviation. Analyzing the other five rotation sequences, proved fruitful. Four other rotation sequences (ZYX, XZY, YXZ and XYZ) yielded low standard deviations. The best rotation sequence would be sequence ZYX (with a total standard deviation of 3,441), in comparison with the chosen equivalent of rotation sequence YXZ (with a total standard deviation of 3,503). It is therefore concluded that, when using an Euler angle representation, all six rotation sequences should be analysed. It is important to actively share knowledge between different research domains. It is recommended to study the axis angle representation and the quaternion representation.
6
References
References
[1] B. D. Bushnell and D. K. Bynum. Malunion of the distal radius.J. Am. Acad. Orthop. Surg.,
15(1):27–40, 2007.
[2] Giuliana Caiti.Take it Personally. PhD thesis, University of Amsterdam, 2019.
[3] W. P. Cooney, J. H. Dobyns, and R. L. Linscheid. Complications of colles’ fractures.J. Bone Joint Surg. Am., 62(4):613–9, 1980.
[4] F. Dunn and I. Parberry. 3D Math Primer for Graphics and Game Development, chapter 8.
CRC Press, 2 edition, 2011.
[5] J. G. G. Dobbe et al. Computer-assisted planning and navigation for corrective distal radius osteotomy, based on pre- and intraoperative imaging. IEEE Trans. Biomed. Eng.,
58(1):182–190, 2011.
[6] J. Fischer, N. W. Thompson, and J. W. K. Harrison. Complications of colles fractures.
Classic Papers in Orthopaedics, 62(4):365–366, 2014.
[7] F. Sebastian Grassia. Practical parameterization of rotations using the exponential map.
The Journal of Graphics Tools, 3.3:29–48, 1998.
[8] GuerrillaCG. Euler (gimbal lock) explained. https://www.youtube.com/watch?v= zc8b2Jo7mno&t=147s, 2009. Last accessed 11 December 2019.
[9] Hashim A. Hashim. Special orthogonal group so(3), euler angles, angle-axis, rodriguez vector and unit-quaternion: Overview, mapping and challenges. arXiv:1909.06669v1,
2019.
[10] Dario Pavllo, Christoph Feichtenhofer, Michael Auli, and David Grangier. Modeling hu-man motion with quaternion-based neural networks. arXiv:1901.07677v2, 2019.
[11] D. Rose. Rotations in three-dimensions: Euler angles and rotation matri-ces. http://danceswithcode.net/engineeringnotes/rotations_in_3d/rotations_ in_3d_part1.html, 2015. Last accessed 19 December 2019.
[12] K. Martin W. Schroeder and B. Lorensen. The VTK: an object-oriented approach to 3D graphics. Pearson Education, Inc, 4 edition, 2006.
[13] G. ’t Hooft, B.Q.P.J. de Wit, and M.J.G. Veltman. Lie groups in physics (online course). 2007.
7
Acknowledgements
After years of following courses, a Bachelors student finishes with a grand finale in the form of an individual research of three months. In those months, you read about a subject, try to form an understanding about the subject, and do a small research of your own. I could not have accomplished this without the weekly meetings with Geert and Iwan. In the beginning, there was only a vague idea of the problem I was to tackle. By ping ponging back and forth every week about what steps would be logical, the problem slowly but surely became more clear. I want to thank Iwan as my daily supervisor for being available by email as well during the other days of the week, and for the constructive criticism on all the previous versions of my thesis. I want to thank Geert as my supervisor for the good ideas during our weekly meetings and for thinking of Gustav as my examiner. I want to thank Gustav for taking the time to attend my presentation and also for being my examiner.
8
Appendices
8.1
Appendix A
Apart from the Z-X-Y rotation order, there are five other rotation orders. All are build from the same three matrix transformations
Mx(φ) = 1 0 0 0 cos(φx) −sin(φx) 0 sin(φx) cos(φx) , My(φ) = cos(φy) 0 sin(φy) 0 1 0) −sin(φy) 0 cos(φy) , Mz(φ) = cos(φz) −sin(φz) 0 sin(φz) cos(φz) 0 0 0 1 . (15)
All rotation orders yield different Euler angles, since different entries from the matrix will be used to calculate these angles.
MZY X= c(φz)c(φy) c(φz)s(φy)s(φx) − c(φx)s(φz) c(φz)c(φx)s(φy) + s(φz)s(φx) c(φy)s(φz) c(φz)c(φx) + s(φz)s(φy)s(φx) −c(φz)s(φx) + c(φx)s(φz)s(φy) −s(φy) c(φy)s(φx) c(φy)c(φx) (16) MY ZX= c(φy)c(φz) −c(φx)c(φy)s(φz) + s(φx)s(φy) c(φx)s(φy) + c(φy)s(φx)s(φz) s(φz) c(φx)c(φz) −c(φz)s(φx) −c(φz)s(φy) c(φx)s(φy)s(φz) + c(φy)s(φx) c(φx)c(φy) − s(φx)s(φy)s(φz) (17) MZXY = c(φy)c(φz) − s(φy)s(φx)s(φz) −s(φz)c(φx) s(φy)c(φz) + c(φy)s(φx)s(φz) c(φy)s(φz) + s(φy)s(φz)c(φx) c(φz)c(φx) s(φz)s(φy) − c(φy)s(φx)c(φz) −s(φy)c(φx) s(φx) c(φy)c(φx) (18) MXZY = c(φy)c(φz) − s(φy)s(φx)s(φz) −s(φz)c(φx) s(φy)c(φz) + c(φy)s(φx)s(φz) c(φy)s(φz) + s(φy)s(φz)c(φx) c(φz)c(φx) s(φz)s(φy) − c(φy)s(φx)c(φz) −s(φy)c(φx) s(φx) c(φy)c(φx) (19) MY XZ= c(φy)c(φz) + s(φx)s(φy)s(φz) −c(φy)s(φz) + c(φz)s(φx)s(φy) c(φx)s(φy) c(φx)s(φz) c(φx)c(φz) −s(φx) c(φy)s(φx)s(φz) − c(φz)s(φy) c(φy)c(φz)s(φx) + s(φy)s(φz) c(φx)c(φy) (20) MXY Z= c(φy)c(φz) −c(φy)s(φz) s(φy) c(φx)s(φz) + c(φz)s(φx)s(φy) c(φx)c(φz) − s(φx)s(φy)s(φz) −c(φy)s(φx) −c(φx)c(φz)s(φy) + s(φx)s(φz) c(φx)s(φy)s(φz) + c(φz)s(φx) c(φx)c(φy) (21)
