squat

Statistics for studying 3D object orientations over time

L. Bellanger

Department of Mathematics Jean Leray, UMR CNRS 6629, Nantes University, Ecole Centrale de Nantes, France

A. Stamm

Department of Mathematics Jean Leray, UMR CNRS 6629, Nantes University, Ecole Centrale de Nantes, France

2024-07-11

Context

The eGait project

Two key numbers in France

Multiple Sclerosis

100,000

Parkinson's disease

160,000

Two key observations

  • Gait impairment: major symptom impacting quality of daily life;
  • Clinical gait assessment mostly qualitative and biased:
    • Only in a constrained environment (as opposed to free-living environment);
    • Mainly based on the expertise of the neurologist;
    • Quantitative measures boil down to timing a given walking distance.

The eGait device

  • An inertial measurement unit (IMU): defines the data (rotations over time),
  • A smartphone application: collects the data,
  • Statistical methods for rotation-valued functional data: analyses the data
(a) Data recording
(b) Smartphone & App
(c) Rotation data
Figure 1: The eGait device (Brard et al. 2022; Drouin 2022; Drouin et al. 2023; Ballante et al. 2023).

Rotation data

Rotation matrix.

Rotation matrix.

Roll-pitch-yaw angles. Courtesy of Wikimedia Commons. CC-BY-SA 4.0.

Roll-pitch-yaw angles. Courtesy of Wikimedia Commons. CC-BY-SA 4.0.

Axis-angle representation.

Axis-angle representation.

Rotation data

Unit quaternions

Unit quaternions are hypercomplex numbers (Hamilton) of unit norm that encode a 3D rotation via:

\[ \mathbf{q} = (q_w, q_x, q_y, q_z) = \left(\cos \frac{\theta}{2}, i v_x \sin \frac{\theta}{2}, j v_y \sin \frac{\theta}{2}, k v_z \sin \frac{\theta}{2} \right), \]

where:

  • \(\mathbf{v} = (v_x, v_y, v_z)^\top\) and \(\theta\) are resp. axis and angle of rotation,
  • \(i\), \(j\) and \(k\) are s.t. \(i^2 = j^2 = k^2 = ijk = −1\) and \(\| q \| = 1\).

Mathematics of quaternions

Algebra of unit quaternion

Usual addition and scalar multiplication do not apply.

\[ \| q_1 + q_2 \| \ne 1 \quad \| \lambda q \| \ne 1 \]

Multiplication follows Hamilton’s rule:

\[ \mathbf{p} \mathbf{q} = \begin{pmatrix} p_w q_w + p_x q_x + p_y q_y + p_z q_z \\ p_w q_x - p_x q_w - p_y q_z + p_z q_y \\ p_w q_y + p_x q_z - p_y q_w - p_z q_x \\ p_w q_z - p_x q_y + p_y q_x - p_z q_w \end{pmatrix} \]

Lie group of unit quaternions

Sola, Deray, and Atchuthan (2018)

The squat package

Using the Eigen library

Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.

Quaternion class

  • Conversion to other matrix representations;
  • Geodesic distance;
  • Conjugate / inverse;
  • Hamilton product;
  • Norm computation and normalization;
  • Dot product;
  • Spherical linear interpolation (SLERP);
  • Random quaternion generation.

How to use Eigen in R?

RcppEigen

QTS class for quaternion time series

The qts class

A tibble with 5 columns:

  • time: Time points at which quaternions were collected;
  • w : 1st coord. of each quaternion;
  • x : 2nd coord. of each quaternion;
  • y : 3rd coord. of each quaternion;
  • z : 4th coord. of each quaternion.

An example

time w x y z
0 0.9942745 0.0797319 0.0698791 0.0133386
1 0.9948330 0.0745706 0.0676322 0.0131288
2 0.9954224 0.0693108 0.0645704 0.0126874
3 0.9960205 0.0639757 0.0609104 0.0118434
4 0.9965184 0.0594915 0.0574216 0.0107045
5 0.9969400 0.0557187 0.0540293 0.0093213
6 0.9972854 0.0527473 0.0507937 0.0077156

Methods for QTS objects

List of methods

