Calculation of 2D Angles

from BIOMCH-L Dec 2002

Dear all,

Many years ago, Ray Smith and I developed a software package for
digitizing marker locations from digital movies of patients. We then
calculated the segment angles and derived the joint angles. It was
published here: Kirtley C & Smith RA (2001) Application of Multimedia to
the study of Human Movement. Multimedia Tools and Applications: 14 (3):
259-268.

I recall that at the time, Ray and I were fairly new to motion analysis
and we agonized for quite some time about how to cope with the
sectorization ('CAST') problem. For example, say you have a shank
segment that's almost vertical. The ATAN function of (Yknee -
Yankle)/(Xknee - Xankle) returns the correct angle of the shank when the
angle with the floor is less than 90 degrees, but returns a negative
value when it's at more than 90 degrees. For example, 120 degrees would
come out -30 degrees.

Eventually, Ray came up with a method of laboriously checking which
sector (quadrant) the segment is in and thereby correcting the angle.
This works fine, but I have always thought that there must be a more
elegant solution. I wonder if anyone can enlighten me?

I may as well take the opportunity to wish everyone Happy holidays...

Chris
--
Dr. Chris Kirtley MD PhD
Associate Professor
Dept. of Biomedical Engineering
Catholic University of America



I've got around this problem in the past by using 'if - then' statements.  For example code for
'if hip angle >180, then 180 + abs(hip angle)' where abs is the absolute value of the negative
angle will give a correct value.  Not sure if that helps in your programming format.

Tony

Anthony J. Blazevich, PhD.
Lecturer in Biomechanics
Department of Sport Sciences
Brunel University, Kingston Lane
Uxbridge UB8 3PH
UK



Chris,

I don't know if my solution is 'elegant' or not, but have a look at the
attached BASIC code.  The 'elegance' lies in the test of segment orientation
found in the gosub routine 'SegmentPosition' (I think).

For j = 1 to NSegments

  For i = InitFrameNumber to LastFrameNumber

    IP       = Iprox(j)
    ID       = Idist(j)
    Gosub SegmentPosition
  Next i

  HorFlag1 = 0:HorFlag2 = 0
  For i = InitFrameNumber to LastFrameNumber

    if i < LastFrameNumber then Gosub HorizontalTest
    theta(i,j) = THold(i) * raddeg

  Next i
  if HorFlag1 = 0 AND HorFlag2 = 0 then goto Continue1
  if HorFlag2 = 0 OR  HorFlag1 > HorFlag2 then
      goto Continue1
  else
      ISeg = j
      Gosub HorizontalFix
  end if

Continue1:

Next j

If ISFLag = 0 then goto SkipTheAboveSection

  if ISFlag > 20 then Note = 20

  For i = InitFrameNumber to LastFrameNumber

    if Note = 20 then

      IP       = ITempSeg(3,ISFlag-Note)
      ID       = ITempSeg(4,ISflag-Note)

    else

      IP       = ITempSeg(1,ISFlag)
      ID       = ITempSeg(2,ISflag)

    end if

    Gosub SegmentPosition

  Next i

  HorFlag1% = 0:HorFlag2% = 0
  For i = InitFrameNumber to LastFrameNumber

    if i < LastFrameNumber then Gosub HorizontalTest
    theta(i,0) = THold(i) * raddeg

  Next i
  if HorFlag1% = 0 AND HorFlag2% = 0 then goto SkipTheAboveSection
  if HorFlag2% = 0 OR  HorFlag1% > HorFlag2% then
      goto SkipTheAboveSection
  else
      ISeg = 0
      Gosub HorizontalFix
  end if

SkipTheAboveSection:

If IDAngleFlag = 0 then goto PrintMenu

For i = 1 to NJoints

  For j = 1 to NSegments

   If Prox$(i) = Segments$(j) then JProx(i) = j
   If Dist$(i) = Segments$(j) then JDist(i) = j

  Next j

Next i

If IDirection = -1 or IDirection = 1 then

  If ISFlag <> 0 then

    If ISFlag < 20 then JProx(ISFlag) = 0 else JDist(ISFlag-20) = 0

  end if

For j = 1 to NJoints

    For i = InitFrameNumber to LastFrameNumber

       rtheta(i,j) = 180 + (IDirection * (theta(i,JProx(j)) -_
                     theta(i,JDist(j)))) + Bias(j)

    Next i

  Next j

Else

  locate 10,4:Print _
  "WARNING, WARNING:  No angle measurement direction provided!!!"
  delay 1
  locate 10,4:Print _
  "                                                             "
  goto PrintMenu

End if

For j = 1 to NJoints
  For i = 2 to NumberOfFrames-1
     romega(i,j) = (rtheta(i+1,j) - rtheta(i-1,j)) * adegrad / TwiceTime
  Next i
Next j

For j = 1 to NJoints
  For i = 3 to NumberOfFrames-2
     ralpha(i,j) = (romega(i+1,j) - romega(i-1,j))/TwiceTime
  Next i
Next j

locate 10,4:Print "Relative Angle Kinematics calculated"

