ÿþ<!doctype html public "//w3c//dtd html 4.0 transitional//en">
<html>
<head>
<meta httpequiv="ContentType" content="text/html; charset=iso88591">
<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 BIOMCHL 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 & Smith RA (2001) Application of Multimedia
to
<br>the study of Human Movement. Multimedia Tools and Applications: 14
(3):
<br>259268.
<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.
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. 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. 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> For i = InitFrameNumber to LastFrameNumber
<p> IP = Iprox(j)
<br> ID = Idist(j)
<br> Gosub SegmentPosition
<br> Next i
<p> HorFlag1 = 0:HorFlag2 = 0
<br> For i = InitFrameNumber to LastFrameNumber
<p> if i < LastFrameNumber then Gosub HorizontalTest
<br> theta(i,j) = THold(i) * raddeg
<p> Next i
<br> if HorFlag1 = 0 AND HorFlag2 = 0 then goto Continue1
<br> if HorFlag2 = 0 OR HorFlag1 > HorFlag2 then
<br> goto Continue1
<br> else
<br> ISeg = j
<br> Gosub HorizontalFix
<br> end if
<p>Continue1:
<p>Next j
<p>If ISFLag = 0 then goto SkipTheAboveSection
<p> if ISFlag > 20 then Note = 20
<p> For i = InitFrameNumber to LastFrameNumber
<p> if Note = 20 then
<p> IP
= ITempSeg(3,ISFlagNote)
<br> ID
= ITempSeg(4,ISflagNote)
<p> else
<p> IP
= ITempSeg(1,ISFlag)
<br> ID
= ITempSeg(2,ISflag)
<p> end if
<p> Gosub SegmentPosition
<p> Next i
<p> HorFlag1% = 0:HorFlag2% = 0
<br> For i = InitFrameNumber to LastFrameNumber
<p> if i < LastFrameNumber then Gosub HorizontalTest
<br> theta(i,0) = THold(i) * raddeg
<p> Next i
<br> if HorFlag1% = 0 AND HorFlag2% = 0 then goto SkipTheAboveSection
<br> if HorFlag2% = 0 OR HorFlag1% > HorFlag2% then
<br> goto SkipTheAboveSection
<br> else
<br> ISeg = 0
<br> Gosub HorizontalFix
<br> end if
<p>SkipTheAboveSection:
<p>If IDAngleFlag = 0 then goto PrintMenu
<p>For i = 1 to NJoints
<p> For j = 1 to NSegments
<p> If Prox$(i) = Segments$(j) then JProx(i) = j
<br> If Dist$(i) = Segments$(j) then JDist(i) = j
<p> Next j
<p>Next i
<p>If IDirection = 1 or IDirection = 1 then
<p> If ISFlag <> 0 then
<p> If ISFlag < 20 then JProx(ISFlag) = 0 else JDist(ISFlag20)
= 0
<p> end if
<p>For j = 1 to NJoints
<p> For i = InitFrameNumber to LastFrameNumber
<p> rtheta(i,j) = 180 + (IDirection
* (theta(i,JProx(j)) _
<br>
theta(i,JDist(j)))) + Bias(j)
<p> Next i
<p> Next j
<p>Else
<p> locate 10,4:Print _
<br> "WARNING, WARNING: No angle measurement direction provided!!!"
<br> delay 1
<br> locate 10,4:Print _
<br> "
"
<br> goto PrintMenu
<p>End if
<p>For j = 1 to NJoints
<br> For i = 2 to NumberOfFrames1
<br> romega(i,j) = (rtheta(i+1,j)  rtheta(i1,j))
* adegrad / TwiceTime
<br> Next i
<br>Next j
<p>For j = 1 to NJoints
<br> For i = 3 to NumberOfFrames2
<br> ralpha(i,j) = (romega(i+1,j)  romega(i1,j))/TwiceTime
<br> Next i
<br>Next j
<p>locate 10,4:Print "Relative Angle Kinematics calculated"
<p>SegmentPosition:
<p> Rise = y(i,ID)  y(i,IP)
<br> Rrun = x(i,ID)  x(i,IP)
<br> THold(i) = 0.0
<p> if Rrun < 0.0 then
' Distal marker is behind proximal marker
<p> THold(i) = pi
+ ATN(Rise/Rrun)
<p> end if
<p> if Rrun > 0.0 then
' Distal marker is ahead of proximal marker
<p> if Rise > 0.0 then
' Distal marker is higher than proximal
<p> THold(i) = ATN(Rise/Rrun)
<p> elseif Rise < 0.0 then
' Distal marker is lower than proximal
<p> THold(i) = pi2 + ATN(Rise/Rrun)
<p> end if
<p> end if
<p> if Rrun = 0.0 then
' Segment is vertical
<p> if Rise < 0 then
' Distal marker is below proximal marker
<p> THold(i) = pi
/ 2 * 1
<p> end if
<p> if Rise > 0 then
' Distal marker is above proximal marker
<p> THold(i) = pi
/ 2
<p> end if
<p> end if
<p>return
<p>HorizontalTest:
<p> DTest = THold(i+1)THold(i)
<br> if ABS(DTest) > pi then
<br> if DTest < 0 then
<br> if HorFlag1%
= 0 then HorFlag1% = i+1
<br>
' Segment crosses Horizontal from
<br> THold(i+1) =
THold(i+1) + pi2 ' below to above
<br> else
<br> if HorFlag2%
= 0 then HorFlag2% = i+1
<br>
' Segment crosses from above to
<br> THold(i+1) =
THold(i+1)  pi2 ' below
<br> end if
<br> end if
<p>return
<p>HorizontalFix:
<p> For i = InitFrameNumber to LastFrameNumber
<p> theta(i,ISeg) = theta(i,ISeg) + 360.0
<p> 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: (01273) 644 166
[Int'l: +44 1273 644166]
<br>FAX:
(01273) 643 652 [Int'l:
+44 1273 643652]
<br>Laboratory: (01273) 643 767
[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<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<0.0) then {angle = (360theta1)}
<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>
<p>angle
<p>Phase angle.
<p>*******************************************************
<br>At Hof
<br>Institute of Human Movement Sciences &
<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>NL9700 AD GRONINGEN
<br>THE NETHERLANDS
<br>Tel: (31) 50 363 2645
<br>Fax: (31) 50 363 3150
<br>email: a.l.hof@med.rug.nl
<br>http://www.ppsw.rug.nl/~ibw/
<br>
<hr WIDTH="100%">
<br>HI Chris,
<br> 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 <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 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>Email: 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, &). This function is called with two arguments instead
of only
<br>one argument.
<br>Use: alpha = ATAN2(Y, X) in the range ( Pi < alpha < Pi)
<br>Instead of: alpha = ATAN(Z); Z = Y/X in the range ( Pi/2 < alpha
< 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>tanfunction 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 tanfunction, 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 atanfunction 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 atanfunction
<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 atan2function 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ø</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. That solves
the 90 degree
<br>problem, so all you need to worry about is the 180 degree problem.
If you
<br>calculate the angles in a local coordinate system, this is never a
problem.
<p>Jim Richards <jimr@UDel.Edu>
<br>
<hr WIDTH="100%">
<br>The ATAN2 function in Matlab eliminates the quadrature problem you
are refering to...
<p>Musolino, Mark <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.
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. Unit vectors i and j are associated with
the X and Y
<br>axes, respectively. Vector S is defined as pointing from the
proximal to
<br>the distal endpoint of the shank; thus q=S/magnitude(S). 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 < 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 < 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>
<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) 5064312
<br>FAX: (310) 5064785
<p>Here is the code in FORTRAN
<p> SUBROUTINE VECTP(V1,V2,V3)
<br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
<br>CCC
<br>C
<br>C SUBROUTINE VECTP COMPUTES THE CROSS PRODUCT OF VECTOR V1 CROSS
VECTOR V2.
<br>C THE ANSWER IS RETURNED IN VECTOR V3.
<br>C
<br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
<br>CCC
<br> IMPLICIT REAL*8 (AH,OZ)
<br> DIMENSION V1(3),V2(3),V3(3)
<br> V3(1)=(V1(2)*V2(3))(V1(3)*V2(2))
<br> V3(2)=(V1(3)*V2(1))(V1(1)*V2(3))
<br> V3(3)=(V1(1)*V2(2))(V1(2)*V2(1))
<br> RETURN
<br> END
<p> SUBROUTINE SCALARP(Z1,Z2,DOTP)
<br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
<br>CCC
<br>C
<br>C SUBROUTINE SCALARP COMPUTES THE SCALAR PRODUCT OF VECTOR Z1
AND VECTOR
<br>Z2.
<br>C THE ANSWER IS RETURNED IN DOTP.
<br>C
<br>CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
<br>CCC
<br> REAL*8 DOTP,Z1(3),Z2(3)
<br> INTEGER*4 K
<br> DOTP=0.0
<br> DO 1 K=1,3
<br> DOTP=DOTP+(Z1(K)*Z2(K))
<br>1 CONTINUE
<br> RETURN
<br> 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
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>
<p>Richard Baker <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 THEN angle = ATN(y/x) ELSE angle
= SGN(y) * pi / 2
<br>IF x < 0 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 606127251, U.S.A.
<br>Phone: Office 3129961504, Lab 3124131137
<br>
Fax 3129964583
<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>
In the third quadrant, both x and y are negative, but because their ratio
<br>
is positive, ATN returns a positive angle (between 0 and pi/2). When pi
is
<br>
added to it (based on x being negative), the result is between pi and
<br>
3pi/2, which is indeed in the third quadrant.
<p>
In the fourth quadrant, x is positive and y is negative. ATN returns a
<br>
negative angle (between 0 and pi/2), which is correct.
<p>
For completeness' sake I might add that in the second quadrant x is
<br>
negative and y is positive. Thus ATN returns a negative angle (between
0
<br>
and pi/2). And when pi is added to it (based on x being negative), the
<br>
result is between pi and +pi/2, which is indeed in the second quadrant.
<p>
Thus, the procedure returns the correct value in all quadrants.
<p>
However, the values returned vary between 0 and 3pi/2 for quadrants 13,
<br>
and between 0 and pi/2 for quadrant 4. There is nothing wrong with this.
<br>
But if one prefers 0 to +  pi ranges, one can add a step in which, if
the
<br>
angle is greater than pi, then 2pi is subtracted from it. Or if one prefers
<br>
the range from 0 to 2pi, one can simply ask that if the angle is negative,
<br>
then 2pi should be added to it.
<p>
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 trigfunctions? 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<0 then answer=answer+2pi"
<p>I assume that using machinecoded 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> Zcomponent of i x q =
<br> 1) sin(theta)
(0 <= theta < pi )
<br> 2) sin(2pi  theta) = sin(theta)
= sin(theta) (pi <= theta
< 2pi)
<p>Thus
<p> Zcompone
<br>nt of i x q = sin(theta)
(0 <= theta < 2pi)
<p>Likewise, use of the dot product provides the cosine value:
<p> cos(theta) = ( i @ q ) / iq = 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> qx = cos(theta)
<br> 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> float x;
<br> float y;
<br>}
<p>There are three functions that are related to the angle computation:
<p> float C2DPoint::AngPos();
<br> void C2DPoint::RotateBy( C2DPoint vc );
<br> C2DPoint C2DPoint::UnitVector();
<p>For example, you have two vectors A and B on the XYplane. If you want
to compute the angular
<br>position of vector A to the Xaxis, it can be done as follows:
<p> C2DPoint A;
<br> // assign values to A here
<br> A =
<br>A.UnitVector(); // make A a unit vector
<br> float angle = A.AngPos(); // '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> C2DPoint A, B;
<br> // Assign values to A and B here
<br> A.RotateBy( B ); // This function
first compputes the unit vector of B and rotate vector
<br>A using B.
<br>
// Thus, the Xaxis rotates to B.
<br> A = A.UnitVector();
<br> float angle = A.AngPos();
<p>The AngPos() function is similar to ATAN2 b
<br>ut has additional features:
<p>float C2DPoint::AngPos()
<br>{
<br> // x = cos, y = sin
<br> if ( x = 99999.f  y ==
99999.f ) // if the vector is not initialized yet
<br>
return 99999.f;
<p> if ( x == 0.f &&
y == 0.f ) // if both components are zeros
<br>
return 99999.f;
<p> if ( x > 0.f ) // Quad
I & IV
<br> {
<br>
return (float)atan( y / x );
<br> }
<br> else if ( x < 0.f )
<br> {
<br>
if ( y > 0.f ) // Quad II: move from Quad IV to II by adding pi
<br>
return (float)atan( y / x ) + 3.141592f;
<br>
else
// Quad III: move from Quad
<br> I to Quad III by subtracting pi
<br>
return (float)atan( y / x )  3.141592f;
<br> }
<br> else //
0.0
<br> {
<br>
if ( y > 0.f ) // pi/2
<br>
return 3.141592f / 2.f;
<br>
else
// pi/2
<br>
return 3.141592f / 2.f;
<br> }
<br>}
<p>Note here that the ATAN function returns angle in the range of pi/2
< theta < 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
< 0). If the cosine is 0
<br>and the sine is not,
<br> the angle is either pi/2 (sine > 0) or pi/2 (sine < 0). The
AngPos() function returns a value
<br>in the range of pi <= theta < pi.
<p>The RotateBy function is simple:
<p>void C2DPoint::RotateBy( C2DPoint &vc )
<br>{
<br> if ( x == 99999.f  y
== 99999.f ) // the current vector is not initialized
<br>
return; // quit
<p> if ( vc.x == 99999.f 
vc.y == 99999.f ) // vc is not initialized
<br>
return; // no rotation
<br> else
<br>
vc = vc.UnitVector(); // make vc a unit vector
<p> // rotation
<br> C2DPoint tp = *this;
/
<br>/ copy the current vector to a temporary vector
<br> // x y
x = cos, y = sin
<br> // y x
<br> x = vc.x * tp.x + vc.y *
tp.y; // rotation
<br> 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> float mag = Magnitude();
// this function computes magnitude of the vector
<br> C2DPoint rtn;
<br> if ( mag != 0.f )
<br> {
<br>
rtn.x = x / mag;
<br>
rtn.y = y / mag;
<br> }
<br> 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 & vc )
<br>{
<br> // check if the vector is
missing
<br> if ( x == 99999.f  y
== 99999.f  vc.x == 99999.f  vc.y == 99999.f )
<br>
return 99999.f;
<br>
<p> if ( ( x == 0.f &&
y == 0.f )  ( vc.x == 0.f && vc.y == 0.f ) )
<br>
return 99999.f;
<p> // compute the unit vectors
<br> C2DPoint A = UnitVector();
<br> C2DPoint B = vc.UnitVector();
<p> // compute the angle using
the acos function
<br> 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> C2DPoint A, B; // A:
thigh vector (hip>knee), B: shank vector (knee>ankle)
<br> // assign values to the vectors here
<br> A = A * (1
<br>.f); // reverse the direction of A (knee>hip)
<br> 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> C2DPoint A, B; // A:
thigh vector (hip>knee), B: shank vector (knee>ankle)
<br> // assign values to the vectors here
<br> A.RotateBy( B );
<br> float angpos = A.AngPos();
// 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 3D 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 2D 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 > UserAngle Issues
<br>Theories > Volume II > Joint Kinematics > Computation of the Orientation
Angles
<p>Merry Christmas!
<p>YoungHoo Kwon
<br>
<br> YoungHoo 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 762045647
<br> Phone: (940) 8982598
<br> Fax: (940) 8982581
<br> Email: ykwon@twu.edu <<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 < 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>
<p>Sincerely,
<br>ZongMing 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... It
<br>seems that 2D analysis is far from obsolete.
<p>Anyway, maybe not all programming languages use the
<br>same order of arguments in the atan2 function. In Matlab,
<br>Fortran, and C, the correct syntax would be
<p> 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<=theta1<=2Pi and Pi<=theta2<=+Pi
<br>For rotations relative to the positive 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 < 0.0 ){ angle = 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 < 0 ){ angle = angle + pi; }
<br>IF( angle < 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 < 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 < 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<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. It seems that everything has been
said and
<br>we should try to end this specific topic. Below is a brief summary
and
<br>some comments on analogous issues in 3D analysis. Please feel
free to
<br>point out any errors, and further contributions on 3D joint angles
are
<br>certainly welcome. Who knows, we may even revive the "helical
angle
<br>wars" that took place on BiomchL 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. This
method
<br>uses the fact that a dot product is proportional to the cosine of the
<br>angle between two vectors. 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 & minus 180 degree range of motion
and
<br>has the advantage that it also works in 3D.
<p>As said by Paolo, and in earlier contributions, ATAN2 can be used
<br>to obtain the absolute orientation of a body segment. It is then
<br>a simple matter to subtract the orientations of two segments to get
<br>the angle between them. This will give the same result as the
"dot
<br>product" method, but can only be used with 2D 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)! I think the y,x order came about because atan2(a,b)
was
<br>thought of as a 360degree version of atan(a/b). I verified the
<br>parameter order for Matlab, Fortran, and C. 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.
As a general
<br>rule, "if" is risky when applied to measured data; an "if" can amplify
noise
<br>when not used carefully. 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 realworld data. 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. I have made such mistakes myself...
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 3D (cardanic) joint angles. 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>YoungHoo 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 3D joint angles (theta, phi, psi) as follows:
<p> cos(th)cos(ps) sin(ph)sin(th)cos(ps)+cos(ph)sin(ps)
<br>cos(ph)sin(th)cos(ps)+sin(ph)sin(ps)
<p>T = cos(th)sin(ps) sin(ph)sin(th)sin(ps)+cos(ph)sin(ps)
<br>cos(ph)sin(th)sin(ps)+sin(ph)cos(ps)
<p> sin(th)
sin(ph)cos(th)
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> ph = atan2(t32,t33)
<p> th = atan2(t31,sqrt(t32*t32+t33*t33))
<p> 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 BiomchL in 1990, consistently uses ATAN2 whereas some
more
<br>recent authors (including Kwon3D and KineMat) do not. 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://biomchl.isbweb.org">http://biomchl.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 BiomchL readership has rediscovered
<br>this wonderful mathematical function ATAN2. 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 BiomchL "helical angle wars" can be found at
<br><a target="_blank" href="http://biomchl.isbweb.org">http://biomchl.isbweb.org</a>
<p>
<p>Ton van den Bogert, BiomchL comoderator
<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 (ND20)
<br>Cleveland, OH 44195, USA
<br>Phone/Fax: (216) 4445566/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>> ph = atan2(t32,t33)
<br>> th = atan2(t31,sqrt(t32*t32+t33*t33))
<br>> 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>
<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>YoungHoo Kwon
<br>
<br> YoungHoo 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. 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. Hence, these
<br>angles are a correct description of the motion that was measured.
<br>The problem that YoungHoo 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). This 180 degree range for the second rotation
<br>is a common convention for cardanic angles, this is needed to obtain
unique
<br>angles. 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 YoungHoo.
<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). 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> You will see a similar case in the ring event of gymnastics.
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. 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. 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>YoungHoo Kwon
<br><a href="mailto:bogert@bme.ri.ccf.org"></a>
</body>
</html>