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
Tony
Anthony J. Blazevich,
PhD.
Lecturer in Biomechanics
Department of Sport Sciences
Brunel University, Kingston Lane
Uxbridge UB8 3PH
UK
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
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
School of Physiotherapy
University of Otago
Dunedin, NZ.
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.
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/
Kevin J Deluzio <kevin.deluzio@dal.ca>
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
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
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>
Musolino, Mark <musolinomc@msx.upmc.edu>
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
Richard Baker <richard.baker@rch.org.au>
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
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
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
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
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
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)
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
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
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
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
> 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.