SegmentPosition:

    Rise     = y(i,ID) - y(i,IP)
    Rrun     = x(i,ID) - x(i,IP)
    THold(i) = 0.0

    if Rrun  < 0.0 then          ' Distal marker is behind proximal marker

          THold(i) = pi + ATN(Rise/Rrun)

    end if

    if Rrun  > 0.0 then          ' Distal marker is ahead of proximal marker

     if Rise > 0.0 then          ' Distal marker is higher than proximal

       THold(i) = ATN(Rise/Rrun)

     elseif Rise < 0.0 then      ' Distal marker is lower than proximal

       THold(i) = pi2 + ATN(Rise/Rrun)

     end if

    end if

    if Rrun = 0.0 then           ' Segment is vertical

       if Rise < 0 then          ' Distal marker is below proximal marker

          THold(i) = pi / 2 * -1

       end if

       if Rise > 0 then          ' Distal marker is above proximal marker

          THold(i) = pi / 2

       end if

    end if

return

HorizontalTest:

      DTest = THold(i+1)-THold(i)
      if ABS(DTest) > pi then
        if DTest < 0 then
          if HorFlag1% = 0 then HorFlag1% = i+1
                                              ' Segment crosses Horizontal from
          THold(i+1) = THold(i+1) + pi2       ' below to above
        else
          if HorFlag2% = 0 then HorFlag2% = i+1
                                              ' Segment crosses from above to
          THold(i+1) = THold(i+1) - pi2       ' below
        end if
      end if

return

HorizontalFix:

    For i = InitFrameNumber to LastFrameNumber

      theta(i,ISeg) = theta(i,ISeg) + 360.0

    Next i

return

Cheers,
Drew

Human Movement Laboratory
Clinical Research Centre
School of Health Professions
University of Brighton
49 Darley Road
Eastbourne, East Sussex
BN20 7UR
United Kingdom

Voice:  (01273) 644 166         [Int'l: +44 1273 644166]
FAX:            (01273) 643 652         [Int'l: +44 1273 643652]
Laboratory:     (01273) 643 767         [Int'l: +44 1273 643767]
http://www.brighton.ac.uk/sohp/index.htm



Dear Chris,

Sorry I missed the making the summary of replies, I
have only just read you email.

The problem is to calculate the segment angle from the
vector joining points (x1,y1) and (x2,y2) (the origins of
two adjacent segments), where the angle is relative to
the positive X axis. An alternative way avoiding using
atan and reducing the number of if() then {} statements
is:
a) define a unit vector from point 1 to point 2 called (x,y).
b) define theta1 = acos(x) and theta2 = asin(y).
c) if(theta2>=0.0) then {angle = theta1}
For calculating the angle as both positive and negative
rotations from the X axis (ie. +-180)
c) if(theta2<0.0) then {angle = (-1.0)*theta1}
For calculating the angle as positive only rotations from
the X axis (ie. +360)
c) if(theta2<0.0) then {angle = (360-theta1)}

Hope the helps

Allan Carman

School of Physiotherapy
University of Otago
Dunedin, NZ.



Dear Chris,

For atan(y/x), you can tell the quadrant where the angle lies by examining the
signs of y and x. For instance, if both y and x are negative, y/x will be in
the third quadrant. If one does not examine the signs of y and x, since y/x
will be positive, atan(y/x) will be mistaken to be in the first quadrant.
Of course, if you do atan(value), i.e. with no knowledge of x and y, there is
no way that you can tell the quadrant. Fortunately, in gait analysis, since
you derive x and y from the marker positions, eg (Yknee - Yankle)/(Xknee -
Xankle), you should have no problem.
In Matlab, there is actually a function atan2 which does the above
automatically. I prefer atan2 to atan so that we don't have to do manual
testing of the signs of x and y.
Does Ray Smith use the same method?

Happy Christmas to you too.

Raymond Lee



Chris,
As already mentioned, the atan2 function does what you want, in
Matlab, Excel etc.. If there are still angle jumps present, e.g. in the
case the angle goes over +/ - pi, you can in Matlab use the
function UNWRAP
>From the manual:
>>>>> unwrap
Unwrap phase angles.

Syntax
p = unwrap(p)
Description
p = unwrap(p) corrects the phase angles in vector p by adding
multiples of 2*pi where needed, to smooth the transitions across
branch cuts. When p is a matrix, unwrap corrects the phase
angles down each column. The phase must be in radians.

The unwrap function is part of the standard MATLAB language.

Limitations
unwrap tries to detect branch cut crossings, but it can be fooled by
sparse, rapidly changing phase values.

.
 

angle

Phase angle.

*******************************************************
At Hof
Institute of Human Movement Sciences &
Laboratory of Human Movement Analysis AZG
University of Groningen
A. Deusinglaan 1, room 321
postal address:
PO Box 196
NL-9700 AD GRONINGEN
THE NETHERLANDS
Tel:   (31) 50 363 2645
Fax:   (31) 50 363 3150
e-mail: a.l.hof@med.rug.nl
http://www.ppsw.rug.nl/~ibw/



HI Chris,
         The best solution that I've come across for this is the atan2
function in matlab. It does the quadrant checking for you.

Kevin J Deluzio <kevin.deluzio@dal.ca>