RcppEigen for interfacing with Eigen:

  • centring();
  • ggplot2::autoplot();
  • graphics::plot();
  • +, -, *, inverse_qts().

An example

vespa64$igp[[1]] * vespa64$igp[[2]]
# A tibble: 101 × 5
    time         w         x         y         z
   <int> <dec:.5!> <dec:.5!> <dec:.5!> <dec:.5!>
 1     0   0.97621   0.16129   0.14303   0.02323
 2     1   0.97856   0.15079   0.13832   0.02344
 3     2   0.98096   0.14017   0.13241   0.02327
 4     3   0.98339   0.12929   0.12539   0.02235
 5     4   0.98557   0.11930   0.11828   0.02086
 6     5   0.98744   0.11058   0.11124   0.01875
 7     6   0.98894   0.10387   0.10465   0.01602
 8     7   0.99017   0.09883   0.09812   0.01261
 9     8   0.99126   0.09487   0.09128   0.00851
10     9   0.99225   0.09158   0.08394   0.00386
# ℹ 91 more rows

QTS sample class

The qts_sample class is a list of qts objects.

List of methods

  • as_qts_sample(), is_qts_sample(), [;
  • append();
  • rnorm_qts();
  • scale();
  • mean(), median();
  • autoplot(), plot().

Manipulating (sets of) QTS

A number of methods which can be applied to both the qts and qts_sample classes.

List of selected methods

  • log() and exp(): maps from and to the Lie algebra;
  • normalize(): ensures that all quaternions in the time series are unit quaternions;
  • resample(): performs uniform resampling using SLERP;
  • smooth(): performs smoothing of a QTS by SLERP interpolation;
  • moving_average(): performs QTS smoothing via moving average;
  • hemispherize(): ensures that there are no discontinuities in QTS due to quaternion flips since two unit quaternions \(\mathbf{q}\) and \(-\mathbf{q}\) encode the same rotation.

An example of analysis with the squat package

Individual Gait Patterns

Cheung et al. (2017)

Cheung et al. (2017)

From IGP space to \(L^2\) space

Original manifold \(\mathbb{S}^3\)

\[ \scriptsize{ \begin{array}{rccc} \mathbf{q}: & [0,1] & \to & \mathbb{S}^3 \\ & s & \mapsto & \mathbf{q}(s) \end{array} } \]

Tangent space \(\mathcal{T}\mathbb{S}^3 \approx \mathbb{R}^3\)

\[ \scriptsize{ \begin{array}{rccc} \mathbf{t}: & [0,1] & \to & \mathbb{R}^3 \\ & s & \mapsto & \log(\mathbf{q}(s)) = (\theta(s) / 2) \mathbf{v}(s) \end{array} } \]

Square-root velocity space \(L^2 \left( [0, 1], \mathbb{R}^3 \right)\) (Srivastava and Klassen 2016)

\[ \scriptsize{ \begin{array}{rccc} \mathbf{g}: & [0,1] & \to & \mathbb{R}^3 \\ & s & \mapsto & \begin{cases} \frac{\mathbf{t}^\prime(s))}{\sqrt{\| \mathbf{t}^\prime(s)) \|}} & \text{if } \mathbf{t}^\prime(s) \neq 0 \\ 0 & \text{otherwise} \end{cases} \end{array} } \]

Statistical analyses

Clustering - Deps: fdacluster, fdasrvf, dbscan (Sangalli et al. 2010; Vantini 2012)

  • S3 impl. of kmeans(), hclust() and dbscan() for qts_sample objects;
  • Return an object of class qtsclust;
  • S3 impl. of autoplot() and plot() for qtsclust objects for visualization.

PCA - Deps: MFPCA (Happ and Greven 2018; Happ-Kurz 2020)

  • S3 impl. of prcomp() for qts_sample objects;
  • Returns an object of class prcomp_qts;
  • S3 impl. of autoplot(), plot() and predict() for prcomp_qts objects for easy visualizations and prediction.

Clustering Example

