<!doctype html public "-//w3c//dtd html 4.0 transitional//en"> <html> <head> <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1"> <meta name="GENERATOR" content="Mozilla/4.79 [en] (Windows NT 5.0; U) [Netscape]"> <meta name="Author" content="Dr. Chris Kirtley"> <title>2D Angle Calculation</title> </head> <body> <h2> Calculation of 2D Angles</h2> <h4> from BIOMCH-L Dec 2002</h4> Dear all, <p>Many years ago, Ray Smith and I developed a software package for <br>digitizing marker locations from digital movies of patients. We then <br>calculated the segment angles and derived the joint angles. It was <br>published here: Kirtley C &amp; Smith RA (2001) Application of Multimedia to <br>the study of Human Movement. Multimedia Tools and Applications: 14 (3): <br>259-268. <p>I recall that at the time, Ray and I were fairly new to motion analysis <br>and we agonized for quite some time about how to cope with the <br>sectorization ('CAST') problem. For example, say you have a shank <br>segment that's almost vertical. The ATAN function of (Yknee - <br>Yankle)/(Xknee - Xankle) returns the correct angle of the shank when the <br>angle with the floor is less than 90 degrees, but returns a negative <br>value when it's at more than 90 degrees. For example, 120 degrees would <br>come out -30 degrees. <p>Eventually, Ray came up with a method of laboriously checking which <br>sector (quadrant) the segment is in and thereby correcting the angle. <br>This works fine, but I have always thought that there must be a more <br>elegant solution. I wonder if anyone can enlighten me? <p>I may as well take the opportunity to wish everyone Happy holidays... <p>Chris <br>-- <br><a href="mailto:kirtleymd@yahoo.com">Dr. Chris Kirtley MD PhD</a> <br>Associate Professor <br>Dept. of Biomedical Engineering <br>Catholic University of America <br> <hr WIDTH="100%"> <br>I've got around this problem in the past by using 'if - then' statements.&nbsp; For example code for <br>'if hip angle >180, then 180 + abs(hip angle)' where abs is the absolute value of the negative <br>angle will give a correct value.&nbsp; Not sure if that helps in your programming format. <p>Tony <p><a href="mailto:Anthony.Blazevich@brunel.ac.uk">Anthony J. Blazevich, PhD.</a> <br>Lecturer in Biomechanics <br>Department of Sport Sciences <br>Brunel University, Kingston Lane <br>Uxbridge UB8 3PH <br>UK <hr WIDTH="100%"> <br>Chris, <p>I don't know if my solution is 'elegant' or not, but have a look at the <br>attached BASIC code.&nbsp; The 'elegance' lies in the test of segment orientation <br>found in the gosub routine 'SegmentPosition' (I think). <p>For j = 1 to NSegments <p>&nbsp; For i = InitFrameNumber to LastFrameNumber <p>&nbsp;&nbsp;&nbsp; IP&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = Iprox(j) <br>&nbsp;&nbsp;&nbsp; ID&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = Idist(j) <br>&nbsp;&nbsp;&nbsp; Gosub SegmentPosition <br>&nbsp; Next i <p>&nbsp; HorFlag1 = 0:HorFlag2 = 0 <br>&nbsp; For i = InitFrameNumber to LastFrameNumber <p>&nbsp;&nbsp;&nbsp; if i &lt; LastFrameNumber then Gosub HorizontalTest <br>&nbsp;&nbsp;&nbsp; theta(i,j) = THold(i) * raddeg <p>&nbsp; Next i <br>&nbsp; if HorFlag1 = 0 AND HorFlag2 = 0 then goto Continue1 <br>&nbsp; if HorFlag2 = 0 OR&nbsp; HorFlag1 > HorFlag2 then <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; goto Continue1 <br>&nbsp; else <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ISeg = j <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Gosub HorizontalFix <br>&nbsp; end if <p>Continue1: <p>Next j <p>If ISFLag = 0 then goto SkipTheAboveSection <p>&nbsp; if ISFlag > 20 then Note = 20 <p>&nbsp; For i = InitFrameNumber to LastFrameNumber <p>&nbsp;&nbsp;&nbsp; if Note = 20 then <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; IP&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = ITempSeg(3,ISFlag-Note) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ID&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = ITempSeg(4,ISflag-Note) <p>&nbsp;&nbsp;&nbsp; else <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; IP&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = ITempSeg(1,ISFlag) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ID&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = ITempSeg(2,ISflag) <p>&nbsp;&nbsp;&nbsp; end if <p>&nbsp;&nbsp;&nbsp; Gosub SegmentPosition <p>&nbsp; Next i <p>&nbsp; HorFlag1% = 0:HorFlag2% = 0 <br>&nbsp; For i = InitFrameNumber to LastFrameNumber <p>&nbsp;&nbsp;&nbsp; if i &lt; LastFrameNumber then Gosub HorizontalTest <br>&nbsp;&nbsp;&nbsp; theta(i,0) = THold(i) * raddeg <p>&nbsp; Next i <br>&nbsp; if HorFlag1% = 0 AND HorFlag2% = 0 then goto SkipTheAboveSection <br>&nbsp; if HorFlag2% = 0 OR&nbsp; HorFlag1% > HorFlag2% then <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; goto SkipTheAboveSection <br>&nbsp; else <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ISeg = 0 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Gosub HorizontalFix <br>&nbsp; end if <p>SkipTheAboveSection: <p>If IDAngleFlag = 0 then goto PrintMenu <p>For i = 1 to NJoints <p>&nbsp; For j = 1 to NSegments <p>&nbsp;&nbsp; If Prox$(i) = Segments$(j) then JProx(i) = j <br>&nbsp;&nbsp; If Dist$(i) = Segments$(j) then JDist(i) = j <p>&nbsp; Next j <p>Next i <p>If IDirection = -1 or IDirection = 1 then <p>&nbsp; If ISFlag &lt;> 0 then <p>&nbsp;&nbsp;&nbsp; If ISFlag &lt; 20 then JProx(ISFlag) = 0 else JDist(ISFlag-20) = 0 <p>&nbsp; end if <p>For j = 1 to NJoints <p>&nbsp;&nbsp;&nbsp; For i = InitFrameNumber to LastFrameNumber <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; rtheta(i,j) = 180 + (IDirection * (theta(i,JProx(j)) -_ <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; theta(i,JDist(j)))) + Bias(j) <p>&nbsp;&nbsp;&nbsp; Next i <p>&nbsp; Next j <p>Else <p>&nbsp; locate 10,4:Print _ <br>&nbsp; "WARNING, WARNING:&nbsp; No angle measurement direction provided!!!" <br>&nbsp; delay 1 <br>&nbsp; locate 10,4:Print _ <br>&nbsp; "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; " <br>&nbsp; goto PrintMenu <p>End if <p>For j = 1 to NJoints <br>&nbsp; For i = 2 to NumberOfFrames-1 <br>&nbsp;&nbsp;&nbsp;&nbsp; romega(i,j) = (rtheta(i+1,j) - rtheta(i-1,j)) * adegrad / TwiceTime <br>&nbsp; Next i <br>Next j <p>For j = 1 to NJoints <br>&nbsp; For i = 3 to NumberOfFrames-2 <br>&nbsp;&nbsp;&nbsp;&nbsp; ralpha(i,j) = (romega(i+1,j) - romega(i-1,j))/TwiceTime <br>&nbsp; Next i <br>Next j <p>locate 10,4:Print "Relative Angle Kinematics calculated" <p>SegmentPosition: <p>&nbsp;&nbsp;&nbsp; Rise&nbsp;&nbsp;&nbsp;&nbsp; = y(i,ID) - y(i,IP) <br>&nbsp;&nbsp;&nbsp; Rrun&nbsp;&nbsp;&nbsp;&nbsp; = x(i,ID) - x(i,IP) <br>&nbsp;&nbsp;&nbsp; THold(i) = 0.0 <p>&nbsp;&nbsp;&nbsp; if Rrun&nbsp; &lt; 0.0 then&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Distal marker is behind proximal marker <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; THold(i) = pi + ATN(Rise/Rrun) <p>&nbsp;&nbsp;&nbsp; end if <p>&nbsp;&nbsp;&nbsp; if Rrun&nbsp; > 0.0 then&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Distal marker is ahead of proximal marker <p>&nbsp;&nbsp;&nbsp;&nbsp; if Rise > 0.0 then&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Distal marker is higher than proximal <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; THold(i) = ATN(Rise/Rrun) <p>&nbsp;&nbsp;&nbsp;&nbsp; elseif Rise &lt; 0.0 then&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Distal marker is lower than proximal <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; THold(i) = pi2 + ATN(Rise/Rrun) <p>&nbsp;&nbsp;&nbsp;&nbsp; end if <p>&nbsp;&nbsp;&nbsp; end if <p>&nbsp;&nbsp;&nbsp; if Rrun = 0.0 then&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Segment is vertical <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if Rise &lt; 0 then&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Distal marker is below proximal marker <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; THold(i) = pi / 2 * -1 <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; end if <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if Rise > 0 then&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Distal marker is above proximal marker <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; THold(i) = pi / 2 <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; end if <p>&nbsp;&nbsp;&nbsp; end if <p>return <p>HorizontalTest: <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DTest = THold(i+1)-THold(i) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ABS(DTest) > pi then <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if DTest &lt; 0 then <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if HorFlag1% = 0 then HorFlag1% = i+1 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Segment crosses Horizontal from <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; THold(i+1) = THold(i+1) + pi2&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' below to above <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if HorFlag2% = 0 then HorFlag2% = i+1 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' Segment crosses from above to <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; THold(i+1) = THold(i+1) - pi2&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ' below <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; end if <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; end if <p>return <p>HorizontalFix: <p>&nbsp;&nbsp;&nbsp; For i = InitFrameNumber to LastFrameNumber <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; theta(i,ISeg) = theta(i,ISeg) + 360.0 <p>&nbsp;&nbsp;&nbsp; Next i <p>return <p>Cheers, <br>Drew <p>Human Movement Laboratory <br>Clinical Research Centre <br>School of Health Professions <br>University of Brighton <br>49 Darley Road <br>Eastbourne, East Sussex <br>BN20 7UR <br>United Kingdom <p>Voice:&nbsp; (01273) 644 166&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [Int'l: +44 1273 644166] <br>FAX:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; (01273) 643 652&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [Int'l: +44 1273 643652] <br>Laboratory:&nbsp;&nbsp;&nbsp;&nbsp; (01273) 643 767&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; [Int'l: +44 1273 643767] <br>http://www.brighton.ac.uk/sohp/index.htm <br> <hr WIDTH="100%"> <br>Dear Chris, <p>Sorry I missed the making the summary of replies, I <br>have only just read you email. <p>The problem is to calculate the segment angle from the <br>vector joining points (x1,y1) and (x2,y2) (the origins of <br>two adjacent segments), where the angle is relative to <br>the positive X axis. An alternative way avoiding using <br>atan and reducing the number of if() then {} statements <br>is: <br>a) define a unit vector from point 1 to point 2 called (x,y). <br>b) define theta1 = acos(x) and theta2 = asin(y). <br>c) if(theta2>=0.0) then {angle = theta1} <br>For calculating the angle as both positive and negative <br>rotations from the X axis (ie. +-180) <br>c) if(theta2&lt;0.0) then {angle = (-1.0)*theta1} <br>For calculating the angle as positive only rotations from <br>the X axis (ie. +360) <br>c) if(theta2&lt;0.0) then {angle = (360-theta1)} <p>Hope the helps <p><a href="mailto:acarman@gandalf.otago.ac.nz">Allan Carman</a> <p>School of Physiotherapy <br>University of Otago <br>Dunedin, NZ. <hr WIDTH="100%"> <br>Dear Chris, <p>For atan(y/x), you can tell the quadrant where the angle lies by examining the <br>signs of y and x. For instance, if both y and x are negative, y/x will be in <br>the third quadrant. If one does not examine the signs of y and x, since y/x <br>will be positive, atan(y/x) will be mistaken to be in the first quadrant. <br>Of course, if you do atan(value), i.e. with no knowledge of x and y, there is <br>no way that you can tell the quadrant. Fortunately, in gait analysis, since <br>you derive x and y from the marker positions, eg (Yknee - Yankle)/(Xknee - <br>Xankle), you should have no problem. <br>In Matlab, there is actually a function atan2 which does the above <br>automatically. I prefer atan2 to atan so that we don't have to do manual <br>testing of the signs of x and y. <br>Does Ray Smith use the same method? <p>Happy Christmas to you too. <p><a href="mailto:rsrlee@polyu.edu.hk">Raymond Lee</a> <hr WIDTH="100%"> <br>Chris, <br>As already mentioned, the atan2 function does what you want, in <br>Matlab, Excel etc.. If there are still angle jumps present, e.g. in the <br>case the angle goes over +/ - pi, you can in Matlab use the <br>function UNWRAP <br>>From the manual: <br>>>>>> unwrap <br>Unwrap phase angles. <p>Syntax <br>p = unwrap(p) <br>Description <br>p = unwrap(p) corrects the phase angles in vector p by adding <br>multiples of 2*pi where needed, to smooth the transitions across <br>branch cuts. When p is a matrix, unwrap corrects the phase <br>angles down each column. The phase must be in radians. <p>The unwrap function is part of the standard MATLAB language. <p>Limitations <br>unwrap tries to detect branch cut crossings, but it can be fooled by <br>sparse, rapidly changing phase values. <p>. <br>&nbsp; <p>angle <p>Phase angle. <p>******************************************************* <br>At Hof <br>Institute of Human Movement Sciences &amp; <br>Laboratory of Human Movement Analysis AZG <br>University of Groningen <br>A. Deusinglaan 1, room 321 <br>postal address: <br>PO Box 196 <br>NL-9700 AD GRONINGEN <br>THE NETHERLANDS <br>Tel:&nbsp;&nbsp; (31) 50 363 2645 <br>Fax:&nbsp;&nbsp; (31) 50 363 3150 <br>e-mail: a.l.hof@med.rug.nl <br>http://www.ppsw.rug.nl/~ibw/ <br> <hr WIDTH="100%"> <br>HI Chris, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; The best solution that I've come across for this is the atan2 <br>function in matlab. It does the quadrant checking for you. <p>Kevin J Deluzio &lt;kevin.deluzio@dal.ca> <hr WIDTH="100%"> <br>Dear Chris, <br>I am familiar with this problem. Many years ago I calculated a planar closed structure with one degree of freedom and made a <br>simulation. I calculated the so called "oriented angle" between the two links - the angle between the first link and the second <br>link between 0 and 360 degrees. Three points are need - one in the first link, one in the second link and the coordinates of the <br>joint centre. I used arctg and arcsin&nbsp; function to check in what quadrant is the angle. I have this program in Matlab code. If you <br>are interested I can send you this program or to explain in details the calculations and the formulas. <br>Merry Christmas and Happy New Year <br>Rosi Raikova <br>-------------------- <br>Assoc. Prof. Rositsa Todorova Raikova <br>CENTRE OF BIOMEDICAL ENGINEERING <br>BULGARIAN ACADEMY OF SCIENCES <br>Acad.G.Bonchev str., Bl.105 <br>1113 Sofia, Bulgaria <br>tel. (+3592) 70 05 27 <br>fax: (+3592) 72 37 87 <br>E-mail: Rosi.Raikova@clbme.bas.bg <hr WIDTH="100%"> <br>Dear Dr. Kirtley, <p>Are you able to use the ATAN2 Function? The ATAN2 function is <br>implemented in most of the standard program languages (C, MatLab, <br>FORTRAN, &amp;). This function is called with two arguments instead of only <br>one argument. <br>Use: alpha = ATAN2(Y, X) in the range ( Pi &lt; alpha &lt; Pi) <br>Instead of: alpha = ATAN(Z); Z = Y/X in the range ( Pi/2 &lt; alpha &lt; Pi/2) <br>However, internally the ATAN2 function does nothing else than your self <br>written tool: determining the particular quadrant from the signs of the <br>two arguments. There is definitely no other way to calculate it. <p>Best regards <br>Ulrich Simon <p>Dr.-Ing. Ulrich Simon <br>Institute of Orthopaedic Research and Biomechanics <br>University of Ulm <br>www.biomechanics.de <hr WIDTH="100%"> <br>Dr Kirtley, <p>As you probably know, your problem stems from the fact that the <br>tan-function is periodic with a period of pi and thereby maps all real <br>numbers to all real numbers with periodicity pi. By definition atan, <br>the inverse of the tan-function, maps all real numbers to (-pi/2,pi/2). <br>If the argument of atan is negative, then it is mapped to (-pi/2,0) <br>and if it is positive then it is mapped to (0,pi/2). <p>Without loss of generalization, let the origin of the shank be in <br>coordinates 0,0 and the knee at x,y. The angle of the knee will then be <br>expressed as atan(y/x). y is always positive, whereas x can be both <br>positive and negative. For positive x, all is well in your case, but for <br>negative x, the fraction y/x becomes negative. atan does not distinguish <br>between the two cases of negative y/positive x and positive y/negative x <br>since all it sees is a single negative number. The atan-function then <br>"chooses" to map the negative number to the area (-pi/2,0). Note: atan <br>of course also doesn't distinguish between x and y both positive or both <br>negative, but that is not a problem in your case since y is always <br>positive. <p>The conclusion is that in order to be able to use the atan-function <br>properly, you have to know the sign of both y and x (i.e. the quadrant). <br>In this case it is quite simple since the only argument that needs to be <br>checked is x: If x is negative, then the value pi should be added to the <br>calculated angle. <p>However, due to this mapping problem, an alternative function has <br>arisen. It is commonly known as the atan2-function and it takes two <br>arguments instead of one: atan(y,x). This function will map all real <br>numbers to the area (-pi,pi). Its implementation is just as described <br>above: Calculate the arctangent of the ratio y/x and check the signs of <br>x and y to determine the quadrant. This function could just as easily <br>have been implemented to map all real numbers to (0,2pi). <p>Hope this helps. <p>Kind regards, <br><a href="mailto:esi@itk.ntnu.no">Einar S. Ids&oslash;</a> <p>PhD student <br>Dept. of Engineering Cybernetics <br>Norwegian University of Science and Technology <hr WIDTH="100%"> <br>Hi Chris, <p>Not sure what language you used, but many also have an ATAN2 function which <br>returns an angle between 0 and pi and 0 and -pi.&nbsp; That solves the 90 degree <br>problem, so all you need to worry about is the 180 degree problem.&nbsp; If you <br>calculate the angles in a local coordinate system, this is never a problem. <p>Jim Richards &lt;jimr@UDel.Edu> <br> <hr WIDTH="100%"> <br>The ATAN2 function in Matlab eliminates the quadrature problem you are refering to... <p>Musolino, Mark &lt;musolinomc@msx.upmc.edu> <br> <hr WIDTH="100%"> <br>Chris <p>If you define a horizontal unit vector (i) and a unit vector that points <br>along the longitudinal axis of the segment (q), the dot (scalar) product <br>function readily provides the angle between the two unit vectors.&nbsp; The sign <br>of the Z component of the cross product of the two vectors then provides all <br>the information necessary to interpret the angle information. <p>For example, assume we have an orthogonal coordinate system with the X axis <br>horizontal and pointing to the right, the Y axis is vertical, and Z is <br>defined as X cross Y.&nbsp; Unit vectors i and j are associated with the X and Y <br>axes, respectively.&nbsp; Vector S is defined as pointing from the proximal to <br>the distal endpoint of the shank; thus q=S/magnitude(S).&nbsp; The angle (theta) <br>between i and q is computed as: <p>Theta = acos(i @ q), where @ is the dot product operation. <p>Define a third vector, A, such that A = i X q, where X is the cross product <br>operation. <p>If the Z component of A is positive, you know that the distal endpoint has a <br>higher (greater) Z coordinate that the proximal endpoint, thus the segment <br>is located in either the first (theta &lt; 90) or second (theta > 90) quadrant. <br>If the Z component of A is negative, you know that the distal endpoint has a <br>lower (smaller) Z coordinate that the proximal endpoint, thus the segment is <br>located in either the fourth (theta &lt; 90) or third (theta > 90) quadrant. <br>At this point, assigning unique angle values based on quadrant location is <br>trivial. <p>Hope this helps, <p>Michael <br>&nbsp; <p>________________________ <br>Michael E. Feltner, Ph.D, FACSM <br>Dept. of Sports Medicine <br>Pepperdine University <br>Malibu, CA 90263 USA <br>EMAIL: michael.feltner@pepperdine.edu <br>WEB: http://faculty.pepperdine.edu/mfeltner/ <br>VOICE: (310) 506-4312 <br>FAX: (310) 506-4785 <p>Here is the code in FORTRAN <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; SUBROUTINE VECTP(V1,V2,V3) <br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC <br>CCC <br>C <br>C&nbsp; SUBROUTINE VECTP COMPUTES THE CROSS PRODUCT OF VECTOR V1 CROSS VECTOR V2. <br>C&nbsp; THE ANSWER IS RETURNED IN VECTOR V3. <br>C <br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC <br>CCC <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; IMPLICIT REAL*8 (A-H,O-Z) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DIMENSION V1(3),V2(3),V3(3) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; V3(1)=(V1(2)*V2(3))-(V1(3)*V2(2)) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; V3(2)=(V1(3)*V2(1))-(V1(1)*V2(3)) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; V3(3)=(V1(1)*V2(2))-(V1(2)*V2(1)) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; RETURN <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; SUBROUTINE SCALARP(Z1,Z2,DOTP) <br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC <br>CCC <br>C <br>C&nbsp; SUBROUTINE SCALARP COMPUTES THE SCALAR PRODUCT OF VECTOR Z1 AND VECTOR <br>Z2. <br>C&nbsp; THE ANSWER IS RETURNED IN DOTP. <br>C <br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC <br>CCC <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; REAL*8 DOTP,Z1(3),Z2(3) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; INTEGER*4 K <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DOTP=0.0 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DO 1 K=1,3 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; DOTP=DOTP+(Z1(K)*Z2(K)) <br>1&nbsp;&nbsp;&nbsp;&nbsp; CONTINUE <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; RETURN <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; END <hr WIDTH="100%"> <br>About 90% of problems can be resolved by using ASIN rather than ATAN. This <br>returns a positive or negative value depending on whether the angle is one <br>side or the other of vertical as opposed to the ATAN function which returns <br>both as positive. It doesn't work beyond 90 degrees either way of&nbsp; course <br>but few joints have a functional range of movement of more than 180 degrees <br>and thus by sensible choice of the reference line about which ASIN is <br>calculated for each joint/segment this problem can generally be avoided. <br>&nbsp; <p>Richard Baker &lt;richard.baker@rch.org.au> <hr WIDTH="100%"> <br>Given x and y, we want to find arctan(y/x) in the correct quadrant. ATN is the stadard arctangent function, and SGN is the sign <br>function. <p>The following simple routine does it: <p>IF x not equal to zero&nbsp; THEN angle = ATN(y/x)&nbsp; ELSE angle = SGN(y) * pi / 2 <br>IF x &lt; 0&nbsp; THEN angle = angle + pi <p>That's it. <p>Ziaul Hasan, Ph.D. <br>College of Applied Health Sciences <br>University of Illinois at Chicago -- Mail Code 898 <br>1919 West Taylor Street, Room 447 <br>Chicago, IL 60612-7251, U.S.A. <br>Phone:&nbsp; Office 312-996-1504,&nbsp; Lab 312-413-1137 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Fax&nbsp;&nbsp;&nbsp;&nbsp; 312-996-4583 <hr WIDTH="100%"> <br>Dear Hasan, <p>What about the thrid and fourth quadrants? <p>Chris <br>-- <br>Dr. Chris Kirtley MD PhD <br>Associate Professor <br>Dept. of Biomedical Engineering <br>Catholic University of America <hr WIDTH="100%"> <p>Dear Dr. Kirtley: <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; In the third quadrant, both x and y are negative, but because their ratio <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; is positive, ATN returns a positive angle (between 0 and pi/2). When pi is <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; added to it (based on x being negative), the result is between pi and <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 3pi/2, which is indeed in the third quadrant. <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; In the fourth quadrant, x is positive and y is negative. ATN returns a <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; negative angle (between 0 and -pi/2), which is correct. <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; For completeness' sake I might add that in the second quadrant x is <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; negative and y is positive. Thus ATN returns a negative angle (between 0 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; and -pi/2). And when pi is added to it (based on x being negative), the <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; result is between pi and +pi/2, which is indeed in the second quadrant. <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Thus, the procedure returns the correct value in all quadrants. <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; However, the values returned vary between 0 and 3pi/2 for quadrants 1-3, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; and between 0 and -pi/2 for quadrant 4. There is nothing wrong with this. <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; But if one prefers 0 to + - pi ranges, one can add a step in which, if the <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; angle is greater than pi, then 2pi is subtracted from it. Or if one prefers <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; the range from 0 to 2pi, one can simply ask that if the angle is negative, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; then 2pi should be added to it. <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Zia Hasan <hr WIDTH="100%"> <br>Hi Chris, <p>First of all best wishes for 2003. <p>It was probably a month ago that you published your message about the <br>'CAST'-problem. At the time I did not get around to reacting. I may have <br>been missing something, but would it not make sense to try to use <br>;'imaginary number' functions, rather than trig-functions? I checked in <br>Mathcad and Matlab, they both have functions for calculating the argument <br>(=angle) of an imaginary number, and yield results between -pi and pi. <br>In MathCad: arg(z) <br>In Matlab: angle(z) <p>If you want numbers between 0 and 2pi: "if answer&lt;0 then answer=answer+2pi" <p>I assume that using machine-coded routines will give optimal performance <br>(speed) <p>Greetings, <p>Edsko Hekman <hr WIDTH="100%"> <br>Dear Chris, <p>Several different methods were posted but they do not have any difference in reality. It <br>basically comes down to the question of which quadrant the vector belongs to. <p>Before I get into my point, one thing I might point out is that Michael Feltner's method of <br>using i x q is no different from checking the sign of the sine function: <p>&nbsp;&nbsp;&nbsp; Z-component of i x q = <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 1) sin(theta)&nbsp;&nbsp; (0 &lt;= theta &lt; pi ) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 2) -sin(2pi - theta) = -sin(-theta) = sin(theta)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; (pi &lt;= theta &lt; 2pi) <p>Thus <p>&nbsp;&nbsp;&nbsp; Z-compone <br>nt of i x q = sin(theta)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; (0 &lt;= theta &lt; 2pi) <p>Likewise, use of the dot product provides the cosine value: <p>&nbsp;&nbsp;&nbsp; cos(theta) = ( i @ q ) / |i||q| = i @ q <p>I don't understand why you have to do the dot product and cross product here because the <br>components of the unit vector q readily give you the cosine and sine values: <p>&nbsp;&nbsp;&nbsp; qx = cos(theta) <br>&nbsp;&nbsp;&nbsp; qy = sin(theta) <p>OK, here is what I do. I have developed a C++ class called C2DPoint. This class has an <br>extensive set of vector functions including overloaded <br>operators. This class has two member variables: <p>class C2DPoint <br>{ <br>public: <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; float x; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; float y; <br>} <p>There are three functions that are related to the angle computation: <p>&nbsp;&nbsp;&nbsp; float C2DPoint::AngPos(); <br>&nbsp;&nbsp;&nbsp; void C2DPoint::RotateBy( C2DPoint vc ); <br>&nbsp;&nbsp;&nbsp; C2DPoint C2DPoint::UnitVector(); <p>For example, you have two vectors A and B on the XY-plane. If you want to compute the angular <br>position of vector A to the X-axis, it can be done as follows: <p>&nbsp;&nbsp;&nbsp; C2DPoint A; <br>&nbsp;&nbsp;&nbsp; // assign values to A here <br>&nbsp;&nbsp;&nbsp; A = <br>A.UnitVector(); // make A a unit vector <br>&nbsp;&nbsp;&nbsp; float angle = A.AngPos();&nbsp;&nbsp; // 'angle' receives the angular position. <p>If you want to compute the relative angular position of A to B, it can be done as follows: <p>&nbsp;&nbsp;&nbsp; C2DPoint A, B; <br>&nbsp;&nbsp;&nbsp; // Assign values to A and B here <br>&nbsp;&nbsp;&nbsp; A.RotateBy( B );&nbsp;&nbsp;&nbsp; // This function first compputes the unit vector of B and rotate vector <br>A using B. <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // Thus, the X-axis rotates to B. <br>&nbsp;&nbsp;&nbsp; A = A.UnitVector(); <br>&nbsp;&nbsp;&nbsp; float angle = A.AngPos(); <p>The AngPos() function is similar to ATAN2 b <br>ut has additional features: <p>float C2DPoint::AngPos() <br>{ <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // x = cos, y = sin <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( x = -99999.f || y == -99999.f )&nbsp;&nbsp;&nbsp; // if the vector is not initialized yet <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return -99999.f; <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( x == 0.f &amp;&amp; y == 0.f )&nbsp;&nbsp;&nbsp;&nbsp; // if both components are zeros <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return -99999.f; <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( x > 0.f )&nbsp; // Quad I &amp; IV <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; { <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return (float)atan( y / x ); <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; } <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else if ( x &lt; 0.f ) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; { <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( y > 0.f )&nbsp; // Quad II: move from Quad IV to II by adding pi <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return (float)atan( y / x ) + 3.141592f; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // Quad III: move from Quad <br>&nbsp;I to Quad III by subtracting pi <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return (float)atan( y / x ) - 3.141592f; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; } <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else&nbsp;&nbsp;&nbsp; // 0.0 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; { <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( y > 0.f )&nbsp; // pi/2 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return 3.141592f / 2.f; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // -pi/2 <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return -3.141592f / 2.f; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; } <br>} <p>Note here that the ATAN function returns angle in the range of -pi/2 &lt; theta &lt; pi/2. The cosine <br>function is positive in both Quads I and IV which coincide with this range. If cosine is <br>negative, the vector is either in Quad II (sine > 0) or Quad III (sine &lt; 0). If the cosine is 0 <br>and the sine is not, <br>&nbsp;the angle is either pi/2 (sine > 0) or -pi/2 (sine &lt; 0). The AngPos() function returns a value <br>in the range of -pi &lt;= theta &lt; pi. <p>The RotateBy function is simple: <p>void C2DPoint::RotateBy( C2DPoint &amp;vc ) <br>{ <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( x == -99999.f || y == -99999.f )&nbsp;&nbsp; // the current vector is not initialized <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // quit <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( vc.x == -99999.f || vc.y == -99999.f )&nbsp;&nbsp;&nbsp;&nbsp; // vc is not initialized <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // no rotation <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; vc = vc.UnitVector();&nbsp;&nbsp; // make vc a unit vector <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // rotation <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; C2DPoint tp = *this;&nbsp;&nbsp;&nbsp; / <br>/ copy the current vector to a temporary vector <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // x&nbsp; y&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; x = cos, y = sin <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // -y x <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; x = vc.x * tp.x + vc.y * tp.y;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // rotation <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; y = (-vc.y) * tp.x + vc.x * tp.y; <br>} <p>The UnitVector() function does not need any explanation: <p>C2DPoint C2DPoint::UnitVector() <br>{ <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; float mag = Magnitude();&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // this function computes magnitude of the vector <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; C2DPoint rtn; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( mag != 0.f ) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; { <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; rtn.x = x / mag; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; rtn.y = y / mag; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; } <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return rtn; <br>} <p>Now, let me add one more issue to your original question: the <br>type of angle. There are two types of angles we need to compute in motion analysis: ANGULAR <br>DISTANCE between two vectors and RELATIVE ANGULAR POSITION of a vector to another. A good <br>example of the former is the joint angle while that for the latter is the joint motion <br>(flexion/extension) angle. For example, the knee joint angle can be computed from the shakn and <br>thigh vectors (proximal to distal). If you define the knee angle as the ANGULAR DISTANCE <br>between these two lines, you will always get an angle smal <br>ler than or equalt to pi rad. You will not be able to differentiate slightly flexed knee angle <br>and hyperextended knee angle. When you use the RELATIVE ANGULAR POSITION approach, you can <br>define the knee flexion angle as positive and knee hyperextension as negative. <p>Another function can be defined here: AngDistFrom(). <p>float C2DPoint::AngDistFrom( C2DPoint &amp; vc ) <br>{ <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // check if the vector is missing <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( x == -99999.f || y == -99999.f&nbsp; || vc.x == -99999.f || vc.y == -99999.f ) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return -99999.f; <br>&nbsp; <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if ( ( x == 0.f &amp;&amp; y == 0.f )&nbsp; || ( vc.x == 0.f &amp;&amp; vc.y == 0.f ) ) <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return -99999.f; <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // compute the unit vectors <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; C2DPoint A = UnitVector(); <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; C2DPoint B = vc.UnitVector(); <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // compute the angle using the acos function <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return acos( A.x * B.x + A.y * B.y ); <br>} <p>The knee angle (angular distance between the thigh and the shank) can be computed as follows: <p>&nbsp;&nbsp;&nbsp; C2DPoint A, B;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // A: thigh vector (hip->knee), B: shank vector (knee->ankle) <br>&nbsp;&nbsp;&nbsp; // assign values to the vectors here <br>&nbsp;&nbsp;&nbsp; A = A * (-1 <br>.f);&nbsp;&nbsp;&nbsp; // reverse the direction of A (knee->hip) <br>&nbsp;&nbsp;&nbsp; float angdist = A.AngDistFrom( B ); // angular distance <p>On the other hand, the relative angular position of the shank vector to the thigh vector can be <br>defined as the knee flexion angle, assuming that the subject is facing to the right (+X): <p>&nbsp;&nbsp;&nbsp; C2DPoint A, B;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // A: thigh vector (hip->knee), B: shank vector (knee->ankle) <br>&nbsp;&nbsp;&nbsp; // assign values to the vectors here <br>&nbsp;&nbsp;&nbsp; A.RotateBy( B ); <br>&nbsp;&nbsp;&nbsp; float angpos = A.AngPos();&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // relative angular position <p>Po <br>sitive angle means the knee is flexed and negative angle means the knee is hyperextended. <p>2D angles can be defined as either angular distance or relative angular position of one to <br>another. The angle formed by two 3-D vectors must be defined as the angular distance unless the <br>vectors are projected to a local or global reference frame. <p>In summary, as far as computation of the angular position of a 2-D vector is concerned, <br>checking which quadrant the vector belongs to based on the cosine and sine functi <br>ons is the most straight and efficient method. Other methods will only complicate the <br>situation. <p>For more information, visit http://kwon3d.com and follow the links to <p>Theories > Volume II > User-Angle Issues <br>Theories > Volume II > Joint Kinematics > Computation of the Orientation Angles <p>Merry Christmas! <p>Young-Hoo Kwon <br>--------------------------------------------------------------------- <br>- Young-Hoo Kwon, Ph.D. <br>- Associate Professor, Dept. of Kinesiology <br>- Director, Biomechanics Laboratory <p>- Texas Woman's University <br>- P.O. Box 425647 <br>- Denton, TX 76204-5647 <br>- Phone: (940) 898-2598 <br>- Fax: (940) 898-2581 <br>- Email: ykwon@twu.edu &lt;<a href="mailto:ykwon@twu.edu">mailto:ykwon@twu.edu</a>> <br>- Homepage: http://kwon3d.com <br>- Korean kwon3d eGroup: <a target="_blank" href="http://kwon3d.com/korean/eGroup_kr.html">http://kwon3d.com/korean/eGroup_kr.html</a> <br>- Int'l kwon3d eGroup: http://kwon3d.com/eGroup_i.html <br> <hr WIDTH="100%"> <br>Dear Dr. Kirtley: <p>Happy New Year! <p>Perhaps it is too late to compete for your prize, but the algorithm can be even simpler using ATAN2: <br>============== <br>ATAN2(x,y) returns an angle between -pi and pi. <br>If angle &lt; 0, angle = angle + 2 pi <br>============== <p>I have checked this algorithm with your raw data (<a href="/faq/2Dangles.xls">XLS file attached</a>). <br>&nbsp; <p>Sincerely, <br>Zong-Ming Li, PhD <br>Musculoskeletal Research Center <br>Department of Orthopaedic Surgery <br>University of Pittsburgh <br>E1641 Biomedical Science Tower <br>210 Lothrop Street <br>Pittsburgh, PA 15213 <br>Phone: 412 648 1494 <br>Fax: 412 648 2001 <br>Email: <a href="mailto:zmli@pitt.edu">zmli@pitt.edu</a> <hr WIDTH="100%"> <br>Hi Chris: <p>I can't believe how much traffic this has generated...&nbsp; It <br>seems that 2-D analysis is far from obsolete. <p>Anyway, maybe not all programming languages use the <br>same order of arguments in the atan2 function.&nbsp; In Matlab, <br>Fortran, and C, the correct syntax would be <p>&nbsp; angle = atan2(y,x) <p><a href="mailto:bogert@bme.ri.ccf.org">Ton van den Bogert</a> <br> <hr WIDTH="100%"> <br>Dear Chris, <p>Thank you for your reply. <br>The unit vector is required, I will endeavour to expand on the <br>method presented and to add some thoughts on the topic. <p>With a 2D vector a given by (x,y) and the axes by i(1,0) and <br>j(0,1). The angle between a and the positive i axis is given by: <br>cos(theta1) = a.i/||a||.||i|| = x/||a|| <br>Similarly, The angle between a and the positive j axis is given by: <br>sin(theta2) = a.j/||a||.||i|| = y/||a|| <p>with ||a|| the length of a, therefore <br>theta1 = acos(x/||a||) <br>theta2 = asin(y/||a||) <p>With 0&lt;=theta1&lt;=2Pi and -Pi&lt;=theta2&lt;=+Pi <br>For rotations relative to the positive&nbsp; i axis, if theta2 is greater <br>than zero then a lies in the first or second quadrant, and <br>Angle = theta1 <br>Otherwise theta2 is less than zero a lies in the third or fourth <br>quadrant and <br>Angle = (2.Pi -theta1) <p>Since sign(theta2) = sign(y), theta2 does not need to be <br>calculated and as we are only interested in the sign of y the <br>component y does not need to be normalized. <br>Therefore the approach can be simplified to: <p>x = x / ||(x,y)||; <br>angle = acos(x); <br>if( y &lt; 0.0 ){ angle&nbsp; = 2.0*Pi  angle; } <p>The method using tan <p>IF( x != 0.0 ) <br>{ angle = ATN(y/x); } <br>else <br>{ angle = SGN(y) * pi/2.0; } <br>IF( x &lt; 0 ){ angle = angle + pi; } <br>IF( angle &lt; 0 ){ angle = angle + 2.0*pi; } <p>Does have the advantage that a unit vector is not required, but <br>requires extra if statements to check for division by zero, and in <br>which quadrant the vector lies. <p>A variation of the first cosine method is: <p>x = x / ||(x,y)||; <br>angle = acos(x); <br>if( y &lt; 0.0 ){ angle = (-1.0)*angle; } <p>Which returns an angle between +Pi and -Pi, excluding -Pi <p>A variation on the first tan method was also given <p>IF( x != 0.0 ) <br>{ angle = ATN(y/x); } <br>else <br>{ angle = SGN(y) * pi/2.0; } <br>IF( x &lt; 0 ){ angle = angle + pi; } <br>If( angle > Pi ){ angle = angle - 2.0*pi; } <p>And has the same effect of returning an angle between +Pi and - <br>Pi, excluding -Pi <p>The function ATAN2(x,y) is functionally the same as the second <br>cosine and tan methods in returning and angle between +Pi and <br> Pi, excluding -Pi., containing is own method to return the <br>required value. <p>The method <p>angle = ATAN2(x,y); <br>if(angle&lt;0.0) { angle = angle + 2 pi; } <p>corrects the output from ATAN2(x,y) to the required angle <br>between 0 and +2Pi. <p>The next logical step is to define your own function ATAN3(x,y), <br>which could contain either the first cosine or tan methods to <br>return an angle between 0 and +2.Pi., the later correction would <br>not be required and only one function call would be needed. <p>I feel there is not a lot of difference between the methods, and I <br>think it is up to personal preference as to which is adopted. <p>Cheers <p><a href="mailto:acarman@GANDALF.OTAGO.AC.NZ">Allan Carman</a> <hr WIDTH="100%"> <br>Dear subscribers, <p>As this months's moderator, I feel it is now time to wrap up and <br>summarize this discussion.&nbsp; It seems that everything has been said and <br>we should try to end this specific topic.&nbsp; Below is a brief summary and <br>some comments on analogous issues in 3-D analysis.&nbsp; Please feel free to <br>point out any errors, and further contributions on 3-D joint angles are <br>certainly welcome.&nbsp; Who knows, we may even revive the "helical angle <br>wars" that took place on Biomch-L between 1989 and 1992. <p>Paolo de Leva's posting yesterday, elaborating on an earlier reply from <br>Michael Feltner, describes an elegant way to compute the angle <br>between two body segments, each defined by two markers.&nbsp; This method <br>uses the fact that a dot product is proportional to the cosine of the <br>angle between two vectors.&nbsp; Michael shared his Fortran code, and you <br>can find Matlab code for the same algorithm in the Kinemat <br>package, see <a target="_blank" href="http://www.isbweb.org/software/movanal/kinemat/angle2d.m">http://www.isbweb.org/software/movanal/kinemat/angle2d.m</a> <br>This technique allows a plus &amp; minus 180 degree range of motion and <br>has the advantage that it also works in 3-D. <p>As said by Paolo, and in earlier contributions, ATAN2 can be used <br>to obtain the absolute orientation of a body segment.&nbsp; It is then <br>a simple matter to subtract the orientations of two segments to get <br>the angle between them.&nbsp; This will give the same result as the "dot <br>product" method, but can only be used with 2-D marker data. <p>Please note that there was an error in Chris Kirtley's last posting. <br>The order of arguments in the atan2 function should be atan2(y,x), and <br>not atan2(x,y)!&nbsp; I think the y,x order came about because atan2(a,b) was <br>thought of as a 360-degree version of atan(a/b).&nbsp; I verified the <br>parameter order for Matlab, Fortran, and C.&nbsp; Please check your <br>documentation if you use other languages. <p>An advantage of using ATAN2, as opposed to ATAN, ASIN and ACOS, is that <br>the code is smaller but also that "if" statements are avoided.&nbsp; As a general <br>rule, "if" is risky when applied to measured data; an "if" can amplify noise <br>when not used carefully.&nbsp; Allan Carman's code using "IF( x != 0.0 )" <br>works, but one should realize that the "if" is really only used when x is <br>exactly zero, which never occurs in real-world data.&nbsp; If you had a bug in <br>the "else" branch of such an "if", you would never know it, if you only <br>processed real data.&nbsp; I have made such mistakes myself...&nbsp; It is always <br>better to use algorithms where no "if" is needed. <p>Finally, I would like to mention that similar issues arise when <br>computing 3-D (cardanic) joint angles.&nbsp; The usual procedure is to first <br>compute attitude matrices for each of the two segments, multiply one matrix <br>by the inverse of the other to obtain the relative attitude, and finally <br>extract Euler/Cardanic angles from this relative attitude matrix. <br>Young-Hoo Kwon described this last step on the Kwon3D site: <br>http://www.kwon3d.com/theory/euler/euler_angles.html <p>Unfortunately, few programmers/authors have taken advantage of the ATAN2 <br>function. <p>Borrowing the example from Kwon3D, let's assume a 3x3 attitude matrix <br>with 3-D joint angles (theta, phi, psi) as follows: <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cos(th)cos(ps)&nbsp;&nbsp;&nbsp; sin(ph)sin(th)cos(ps)+cos(ph)sin(ps) <br>-cos(ph)sin(th)cos(ps)+sin(ph)sin(ps) <p>T =&nbsp; -cos(th)sin(ps)&nbsp;&nbsp; -sin(ph)sin(th)sin(ps)+cos(ph)sin(ps) <br>cos(ph)sin(th)sin(ps)+sin(ph)cos(ps) <p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; sin(th)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; -sin(ph)cos(th)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cos(ph)cos(th) <p>As in Kwon3D, we label these matrix elements as t11, t12, t13, etc. <br>Using such an attitude matrix, obtained from measurements, the joint <br>angles can be easily extracted as follows: <p>&nbsp; ph = atan2(-t32,t33) <p>&nbsp; th = atan2(t31,sqrt(t32*t32+t33*t33)) <p>&nbsp; ps = atan2(-t21,t11) <p>(..sqrt is the square root function) <p>Note that 4 of the 9 matrix elements were not used, but if T is a true <br>rotation matrix, i.e. orthonormal, using the other elements would give <br>identical results. <p>It is interesting to note that Herman Woltring's PRP.FORTRAN code, first <br>announced on Biomch-L in 1990, consistently uses ATAN2 whereas some more <br>recent authors (including Kwon3D and KineMat) do not.&nbsp; Herman would have <br>had much to say about the current discussion, but as you know he died in <br>November 1992. <p>See: <br><a target="_blank" href="http://biomch-l.isbweb.org">http://biomch-l.isbweb.org</a> <br><a target="_blank" href="http://www.isbweb.org/software/movanal/prp.fortran">http://www.isbweb.org/software/movanal/prp.fortran</a> <p>So I am glad that now the entire Biomch-L readership has rediscovered <br>this wonderful mathematical function ATAN2.&nbsp; Thanks to Chris Kirtley <br>for raising the question! <p>If you are writing any kind of movement analysis software, be sure to <br>visit the movement analysis section of the ISB software pages: <br><a target="_blank" href="http://www.isbweb.org/software/movanal.html">http://www.isbweb.org/software/movanal.html</a> . <p>An archive of the Biomch-L "helical angle wars" can be found at <br><a target="_blank" href="http://biomch-l.isbweb.org">http://biomch-l.isbweb.org</a> <p>-- <p>Ton van den Bogert, Biomch-L co-moderator <p>-- <p><a href="mailto:bogert@bme.ri.ccf.org">A.J. (Ton) van den Bogert, PhD</a> <br>Department of Biomedical Engineering <br>Cleveland Clinic Foundation <br>9500 Euclid Avenue (ND-20) <br>Cleveland, OH 44195, USA <br>Phone/Fax: (216) 444-5566/9198 <hr WIDTH="100%"> <br>Dear Ton and all, <p>I'd like to point out a mistake in your summary. <p>> As in Kwon3D, we label these matrix elements as t11, t12, t13, etc. <br>> Using such an attitude matrix, obtained from measurements, the joint <br>> angles can be easily extracted as follows: <br>> <br>>&nbsp;&nbsp; ph = atan2(-t32,t33) <br>>&nbsp;&nbsp; th = atan2(t31,sqrt(t32*t32+t33*t33)) <br>>&nbsp;&nbsp; ps = atan2(-t21,t11) <br>> (..sqrt is the square root function) <br>> <p>ph = atan2(-t32,t33) <p>is not right because <p>atan2(-t32,t33) = atan2(sin(ph)cos(th),cos(ph)cos(th)) <br>&nbsp; <p>is not equivalent to <p>atan2(sin(ph),cos(ph)). <p>It depends on the sign of cos(th). If this is negative, you will have <p>atan2(-sin(ph),-cos(ph)) <p>instead. Likewise, <p>atan2(t31,sqrt(t32*t32+t33*t33)) = atan2(sin(th),sqrt(cos(th)^2)) <p>is not equivalent to <p>atan2(sin(th),cos(th)) <p>because sqrt(a2) can be either +a or -a, etc. <p>I agree that the use of ATAN2 is relatively more convenient. I have used FORTRAN, QuickBASIC, <br>C, C++, and now Visual C++. The function has evolved accordingly and atan <br>2 was forgotten in the process somehow. I will update my function using atan2. (^_^) <p>Young-Hoo Kwon <br>------------------------------------------------------ <br>- Young-Hoo Kwon, Ph.D. <br>- Biomechanics Lab, Texas Woman's University <br>- <a href="mailto:ykwon@twu.edu">ykwon@twu.edu</a> <br>- <a target="_blank" href="http://kwon3d.com">http://kwon3d.com</a> <br> <hr WIDTH="100%"> <br>I think my equations are still OK.&nbsp; These three equations produce <br>values for ph,th, and ps that, when plugged into the big equation <br>for T, give the attitude matrix that was measured.&nbsp; Hence, these <br>angles are a correct description of the motion that was measured. <br>The problem that Young-Hoo mentions is related to uniqueness. <p>By using the sqrt() function when computing theta, I restrict <br>theta to a range between -90 and +90 degrees (since the second argument <br>of ATAN2 is positive).&nbsp; This 180 degree range for the second rotation <br>is a common convention for cardanic angles, this is needed to obtain unique <br>angles.&nbsp; If theta is not restricted like this, there would be two <br>sets of angles that produce the same attitude matrix, which is essentially <br>the problem described by Young-Hoo. <p>Once it is understood that theta is always between -90 and +90 degrees, <br>there is no longer a problem. <p><a href="mailto:bogert@bme.ri.ccf.org">A.J. (Ton) van den Bogert, PhD</a> <hr WIDTH="100%"> <br>Dear Ton and all, <p>> By using the sqrt() function when computing theta, I restrict <br>> theta to a range between -90 and +90 degrees (since the second argument <br>> of ATAN2 is positive).&nbsp; This 180 degree range for the second rotation <br>> is a common convention for cardanic angles, this is needed to <br>> obtain unique angles. <p>Restricting theta between -90 and +90 degrees is generally acceptable, but not in the shoulder. <br>For example, flex your shoulder by 45 degrees (phi = 45 degrees) and start abduction (incr <br>ease theta alone). Psi here is 0. When theta passes 90 degree point, your approach will produce <br>a phi of 225 (45 + 180 degrees) degrees since theta should be smaller than 90 degrees all the <br>time. This will also force psi to become 180 degrees. In this continuous motion, phi suddenly <br>jumps from 45 to 225 degrees and psi jumps from 0 to 180 degrees when the arm passes the 90 <br>degree abduction point. Theta will increases to 90 degrees and then starts decreasing. This <br>kind of situation is common in the shoulder. <br>&nbsp;You will see a similar case in the ring event of gymnastics.&nbsp; A Psi of 180 degrees is not <br>acceptable because it means the shoulder is dislocated. (See Fig 3 on <br>http://kwon3d.com/theory/euler/orient.html for the illustration of the case.) <p>My approach will provide two sets of angles. But identifying the unique angle set is not <br>difficult. You just need to go one more step. By checking the discontinuity in phi and psi, you <br>can determine whether theta must be larger than 90 degrees or not. As I mentioned ear <br>lier, when theta must be larger than 90 degrees but we force it to be smaller, phi and psi <br>shows discontinuity by 180 degrees. The bottom line is that we have to consider the nature of <br>the joint motion. <p>This discontinuity by 180 degrees must not be confused with the regular angle discontinuity by <br>360 degrees. For example, if you describe the orientation of a gymnast in the air using the <br>orientation angles, phi gives us the somersault angle. Any airborne maneuvers that have more <br>than a full somersault wil <br>l result in a range of phi of larger than 360 degrees. The somersault angle must be continuous <br>but it can have discontinuities at 360 * n + 180 degrees (n = ..., -1, 0, 1, ...) if your angle <br>computation function provides the angles in the range of -180 to 180 degrees. The discontinuity <br>in the orientation angles by 360 degrees is due to the angle computation range while the <br>discontinuity by 180 degrees is because theta exceeds 90 degrees. <p>> I think my equations are still OK.&nbsp; These three equations produce <p>> values for ph,th, and ps that, when plugged into the big equation <br>> for T, give the attitude matrix that was measured.&nbsp; Hence, these <br>> angles are a correct description of the motion that was measured. <p>Checking T does not tell you that your choice of angles is OK. If you compute T from the <br>angles, it must be the one you started with. But the point is not the attitute, but the way to <br>reach the attitude of interest. We have two different ways to reach the same attitude or two <br>sets of orientation angle <br>s. (Both angle sets will give you the same attitude matrix.) Which way is the correct one is <br>the main issue. Again, we have to pay attention to the nature and continuity of the movement <br>and angles. <p>Young-Hoo Kwon <br><a href="mailto:bogert@bme.ri.ccf.org"></a>&nbsp; </body> </html>