Dear Chris,
I am familiar with this problem. Many years ago I calculated a planar closed structure with one degree of freedom and made a
simulation. I calculated the so called "oriented angle" between the two links - the angle between the first link and the second
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
joint centre. I used arctg and arcsin  function to check in what quadrant is the angle. I have this program in Matlab code. If you
are interested I can send you this program or to explain in details the calculations and the formulas.
Merry Christmas and Happy New Year
Rosi Raikova
--------------------
Assoc. Prof. Rositsa Todorova Raikova
CENTRE OF BIOMEDICAL ENGINEERING
BULGARIAN ACADEMY OF SCIENCES
Acad.G.Bonchev str., Bl.105
1113 Sofia, Bulgaria
tel. (+3592) 70 05 27
fax: (+3592) 72 37 87
E-mail: Rosi.Raikova@clbme.bas.bg

Dear Dr. Kirtley,

Are you able to use the ATAN2 Function? The ATAN2 function is
implemented in most of the standard program languages (C, MatLab,
FORTRAN, &). This function is called with two arguments instead of only
one argument.
Use: alpha = ATAN2(Y, X) in the range ( Pi < alpha < Pi)
Instead of: alpha = ATAN(Z); Z = Y/X in the range ( Pi/2 < alpha < Pi/2)
However, internally the ATAN2 function does nothing else than your self
written tool: determining the particular quadrant from the signs of the
two arguments. There is definitely no other way to calculate it.

Best regards
Ulrich Simon

Dr.-Ing. Ulrich Simon
Institute of Orthopaedic Research and Biomechanics
University of Ulm
www.biomechanics.de



Dr Kirtley,

As you probably know, your problem stems from the fact that the
tan-function is periodic with a period of pi and thereby maps all real
numbers to all real numbers with periodicity pi. By definition atan,
the inverse of the tan-function, maps all real numbers to (-pi/2,pi/2).
If the argument of atan is negative, then it is mapped to (-pi/2,0)
and if it is positive then it is mapped to (0,pi/2).

Without loss of generalization, let the origin of the shank be in
coordinates 0,0 and the knee at x,y. The angle of the knee will then be
expressed as atan(y/x). y is always positive, whereas x can be both
positive and negative. For positive x, all is well in your case, but for
negative x, the fraction y/x becomes negative. atan does not distinguish
between the two cases of negative y/positive x and positive y/negative x
since all it sees is a single negative number. The atan-function then
"chooses" to map the negative number to the area (-pi/2,0). Note: atan
of course also doesn't distinguish between x and y both positive or both
negative, but that is not a problem in your case since y is always
positive.

The conclusion is that in order to be able to use the atan-function
properly, you have to know the sign of both y and x (i.e. the quadrant).
In this case it is quite simple since the only argument that needs to be
checked is x: If x is negative, then the value pi should be added to the
calculated angle.

However, due to this mapping problem, an alternative function has
arisen. It is commonly known as the atan2-function and it takes two
arguments instead of one: atan(y,x). This function will map all real
numbers to the area (-pi,pi). Its implementation is just as described
above: Calculate the arctangent of the ratio y/x and check the signs of
x and y to determine the quadrant. This function could just as easily
have been implemented to map all real numbers to (0,2pi).

Hope this helps.

Kind regards,
Einar S. Idsø

PhD student
Dept. of Engineering Cybernetics
Norwegian University of Science and Technology



Hi Chris,

Not sure what language you used, but many also have an ATAN2 function which
returns an angle between 0 and pi and 0 and -pi.  That solves the 90 degree
problem, so all you need to worry about is the 180 degree problem.  If you
calculate the angles in a local coordinate system, this is never a problem.

Jim Richards <jimr@UDel.Edu>



The ATAN2 function in Matlab eliminates the quadrature problem you are refering to...

Musolino, Mark <musolinomc@msx.upmc.edu>



Chris

If you define a horizontal unit vector (i) and a unit vector that points
along the longitudinal axis of the segment (q), the dot (scalar) product
function readily provides the angle between the two unit vectors.  The sign
of the Z component of the cross product of the two vectors then provides all
the information necessary to interpret the angle information.

For example, assume we have an orthogonal coordinate system with the X axis
horizontal and pointing to the right, the Y axis is vertical, and Z is
defined as X cross Y.  Unit vectors i and j are associated with the X and Y
axes, respectively.  Vector S is defined as pointing from the proximal to
the distal endpoint of the shank; thus q=S/magnitude(S).  The angle (theta)
between i and q is computed as:

Theta = acos(i @ q), where @ is the dot product operation.

Define a third vector, A, such that A = i X q, where X is the cross product
operation.

If the Z component of A is positive, you know that the distal endpoint has a
higher (greater) Z coordinate that the proximal endpoint, thus the segment
is located in either the first (theta < 90) or second (theta > 90) quadrant.
If the Z component of A is negative, you know that the distal endpoint has a
lower (smaller) Z coordinate that the proximal endpoint, thus the segment is
located in either the fourth (theta < 90) or third (theta > 90) quadrant.
At this point, assigning unique angle values based on quadrant location is
trivial.

Hope this helps,

Michael
 