km <- kmeans(
  x = vespa64$igp, # object of class qts_sample
  n_clusters = 2, # number of clusters to find
  seeding_strategy = "kmeans++", # various seeding strategies
  is_domain_interval = TRUE, # are curves defined on a fixed common interval?
  transformation = "srsf", # transformation to apply to the curves
  warping_class = "bpd", # warping class to use
  metric = "l2" # distance metric to use
)
plot(km)
plot(km$best_clustering, type = "amplitude") # S3 specialization for `caps` objects
plot(km$best_clustering, type = "phase") # S3 specialization for `caps` objects

PCA Example

pca <- prcomp(
  x = vespa64$igp, # object of class qts_sample
  M = 5L, # number of principal components to keep
  fit = FALSE # should a reconstruction of the sample from the retained PCs be stored?
)
screeplot(pca)
plot(pca, what = "PC1")
autoplot(pca, what = "scores") + geom_point(aes(color = vespa64$V))

Wrappin’ up

The squat package currently allows you to analyze time series of 3D object orientations

  • seamlessly for clustering and PCA;
  • with a bit of extra work for other analyses (e.g. regression/prediction) via the provided exp() and log() maps.

Do not hesitate to create an issue on the GitHub repository if you want other statistical analyses to be implemented in the same way as the clustering and PCA.

References

userR! 2024 - Salzburg -

Ballante, Elena, Lise Bellanger, Pierre Drouin, Silvia Figini, and Aymeric Stamm. 2023. “Smoothing Method for Unit Quaternion Time Series in a Classification Problem: An Application to Motion Data.” Scientific Reports 13 (1): 9366.
Brard, Raphaël, Lise Bellanger, Laurent Chevreuil, Fanny Doistau, Pierre Drouin, and Aymeric Stamm. 2022. “A Novel Walking Activity Recognition Model for Rotation Time Series Collected by a Wearable Sensor in a Free-Living Environment.” Sensors 22 (9): 3555.
Cheung, Ada S, Hans Gray, Anthony G Schache, Rudolf Hoermann, Daryl Lim Joon, Jeffrey D Zajac, Marcus G Pandy, and Mathis Grossmann. 2017. “Androgen Deprivation Causes Selective Deficits in the Biomechanical Leg Muscle Function of Men During Walking: A Prospective Case–Control Study.” Journal of Cachexia, Sarcopenia and Muscle 8 (1): 102–12.
Drouin, Pierre. 2022. “Amélioration Du Suivi Des Patients Atteints de Maladies Neuro-dégénératives à l’aide d’objets Connectés.” PhD thesis, Nantes Université.
Drouin, Pierre, Aymeric Stamm, Laurent Chevreuil, Vincent Graillot, Laetitia Barbin, Pierre-Antoine Gourraud, David-Axel Laplaud, and Lise Bellanger. 2023. “Semi-Supervised Clustering of Quaternion Time Series: Application to Gait Analysis in Multiple Sclerosis Using Motion Sensor Data.” Statistics in Medicine 42 (4): 433–56.
Happ, Clara, and Sonja Greven. 2018. “Multivariate Functional Principal Component Analysis for Data Observed on Different (Dimensional) Domains.” Journal of the American Statistical Association 113 (522): 649–59.
Happ-Kurz, Clara. 2020. “Object-Oriented Software for Functional Data.” Journal of Statistical Software 93 (5): 1–38. https://doi.org/10.18637/jss.v093.i05.
Kurtek, Sebastian, Anuj Srivastava, Eric Klassen, and Zhaohua Ding. 2012. “Statistical Modeling of Curves Using Shapes and Related Features.” Journal of the American Statistical Association 107 (499): 1152–65.
Sangalli, Laura M, Piercesare Secchi, Simone Vantini, and Valeria Vitelli. 2010. “K-Mean Alignment for Curve Clustering.” Computational Statistics & Data Analysis 54 (5): 1219–33.
Sola, Joan, Jeremie Deray, and Dinesh Atchuthan. 2018. “A Micro Lie Theory for State Estimation in Robotics.” arXiv Preprint arXiv:1812.01537.
Srivastava, Anuj, and Eric P Klassen. 2016. Functional and Shape Data Analysis. Vol. 1. Springer.
Vantini, Simone. 2012. “On the Definition of Phase and Amplitude Variability in Functional Data Analysis.” Test 21 (4): 676–96.