Marcos Duarte
Laboratory of Biomechanics and Motor Control (http://pesquisa.ufabc.edu.br/bmclab)
Federal University of ABC, Brazil# Angular kinematics in a plane (2D)
Human motion is a combination of linear and angular movement and occurs in the three-dimensional (3D) space. For certain movements and depending on the desired or needed degree of detail for the motion analysis, it's possible to perform a two-dimensional (2D, planar) analysis at the main plane of movement. Such simplification is appreciated because the instrumentation and analysis are much more complicated in order to measure 3D motion than for the 2D case.
# Import the necessary libraries
import numpy as np
from IPython.display import display, Latex
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set_context("notebook", font_scale=1.2, rc={"lines.linewidth": 2, "lines.markersize": 8})
For the planar case, the calculation of angles is reduced to the application of trigonometry to the kinematic data. For instance, given the coordinates in a plane of markers on a segment as shown in the figure below, the angle of the segment can be calculated using the inverse function of $\sin,\cos$ , or $\tan$ .
For better numerical accuracy (and also to distinguish values in the whole quadrant), the inverse function of $\tan$ is preferred.
For the data shown in the previous figure:
In computer programming (here, Python/Numpy) this is calculated using: numpy.arctan((y2-y1)/(x2-x1))
. However, for the previous figure the function arctan
can not distinguish if the segment is at 45o or at 225o, arctan
will return the same value. Because this, the function numpy.arctan2(y, x)
is used, but be aware that arctan2
will return angles between $[-\pi,\pi]$:
x1, y1 = 0, 0
x2, y2 = 1, 1
display(Latex(r'Segment\;at\;45^o:'))
angs = [ np.arctan((1-0)/(1-0))*180/np.pi, np.arctan2(1-0, 1-0)*180/np.pi ]
display(Latex(r'Using\;arctan: '+str(angs[0])+r'\;Using\;arctan2: '+str(angs[1])))
display(Latex(r'Segment\;at\;225^o:'))
angs = [ np.arctan((-1-0)/(-1-0))*180/np.pi, np.arctan2(-1-0, -1-0)*180/np.pi ]
display(Latex(r'Using\;arctan: '+str(angs[0])+r'\;Using\;arctan2: '+str(angs[1])))
And Numpy has a function to convert an angle in rad to degrees and the other way around:
print('np.rad2deg(np.pi/2) =', np.rad2deg(np.pi))
print('np.deg2rad(180) =', np.deg2rad(180))
np.rad2deg(np.pi/2) = 180.0 np.deg2rad(180) = 3.141592653589793
Let's simulate a 2D motion of the arm performing two complete turns around the shoulder to exemplify the use of arctan2
:
t = np.arange(0, 2, 0.01)
x = np.cos(2*np.pi*t)
y = np.sin(2*np.pi*t)
ang = np.arctan2(y, x)*180/np.pi
plt.figure(figsize=(12, 4))
hax1 = plt.subplot2grid((2, 3), (0, 0), rowspan=2)
hax1.plot(x, y, 'go')
hax1.plot(0, 0, 'go')
hax1.set_xlabel('x [m]')
hax1.set_ylabel('y [m]')
hax1.set_xlim([-1.1, 1.1])
hax1.set_ylim([-1.1, 1.1])
hax2 = plt.subplot2grid((2, 3), (0, 1), colspan=2)
hax2.plot(t, x, 'bo', label='x')
hax2.plot(t, y, 'ro', label='y')
hax2.legend(numpoints=1, frameon=True, framealpha=.8)
hax2.set_ylabel('Position [m]')
hax2.set_ylim([-1.1, 1.1])
hax3 = plt.subplot2grid((2, 3), (1, 1), colspan=2)
hax3.plot(t, ang, 'go')
hax3.set_yticks(np.arange(-180, 181, 90))
hax3.set_xlabel('Time [s]')
hax3.set_ylabel('Angle [ o]')
plt.tight_layout()
Because the output of the arctan2
is bounded to $[-\pi,\pi]$, the angle measured appears chopped in the figure. This problem can be solved using the function numpy.unwrap
, which detects sudden jumps in the angle and corrects that:
ang = np.unwrap(np.arctan2(y, x))*180/np.pi
hfig, hax = plt.subplots(1,1, figsize=(8,3))
hax.plot(t, ang, 'go')
hax.set_yticks(np.arange(start=0, stop=721, step=90))
hax.set_xlabel('Time [s]')
hax.set_ylabel('Angle [ o]')
plt.tight_layout()
If now we want to measure the angle of a joint (i.e., the angle of a segment in relation to other segment) we just have to subtract the two segment angles (but this is correct only if the angles are at the same plane):
x1, y1 = 0.0, 0.0
x2, y2 = 1.0, 1.0
x3, y3 = 1.1, 1.0
x4, y4 = 2.1, 0.0
hfig, hax = plt.subplots(1,1, figsize=(8,3))
hax.plot((x1,x2), (y1,y2), 'b-', (x1,x2), (y1,y2), 'ro', linewidth=3, markersize=12)
hax.add_patch(matplotlib.patches.FancyArrowPatch(posA=(x1+np.sqrt(2)/3, y1),
posB=(x2/3, y2/3),\
arrowstyle='->,head_length=10,head_width=5', connectionstyle='arc3,rad=0.3'))
plt.text(1/2, 1/5, '$\\theta_1$', fontsize=24)
hax.plot((x3,x4), (y3,y4), 'b-', (x3,x4), (y3,y4), 'ro', linewidth=3, markersize=12)
hax.add_patch(matplotlib.patches.FancyArrowPatch(posA=(x4+np.sqrt(2)/3, y4),
posB=(x4-1/3, y4+1/3),\
arrowstyle='->,head_length=10,head_width=5', connectionstyle='arc3,rad=0.3'))
hax.xaxis.set_ticks((x1,x2,x3,x4))
hax.yaxis.set_ticks((y1,y2,y3,y4))
hax.xaxis.set_ticklabels(('x1','x2','x3','x4'), fontsize=20)
#hax.yaxis.set_ticklabels(('$y_1,\,y_4′,′y_2,\,y_3$'), fontsize=20)
plt.text(x4+.2,y4+.3,'$\\theta_2$', fontsize=24)
hax.add_patch(matplotlib.patches.FancyArrowPatch(posA=(x2-1/3, y2-1/3),
posB=(x3+1/3, y3-1/3),\
arrowstyle='->,head_length=10,head_width=5', connectionstyle='arc3,rad=0.3'))
plt.text(x1+.8,y1+.35,'$\\theta_J=\\theta_2-\\theta_1$', fontsize=24)
hax.set_xlim(min([x1,x2,x3,x4])-0.1, max([x1,x2,x3,x4])+0.5)
hax.set_ylim(min([y1,y2,y3,y4])-0.1, max([y1,y2,y3,y4])+0.1)
hax.grid(xdata=(0,1), ydata=(0,1))
plt.tight_layout()
The joint angle shown above is simply the difference between the adjacent segment angles:
x1, y1, x2, y2 = 0, 0, 1, 1
x3, y3, x4, y4 = 1.1, 1, 2.1, 0
ang1 = np.arctan2(y2-y1, x2-x1)*180/np.pi
ang2 = np.arctan2(y3-y4, x3-x4)*180/np.pi
#print('Angle 1:', ang1, '\nAngle 2:', ang2, '\nJoint angle:', ang2-ang1)
display(Latex(r'\theta_1=\;' + str(ang1) + '^o'))
display(Latex(r'\theta_2=\;' + str(ang2) + '^o'))
display(Latex(r'\theta_J=\;' + str(ang2-ang1)+'^o'))
The following convention is commonly used to describe the knee and ankle joint angles at the sagittal plane (figure from Winter 2005):
In certain cases, we have access to the 3D coordinates of markers but we just care for the angle between segments in the plane defined by these segments (but if there is considerable movement in different planes, this simple 2D angle might give unexpected results).
Consider that p1
and p2
are the 3D coordinates of markers placed on segment 1 and p2
and p3
are the 3D coordinates of the markers on segment 2.
To determine the 2D angle between the segments, one can use the definition of the dot product:
$$ \mathbf{a} \cdot \mathbf{b} = ||\mathbf{a}||\:||\mathbf{b}||\:cos(\theta)\;\;\; \Rightarrow \;\;\; angle = arccos\left(\frac{dot(p2-p1,\;p3-p2)}{norm(p2-p1)*norm(p3-p2)\;\;\;\;\;} \right) $$Or using the definition of the cross product:
$$ \mathbf{a} \times \mathbf{b} = ||\mathbf{a}||\:||\mathbf{b}||\:sin(\theta) \;\; \Rightarrow \;\; angle = arcsin\left(\frac{cross(p2-p1,\;p3-p2)}{norm(p2-p1)*norm(p3-p2)\;\;\;\;\;} \right) $$But because arctan2
has a better numerical accuracy, combine the dot and cross products, and in Python notation:
angle = np.arctan2(np.linalg.norm(np.cross(p1-p2, p4-p3)), np.dot(p1-p2, p4-p3))
See this notebook for a review on the mathematical functions cross product and scalar product.
We can use the formula above for the angle between two 3D vectors to calculate the joint angle even with the 2D vectors we calculated before:
p1, p2 = np.array([0, 0]), np.array([1, 1]) # segment 1
p3, p4 = np.array([1.1, 1]), np.array([2.1, 0]) # segment 2
angle = np.arctan2(np.linalg.norm(np.cross(p1-p2, p4-p3)), np.dot(p1-p2, p4-p3))*180/np.pi
print('Joint angle:', '{0:.1f}'.format(angle))
Joint angle: 90.0
As expected, the same result.
In Numpy, if the third components of vectors are zero, we don't even need to type them; Numpy takes care of adding zero as the third component for the cross product.
The angular position is a vector, its direction is given by the perpendicular axis to the plane where the angular position is described, and the motion if it occurs it's said to occur around this axis.
Angular velocity is the rate (with respect to time) of change of the angular position:
$$ \mathbf{\omega}(t) = \frac{\mathbf{\theta}(t_2)-\mathbf{\theta}(t_1)}{t_2-t_1} = \frac{\Delta \mathbf{\theta}}{\Delta t}$$ $$ \mathbf{\omega}(t) = \frac{d\mathbf{\theta}(t)}{dt} $$Angular acceleration is the rate (with respect to time) of change of the angular velocity, which can also be given by the second-order rate of change of the angular position:
$$ \mathbf{\alpha}(t) = \frac{\mathbf{\omega}(t_2)-\mathbf{\omega}(t_1)}{t_2-t_1} = \frac{\Delta \mathbf{\omega}}{\Delta t}$$Likewise, angular acceleration is the first-order derivative of the angular velocity or the second-order derivative of the angular position vector:
$$ \mathbf{\alpha}(t) = \frac{d\mathbf{\omega}(t)}{dt} = \frac{d^2\mathbf{\theta}(t)}{dt^2} $$The direction of the angular velocity and acceleration vectors is the same as the angular position (perpendicular to the plane of rotation) and the sense is given by the right-hand rule.
As the angular acceleration is the derivative of the angular velocity which is the derivative of angular position, the inverse mathematical operation is the antiderivative (or integral):
$$ \mathbf{\theta}(t) = \mathbf{\theta}_0 + \int \mathbf{\omega}(t) dt $$ $$ \mathbf{\omega}(t) = \mathbf{\omega}_0 + \int \mathbf{\alpha}(t) dt $$Consider a particle rotating around a point at a fixed distance r
(circular motion), as the particle moves along the circle, it travels an arc of length s
.
The angular position of the particle is:
Which is in fact similar to the definition of the angular measure radian:
Then, the distance travelled by the particle is the arc length:
$$ s = r\theta $$As the radius is constant, the relation between linear and angular velocity and acceleration is straightfoward:
$$ v = \frac{ds}{dt} = r\frac{d\theta}{dt} = r\omega $$ $$ a = \frac{dv}{dt} = r\frac{d\omega}{dt} = r\alpha $$a) Calculate the angular and lineat velocity of the gymnast's center of mass at the point of release.
b) Calculate the horizontal distance travelled by the gymnast's center of mass.
a. Calculate and plot the angles of the foot, leg, and thigh segments.
b. Calculate and plot the angles of the ankle, knee, and hip joint.
c. Calculate and plot the velocities and accelerations for the angles calculated in B.
d. Compare the ankle angle using the two different conventions described by Winter (2009), that is, defining the foot segment with the MT5 or the TOE marker.
e. Knowing that a stride period corresponds to the data between frames 1 and 70 (two subsequents toe-off by the right foot), can you suggest a possible candidate for automatic determination of a stride? Hint: look at the vertical displacement and acceleration of the heel marker.
Clik here for the data from Table A.1 (Winter, 2009) from Winter's book student site.
Example: load data and plot the markers' positions:
# load data file
# use Pandas just to read data from internet
data = pd.read_csv('https://raw.githubusercontent.com/BMClab/BMC/master/data/WinterTableA1.txt',
sep=' ', header=None, skiprows=2).to_numpy()
markers = ['RIB CAGE', 'HIP', 'KNEE', 'FIBULA', 'ANKLE', 'HEEL', 'MT5', 'TOE']
fig = plt.figure(figsize=(10, 6))
ax = plt.subplot2grid((2,2),(0, 0))
ax.plot(data[: ,1], data[:, 2::2])
ax.set_xlabel('Time [s]', fontsize=14)
ax.set_ylabel('Horizontal [cm]', fontsize=14)
ax = plt.subplot2grid((2, 2),(0, 1))
ax.plot(data[: ,1], data[:, 3::2])
ax.set_xlabel('Time [s]', fontsize=14)
ax.set_ylabel('Vertical [cm]', fontsize=14)
ax = plt.subplot2grid((2, 2), (1, 0), colspan=2)
ax.plot(data[:, 2::2], data[:, 3::2])
ax.set_xlabel('Horizontal [cm]', fontsize=14)
ax.set_ylabel('Vertical [cm]', fontsize=14)
plt.suptitle('Table A.1 (Winter, 2009): female, 22 yrs, 55.7 kg, 156 cm, ' \
'fast cadence (115 steps/min)', y=1.02, fontsize=14)
ax.legend(markers, loc="upper right", title='Markers')
plt.tight_layout()