________________________
Michael E. Feltner, Ph.D, FACSM
Dept. of Sports Medicine
Pepperdine University
Malibu, CA 90263 USA
EMAIL: michael.feltner@pepperdine.edu
WEB: http://faculty.pepperdine.edu/mfeltner/
VOICE: (310) 506-4312
FAX: (310) 506-4785

Here is the code in FORTRAN

      SUBROUTINE VECTP(V1,V2,V3)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC
C
C  SUBROUTINE VECTP COMPUTES THE CROSS PRODUCT OF VECTOR V1 CROSS VECTOR V2.
C  THE ANSWER IS RETURNED IN VECTOR V3.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION V1(3),V2(3),V3(3)
      V3(1)=(V1(2)*V2(3))-(V1(3)*V2(2))
      V3(2)=(V1(3)*V2(1))-(V1(1)*V2(3))
      V3(3)=(V1(1)*V2(2))-(V1(2)*V2(1))
      RETURN
      END

      SUBROUTINE SCALARP(Z1,Z2,DOTP)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC
C
C  SUBROUTINE SCALARP COMPUTES THE SCALAR PRODUCT OF VECTOR Z1 AND VECTOR
Z2.
C  THE ANSWER IS RETURNED IN DOTP.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC
      REAL*8 DOTP,Z1(3),Z2(3)
      INTEGER*4 K
      DOTP=0.0
      DO 1 K=1,3
        DOTP=DOTP+(Z1(K)*Z2(K))
1     CONTINUE
      RETURN
      END



About 90% of problems can be resolved by using ASIN rather than ATAN. This
returns a positive or negative value depending on whether the angle is one
side or the other of vertical as opposed to the ATAN function which returns
both as positive. It doesn't work beyond 90 degrees either way of  course
but few joints have a functional range of movement of more than 180 degrees
and thus by sensible choice of the reference line about which ASIN is
calculated for each joint/segment this problem can generally be avoided.
 

Richard Baker <richard.baker@rch.org.au>



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
function.

The following simple routine does it:

IF x not equal to zero  THEN angle = ATN(y/x)  ELSE angle = SGN(y) * pi / 2
IF x < 0  THEN angle = angle + pi

That's it.

Ziaul Hasan, Ph.D.
College of Applied Health Sciences
University of Illinois at Chicago -- Mail Code 898
1919 West Taylor Street, Room 447
Chicago, IL 60612-7251, U.S.A.
Phone:  Office 312-996-1504,  Lab 312-413-1137
               Fax     312-996-4583



Dear Hasan,

What about the thrid and fourth quadrants?

Chris
--
Dr. Chris Kirtley MD PhD
Associate Professor
Dept. of Biomedical Engineering
Catholic University of America


Dear Dr. Kirtley:

                      In the third quadrant, both x and y are negative, but because their ratio
                      is positive, ATN returns a positive angle (between 0 and pi/2). When pi is
                      added to it (based on x being negative), the result is between pi and
                      3pi/2, which is indeed in the third quadrant.

                      In the fourth quadrant, x is positive and y is negative. ATN returns a
                      negative angle (between 0 and -pi/2), which is correct.

                      For completeness' sake I might add that in the second quadrant x is
                      negative and y is positive. Thus ATN returns a negative angle (between 0
                      and -pi/2). And when pi is added to it (based on x being negative), the
                      result is between pi and +pi/2, which is indeed in the second quadrant.

                      Thus, the procedure returns the correct value in all quadrants.

                      However, the values returned vary between 0 and 3pi/2 for quadrants 1-3,
                      and between 0 and -pi/2 for quadrant 4. There is nothing wrong with this.
                      But if one prefers 0 to + - pi ranges, one can add a step in which, if the
                      angle is greater than pi, then 2pi is subtracted from it. Or if one prefers
                      the range from 0 to 2pi, one can simply ask that if the angle is negative,
                      then 2pi should be added to it.

                      Zia Hasan



Hi Chris,

First of all best wishes for 2003.

It was probably a month ago that you published your message about the
'CAST'-problem. At the time I did not get around to reacting. I may have
been missing something, but would it not make sense to try to use
;'imaginary number' functions, rather than trig-functions? I checked in
Mathcad and Matlab, they both have functions for calculating the argument
(=angle) of an imaginary number, and yield results between -pi and pi.
In MathCad: arg(z)
In Matlab: angle(z)

If you want numbers between 0 and 2pi: "if answer<0 then answer=answer+2pi"

I assume that using machine-coded routines will give optimal performance
(speed)

Greetings,

Edsko Hekman



Dear Chris,

Several different methods were posted but they do not have any difference in reality. It
basically comes down to the question of which quadrant the vector belongs to.

Before I get into my point, one thing I might point out is that Michael Feltner's method of
using i x q is no different from checking the sign of the sine function:

    Z-component of i x q =
        1) sin(theta)   (0 <= theta < pi )
        2) -sin(2pi - theta) = -sin(-theta) = sin(theta)        (pi <= theta < 2pi)

Thus

    Z-compone
nt of i x q = sin(theta)        (0 <= theta < 2pi)

Likewise, use of the dot product provides the cosine value:

    cos(theta) = ( i @ q ) / |i||q| = i @ q

I don't understand why you have to do the dot product and cross product here because the
components of the unit vector q readily give you the cosine and sine values:

    qx = cos(theta)
    qy = sin(theta)

OK, here is what I do. I have developed a C++ class called C2DPoint. This class has an
extensive set of vector functions including overloaded
operators. This class has two member variables:

class C2DPoint
{
public:
        float x;
        float y;
}

There are three functions that are related to the angle computation:

    float C2DPoint::AngPos();
    void C2DPoint::RotateBy( C2DPoint vc );
    C2DPoint C2DPoint::UnitVector();

For example, you have two vectors A and B on the XY-plane. If you want to compute the angular
position of vector A to the X-axis, it can be done as follows:

    C2DPoint A;
    // assign values to A here
    A =
A.UnitVector(); // make A a unit vector
    float angle = A.AngPos();   // 'angle' receives the angular position.

If you want to compute the relative angular position of A to B, it can be done as follows:

    C2DPoint A, B;
    // Assign values to A and B here
    A.RotateBy( B );    // This function first compputes the unit vector of B and rotate vector
A using B.
                        // Thus, the X-axis rotates to B.
    A = A.UnitVector();
    float angle = A.AngPos();

The AngPos() function is similar to ATAN2 b
ut has additional features:

float C2DPoint::AngPos()
{
        // x = cos, y = sin
        if ( x = -99999.f || y == -99999.f )    // if the vector is not initialized yet
                return -99999.f;

        if ( x == 0.f && y == 0.f )     // if both components are zeros
                return -99999.f;

        if ( x > 0.f )  // Quad I & IV
        {
                return (float)atan( y / x );
        }
        else if ( x < 0.f )
        {
                if ( y > 0.f )  // Quad II: move from Quad IV to II by adding pi
                        return (float)atan( y / x ) + 3.141592f;
                else            // Quad III: move from Quad
 I to Quad III by subtracting pi
                        return (float)atan( y / x ) - 3.141592f;
        }
        else    // 0.0
        {
                if ( y > 0.f )  // pi/2
                        return 3.141592f / 2.f;
                else            // -pi/2
                        return -3.141592f / 2.f;
        }
}

Note here that the ATAN function returns angle in the range of -pi/2 < theta < pi/2. The cosine
function is positive in both Quads I and IV which coincide with this range. If cosine is
negative, the vector is either in Quad II (sine > 0) or Quad III (sine < 0). If the cosine is 0
and the sine is not,
 the angle is either pi/2 (sine > 0) or -pi/2 (sine < 0). The AngPos() function returns a value
in the range of -pi <= theta < pi.

The RotateBy function is simple:

void C2DPoint::RotateBy( C2DPoint &vc )
{
        if ( x == -99999.f || y == -99999.f )   // the current vector is not initialized
                return;         // quit

        if ( vc.x == -99999.f || vc.y == -99999.f )     // vc is not initialized
                return;         // no rotation
        else
                vc = vc.UnitVector();   // make vc a unit vector

        // rotation
        C2DPoint tp = *this;    /
/ copy the current vector to a temporary vector
        // x  y         x = cos, y = sin
        // -y x
        x = vc.x * tp.x + vc.y * tp.y;          // rotation
        y = (-vc.y) * tp.x + vc.x * tp.y;
}

The UnitVector() function does not need any explanation:

C2DPoint C2DPoint::UnitVector()
{
        float mag = Magnitude();        // this function computes magnitude of the vector
        C2DPoint rtn;
        if ( mag != 0.f )
        {
                rtn.x = x / mag;
                rtn.y = y / mag;
        }
        return rtn;
}

Now, let me add one more issue to your original question: the
type of angle. There are two types of angles we need to compute in motion analysis: ANGULAR
DISTANCE between two vectors and RELATIVE ANGULAR POSITION of a vector to another. A good
example of the former is the joint angle while that for the latter is the joint motion
(flexion/extension) angle. For example, the knee joint angle can be computed from the shakn and
thigh vectors (proximal to distal). If you define the knee angle as the ANGULAR DISTANCE
between these two lines, you will always get an angle smal
ler than or equalt to pi rad. You will not be able to differentiate slightly flexed knee angle
and hyperextended knee angle. When you use the RELATIVE ANGULAR POSITION approach, you can
define the knee flexion angle as positive and knee hyperextension as negative.

Another function can be defined here: AngDistFrom().

float C2DPoint::AngDistFrom( C2DPoint & vc )
{
        // check if the vector is missing
        if ( x == -99999.f || y == -99999.f  || vc.x == -99999.f || vc.y == -99999.f )
                return -99999.f;
 

        if ( ( x == 0.f && y == 0.f )  || ( vc.x == 0.f && vc.y == 0.f ) )
                return -99999.f;

        // compute the unit vectors
        C2DPoint A = UnitVector();
        C2DPoint B = vc.UnitVector();

        // compute the angle using the acos function
        return acos( A.x * B.x + A.y * B.y );
}

The knee angle (angular distance between the thigh and the shank) can be computed as follows:

    C2DPoint A, B;      // A: thigh vector (hip->knee), B: shank vector (knee->ankle)
    // assign values to the vectors here
    A = A * (-1
.f);    // reverse the direction of A (knee->hip)
    float angdist = A.AngDistFrom( B ); // angular distance

On the other hand, the relative angular position of the shank vector to the thigh vector can be
defined as the knee flexion angle, assuming that the subject is facing to the right (+X):

    C2DPoint A, B;      // A: thigh vector (hip->knee), B: shank vector (knee->ankle)
    // assign values to the vectors here
    A.RotateBy( B );
    float angpos = A.AngPos();          // relative angular position

Po
sitive angle means the knee is flexed and negative angle means the knee is hyperextended.

2D angles can be defined as either angular distance or relative angular position of one to
another. The angle formed by two 3-D vectors must be defined as the angular distance unless the
vectors are projected to a local or global reference frame.

In summary, as far as computation of the angular position of a 2-D vector is concerned,
checking which quadrant the vector belongs to based on the cosine and sine functi
ons is the most straight and efficient method. Other methods will only complicate the
situation.

For more information, visit http://kwon3d.com and follow the links to

Theories > Volume II > User-Angle Issues
Theories > Volume II > Joint Kinematics > Computation of the Orientation Angles

Merry Christmas!

Young-Hoo Kwon
---------------------------------------------------------------------
- Young-Hoo Kwon, Ph.D.
- Associate Professor, Dept. of Kinesiology
- Director, Biomechanics Laboratory

- Texas Woman's University
- P.O. Box 425647
- Denton, TX 76204-5647
- Phone: (940) 898-2598
- Fax: (940) 898-2581
- Email: ykwon@twu.edu <mailto:ykwon@twu.edu>
- Homepage: http://kwon3d.com
- Korean kwon3d eGroup: http://kwon3d.com/korean/eGroup_kr.html
- Int'l kwon3d eGroup: http://kwon3d.com/eGroup_i.html



Dear Dr. Kirtley:

Happy New Year!

Perhaps it is too late to compete for your prize, but the algorithm can be even simpler using ATAN2:
==============
ATAN2(x,y) returns an angle between -pi and pi.
If angle < 0, angle = angle + 2 pi
==============

I have checked this algorithm with your raw data (XLS file attached).
 

Sincerely,
Zong-Ming Li, PhD
Musculoskeletal Research Center
Department of Orthopaedic Surgery
University of Pittsburgh
E1641 Biomedical Science Tower
210 Lothrop Street
Pittsburgh, PA 15213
Phone: 412 648 1494
Fax: 412 648 2001
Email: zmli@pitt.edu



Hi Chris:

I can't believe how much traffic this has generated...  It
seems that 2-D analysis is far from obsolete.

Anyway, maybe not all programming languages use the
same order of arguments in the atan2 function.  In Matlab,
Fortran, and C, the correct syntax would be

  angle = atan2(y,x)

Ton van den Bogert



Dear Chris,

Thank you for your reply.
The unit vector is required, I will endeavour to expand on the
method presented and to add some thoughts on the topic.

With a 2D vector a given by (x,y) and the axes by i(1,0) and
j(0,1). The angle between a and the positive i axis is given by:
cos(theta1) = a.i/||a||.||i|| = x/||a||
Similarly, The angle between a and the positive j axis is given by:
sin(theta2) = a.j/||a||.||i|| = y/||a||

with ||a|| the length of a, therefore
theta1 = acos(x/||a||)
theta2 = asin(y/||a||)

With 0<=theta1<=2Pi and -Pi<=theta2<=+Pi
For rotations relative to the positive  i axis, if theta2 is greater
than zero then a lies in the first or second quadrant, and
Angle = theta1
Otherwise theta2 is less than zero a lies in the third or fourth
quadrant and
Angle = (2.Pi -theta1)

Since sign(theta2) = sign(y), theta2 does not need to be
calculated and as we are only interested in the sign of y the
component y does not need to be normalized.
Therefore the approach can be simplified to:

x = x / ||(x,y)||;
angle = acos(x);
if( y < 0.0 ){ angle  = 2.0*Pi – angle; }

The method using tan

IF( x != 0.0 )
{ angle = ATN(y/x); }
else
{ angle = SGN(y) * pi/2.0; }
IF( x < 0 ){ angle = angle + pi; }
IF( angle < 0 ){ angle = angle + 2.0*pi; }

Does have the advantage that a unit vector is not required, but
requires extra if statements to check for division by zero, and in
which quadrant the vector lies.

A variation of the first cosine method is:

x = x / ||(x,y)||;
angle = acos(x);
if( y < 0.0 ){ angle = (-1.0)*angle; }

Which returns an angle between +Pi and -Pi, excluding -Pi

A variation on the first tan method was also given

IF( x != 0.0 )
{ angle = ATN(y/x); }
else
{ angle = SGN(y) * pi/2.0; }
IF( x < 0 ){ angle = angle + pi; }
If( angle > Pi ){ angle = angle - 2.0*pi; }

And has the same effect of returning an angle between +Pi and -
Pi, excluding -Pi

The function ATAN2(x,y) is functionally the same as the second
cosine and tan methods in returning and angle between +Pi and
–Pi, excluding -Pi., containing is own method to return the
required value.

The method

angle = ATAN2(x,y);
if(angle<0.0) { angle = angle + 2 pi; }

corrects the output from ATAN2(x,y) to the required angle
between 0 and +2Pi.

The next logical step is to define your own function ATAN3(x,y),
which could contain either the first cosine or tan methods to
return an angle between 0 and +2.Pi., the later correction would
not be required and only one function call would be needed.

I feel there is not a lot of difference between the methods, and I
think it is up to personal preference as to which is adopted.

Cheers

Allan Carman



Dear subscribers,

As this months's moderator, I feel it is now time to wrap up and
summarize this discussion.  It seems that everything has been said and
we should try to end this specific topic.  Below is a brief summary and
some comments on analogous issues in 3-D analysis.  Please feel free to
point out any errors, and further contributions on 3-D joint angles are
certainly welcome.  Who knows, we may even revive the "helical angle
wars" that took place on Biomch-L between 1989 and 1992.

Paolo de Leva's posting yesterday, elaborating on an earlier reply from
Michael Feltner, describes an elegant way to compute the angle
between two body segments, each defined by two markers.  This method
uses the fact that a dot product is proportional to the cosine of the
angle between two vectors.  Michael shared his Fortran code, and you
can find Matlab code for the same algorithm in the Kinemat
package, see http://www.isbweb.org/software/movanal/kinemat/angle2d.m
This technique allows a plus & minus 180 degree range of motion and
has the advantage that it also works in 3-D.

As said by Paolo, and in earlier contributions, ATAN2 can be used
to obtain the absolute orientation of a body segment.  It is then
a simple matter to subtract the orientations of two segments to get
the angle between them.  This will give the same result as the "dot
product" method, but can only be used with 2-D marker data.

Please note that there was an error in Chris Kirtley's last posting.
The order of arguments in the atan2 function should be atan2(y,x), and
not atan2(x,y)!  I think the y,x order came about because atan2(a,b) was
thought of as a 360-degree version of atan(a/b).  I verified the
parameter order for Matlab, Fortran, and C.  Please check your
documentation if you use other languages.

An advantage of using ATAN2, as opposed to ATAN, ASIN and ACOS, is that
the code is smaller but also that "if" statements are avoided.  As a general
rule, "if" is risky when applied to measured data; an "if" can amplify noise
when not used carefully.  Allan Carman's code using "IF( x != 0.0 )"
works, but one should realize that the "if" is really only used when x is
exactly zero, which never occurs in real-world data.  If you had a bug in
the "else" branch of such an "if", you would never know it, if you only
processed real data.  I have made such mistakes myself...  It is always
better to use algorithms where no "if" is needed.

Finally, I would like to mention that similar issues arise when
computing 3-D (cardanic) joint angles.  The usual procedure is to first
compute attitude matrices for each of the two segments, multiply one matrix
by the inverse of the other to obtain the relative attitude, and finally
extract Euler/Cardanic angles from this relative attitude matrix.
Young-Hoo Kwon described this last step on the Kwon3D site:
http://www.kwon3d.com/theory/euler/euler_angles.html

Unfortunately, few programmers/authors have taken advantage of the ATAN2
function.

Borrowing the example from Kwon3D, let's assume a 3x3 attitude matrix
with 3-D joint angles (theta, phi, psi) as follows:

      cos(th)cos(ps)    sin(ph)sin(th)cos(ps)+cos(ph)sin(ps)
-cos(ph)sin(th)cos(ps)+sin(ph)sin(ps)

T =  -cos(th)sin(ps)   -sin(ph)sin(th)sin(ps)+cos(ph)sin(ps)
cos(ph)sin(th)sin(ps)+sin(ph)cos(ps)

         sin(th)                  -sin(ph)cos(th)                        cos(ph)cos(th)

As in Kwon3D, we label these matrix elements as t11, t12, t13, etc.
Using such an attitude matrix, obtained from measurements, the joint
angles can be easily extracted as follows:

  ph = atan2(-t32,t33)

  th = atan2(t31,sqrt(t32*t32+t33*t33))

  ps = atan2(-t21,t11)

(..sqrt is the square root function)

Note that 4 of the 9 matrix elements were not used, but if T is a true
rotation matrix, i.e. orthonormal, using the other elements would give
identical results.

It is interesting to note that Herman Woltring's PRP.FORTRAN code, first
announced on Biomch-L in 1990, consistently uses ATAN2 whereas some more
recent authors (including Kwon3D and KineMat) do not.  Herman would have
had much to say about the current discussion, but as you know he died in
November 1992.

See:
http://biomch-l.isbweb.org
http://www.isbweb.org/software/movanal/prp.fortran

So I am glad that now the entire Biomch-L readership has rediscovered
this wonderful mathematical function ATAN2.  Thanks to Chris Kirtley
for raising the question!

If you are writing any kind of movement analysis software, be sure to
visit the movement analysis section of the ISB software pages:
http://www.isbweb.org/software/movanal.html .

An archive of the Biomch-L "helical angle wars" can be found at
http://biomch-l.isbweb.org

--

Ton van den Bogert, Biomch-L co-moderator

--

A.J. (Ton) van den Bogert, PhD
Department of Biomedical Engineering
Cleveland Clinic Foundation
9500 Euclid Avenue (ND-20)
Cleveland, OH 44195, USA
Phone/Fax: (216) 444-5566/9198



Dear Ton and all,

I'd like to point out a mistake in your summary.

> As in Kwon3D, we label these matrix elements as t11, t12, t13, etc.
> Using such an attitude matrix, obtained from measurements, the joint
> angles can be easily extracted as follows:
>
>   ph = atan2(-t32,t33)
>   th = atan2(t31,sqrt(t32*t32+t33*t33))
>   ps = atan2(-t21,t11)
> (..sqrt is the square root function)
>

ph = atan2(-t32,t33)

is not right because

atan2(-t32,t33) = atan2(sin(ph)cos(th),cos(ph)cos(th))
 

is not equivalent to

atan2(sin(ph),cos(ph)).

It depends on the sign of cos(th). If this is negative, you will have

atan2(-sin(ph),-cos(ph))

instead. Likewise,

atan2(t31,sqrt(t32*t32+t33*t33)) = atan2(sin(th),sqrt(cos(th)^2))

is not equivalent to

atan2(sin(th),cos(th))

because sqrt(a2) can be either +a or -a, etc.

I agree that the use of ATAN2 is relatively more convenient. I have used FORTRAN, QuickBASIC,
C, C++, and now Visual C++. The function has evolved accordingly and atan
2 was forgotten in the process somehow. I will update my function using atan2. (^_^)

Young-Hoo Kwon
------------------------------------------------------
- Young-Hoo Kwon, Ph.D.
- Biomechanics Lab, Texas Woman's University
- ykwon@twu.edu
- http://kwon3d.com



I think my equations are still OK.  These three equations produce
values for ph,th, and ps that, when plugged into the big equation
for T, give the attitude matrix that was measured.  Hence, these
angles are a correct description of the motion that was measured.
The problem that Young-Hoo mentions is related to uniqueness.

By using the sqrt() function when computing theta, I restrict
theta to a range between -90 and +90 degrees (since the second argument
of ATAN2 is positive).  This 180 degree range for the second rotation
is a common convention for cardanic angles, this is needed to obtain unique
angles.  If theta is not restricted like this, there would be two
sets of angles that produce the same attitude matrix, which is essentially
the problem described by Young-Hoo.

Once it is understood that theta is always between -90 and +90 degrees,
there is no longer a problem.

A.J. (Ton) van den Bogert, PhD



Dear Ton and all,

> By using the sqrt() function when computing theta, I restrict
> theta to a range between -90 and +90 degrees (since the second argument
> of ATAN2 is positive).  This 180 degree range for the second rotation
> is a common convention for cardanic angles, this is needed to
> obtain unique angles.

Restricting theta between -90 and +90 degrees is generally acceptable, but not in the shoulder.
For example, flex your shoulder by 45 degrees (phi = 45 degrees) and start abduction (incr
ease theta alone). Psi here is 0. When theta passes 90 degree point, your approach will produce
a phi of 225 (45 + 180 degrees) degrees since theta should be smaller than 90 degrees all the
time. This will also force psi to become 180 degrees. In this continuous motion, phi suddenly
jumps from 45 to 225 degrees and psi jumps from 0 to 180 degrees when the arm passes the 90
degree abduction point. Theta will increases to 90 degrees and then starts decreasing. This
kind of situation is common in the shoulder.
 You will see a similar case in the ring event of gymnastics.  A Psi of 180 degrees is not
acceptable because it means the shoulder is dislocated. (See Fig 3 on
http://kwon3d.com/theory/euler/orient.html for the illustration of the case.)

My approach will provide two sets of angles. But identifying the unique angle set is not
difficult. You just need to go one more step. By checking the discontinuity in phi and psi, you
can determine whether theta must be larger than 90 degrees or not. As I mentioned ear
lier, when theta must be larger than 90 degrees but we force it to be smaller, phi and psi
shows discontinuity by 180 degrees. The bottom line is that we have to consider the nature of
the joint motion.

This discontinuity by 180 degrees must not be confused with the regular angle discontinuity by
360 degrees. For example, if you describe the orientation of a gymnast in the air using the
orientation angles, phi gives us the somersault angle. Any airborne maneuvers that have more
than a full somersault wil
l result in a range of phi of larger than 360 degrees. The somersault angle must be continuous
but it can have discontinuities at 360 * n + 180 degrees (n = ..., -1, 0, 1, ...) if your angle
computation function provides the angles in the range of -180 to 180 degrees. The discontinuity
in the orientation angles by 360 degrees is due to the angle computation range while the
discontinuity by 180 degrees is because theta exceeds 90 degrees.

> I think my equations are still OK.  These three equations produce

> values for ph,th, and ps that, when plugged into the big equation
> for T, give the attitude matrix that was measured.  Hence, these
> angles are a correct description of the motion that was measured.

Checking T does not tell you that your choice of angles is OK. If you compute T from the
angles, it must be the one you started with. But the point is not the attitute, but the way to
reach the attitude of interest. We have two different ways to reach the same attitude or two
sets of orientation angle
s. (Both angle sets will give you the same attitude matrix.) Which way is the correct one is
the main issue. Again, we have to pay attention to the nature and continuity of the movement
and angles.

Young-Hoo Kwon