Reframing Three-Dimensional Morphometrics Through Functional Data Innovations

TL;DR

This paper introduces innovative functional data analysis methods, including SRVF and arc-length parameterisation, to enhance 3D morphometric analysis, demonstrated through simulations and classification of kangaroo dietary categories.

Contribution

It develops seven new pipelines integrating functional data analysis with geometric morphometrics, improving 3D shape analysis by considering curvature and elastic alignment.

Findings

Functional data analysis improves 3D morphometric robustness

Arc-length and SRVF approaches enhance shape differentiation

New pipelines outperform standard geometric morphometrics in classification tasks

Abstract

This study innovates geometric morphometrics by incorporating functional data analysis, the square-root velocity function (SRVF), and arc-length parameterisation for 3D morphometric data, leading to the development of seven new pipelines in addition to the standard geometric morphometrics (GM) approach.. This enables three-dimensional images to be examined from perspectives that do not neglect curvature, through the combined use of arc-length parameterisation, soft-alignment, and elastic-alignment. A simulation study was conducted to demonstrate the general effectiveness of eight pipelines: geometric morphometrics (GM, baseline), arc-GM, functional data morphometrics (FDM), arc-FDM, soft-SRV-FDM, arc-soft-SRV-FDM, elastic-SRV-FDM, and arc-elastic-SRV-FDM. These pipelines were also applied to distinguish dietary categories of kangaroos (omnivores, mixed feeders, browsers, and grazers) using cranial landmarks obtained from 41 extant species. Principal component analysis was conducted, followed by classification analysis using linear discriminant analysis, multinomial regression and support vector machines with a linear kernel. The results highlight the effectiveness of functional data analysis, together with arc-length and SRVF-based approaches, in opening the door to more robust perspectives for analysing three-dimensional morphometrics, while establishing geometric morphometrics as the baseline for comparison.

Tables

Table 1. Table 1: Number of principal components required to capture 95% of cumulative variance across pipelines

Pipeline	Number of Principal Components
GM	60
FDM	14
Arc-GM	41
Arc-FDM	15
Elastic-SRV-FDM	58
Soft-SRV-FDM	13
Arc-Elastic-SRV-FDM	51
Arc-Soft-SRV-FDM	13

Table 2. Table 2: Mean squared error of raw coordinates across pipelines (mean with standard deviation in parentheses).

Pipeline	MSE of Raw Coordinates
GM	0.01994 (0.00034)
FDM	0.00016 (0.00081)
Arc-GM	0.02089 (0.00363)
Arc-FDM	0.00425 (0.00439)
Elastic-SRV-FDM	0.00261 (0.00112)
Soft-SRV-FDM	0.00204 (0.00081)
Arc-Elastic-SRV-FDM	0.00197 (0.00126)
Arc-Soft-SRV-FDM	0.00192 (0.00120)

Table 3. Table 3: Number of principal components required to capture 95% of cumulative variance across pipelines

Pipeline	Number of principal Components
GM	9
FDM	3
Arc-GM	13
Arc-FDM	5
Elastic-SRV-FDM	4
Soft-SRV-FDM	3
Arc-Elastic-SRV-FDM	15
Arc-Soft-SRV-FDM	4

Table 4. Table 4: Mean squared error of kangaroo cranial coordinates across pipelines (mean with standard deviation in parentheses)

Pipeline	MSE of Raw Coordinates
GM	0.01581 (0.00001)
FDM	0.00299 (0.00010)
Arc-GM	0.01385 (0.00068)
Arc-FDM	0.01345 (0.00069)
Elastic-SRV-FDM	0.00424 (0.00033)
Soft-SRV-FDM	0.00481 (0.00018)
Arc-Elastic-SRV-FDM	0.00841 (0.00056)
Arc-Soft-SRV-FDM	0.00794 (0.00042)

Equations

$$\begin{gathered} {\mathbf{x}}_k=(x_{k,{\tau }_1}, \\ \dots , \\ {x}_{k,{\tau }_N}),{\mathbf{y}}_k=(y_{k,{\tau }_1}, \\ \dots , \\ {y}_{k,{\tau }_N}),{\mathbf{z}}_k=(z_{k,{\tau }_1}, \\ \dots , \\ {z}_{k,{\tau }_N}), \end{gathered}$$

$$s(t)={\int }_0^t\Vert \dot{f}(\tau )\Vert d\tau ,$$

$$g(u)=f(t(u)),\text{with}\Vert \frac{dg}{du}\Vert =1.$$

$$\{X^{(1)}(\cdot ),\dots ,X^{(n)}(\cdot )\},$$

$$\{Y^{(1)}(\cdot ),\dots ,Y^{(n)}(\cdot )\}\text{and}\{Z^{(1)}(\cdot ),\dots ,Z^{(n)}(\cdot )\},$$

$$\mathbf{X}=(X,Y,Z{)}^{\top },$$

$$X:I_x\to \mathbb{R},X\in L^2(I_x),$$

$$I:=I_x\times I_y\times I_z,$$

$$H:=L^2(I_x)\times L^2(I_y)\times L^2(I_z).$$

$$⟨⟨f,g⟩⟩:=\sum _{p\in \{x,y,z\}}⟨f_p,g_p⟩=\sum _{p\in \{x,y,z\}}{\int }_{I_p}{f}_p(t_p)g_p(t_p)d{t}_p,$$

$$\mathbb{E}[\mathbf{X}(t)]:=(\mathbb{E}[X(t_x)],\mathbb{E}[Y(t_y)],\mathbb{E}[Z(t_z)]{)}^{\top }=0,\forall t=(t_x,t_y,t_z{)}^{\top }\in I.$$

$$C(s,t)=\mathbb{E}[\mathbf{X}(s)\mathbf{X}(t{)}^{\top }],$$

$$C_{p,q}(s_p,t_q)=\mathbb{E}[X_p(s_p)X_q(t_q)]=\mathrm{Cov}(X_p(s_p),X_q(t_q)),$$

$$(\Gamma f{)}^{(q)}(t_q):=⟨⟨C_{\cdot ,q}(\cdot ,t_q),f(\cdot )⟩⟩=\sum _{p=1}^P{\int }_{I_p}{C}_{p,q}(s_p,t_q)f_p(s_p)d{s}_p.$$

$$\Gamma {\phi }_j={\lambda }_j{\phi }_j,{\lambda }_j\to 0\text{as}j\to \infty .$$

$$\mathbf{X}(t)=\sum _{j=1}^{\infty }{\xi }_j{\phi }_j(t),t\in I,$$

$${\mathbf{X}}_{[J]}(t)=\sum _{j=1}^J{\xi }_j{\phi }_j(t),$$

$$X_{p,[J_p]}(t_p)=\sum _{j=1}^{J_p}{\xi }_{p,j}{\phi }_{p,j}(t_p),$$

$${\hat{K}}_p(s,t)=\frac{1}{n-1}\sum _{i=1}^n{X}_p^{(i)}(s)X_p^{(i)}(t),$$

$$\hat{Z}=\frac{1}{n-1}{\Xi }^{\top }\Xi .$$

$${\hat{\phi }}_{p,j}(t_p)=\sum _{k=1}^{J_p}[{\hat{v}}_j{]}_{(p,k)}{\hat{\phi }}_{p,k}(t_p).$$

$$f^{(i)}:T\to {\mathbb{R}}^3,$$

$$f^{(i)}(t)=(x^{(i)}(t),y^{(i)}(t),z^{(i)}(t){)}^{\top }.$$

$$q(t)=\frac{\dot{f}(t)}{\sqrt{\Vert \dot{f}(t)\Vert }},t\in [0,1].$$

$$Q^{(i)}(t)=(q_x^{(i)}(t),q_y^{(i)}(t),q_z^{(i)}(t){)}^{\top }.$$

$$t_i=\frac{i-1}{n_{\mathrm{points}}-1},i=1,\dots ,n_{\mathrm{points}}.$$

$$x(t)$$

$$y(t)$$

$$z(t)$$

$$x(t)$$

$$y(t)$$

$$z(t)$$

Full text

a Institute of Mathematical Sciences, Faculty of Science, Universiti Malaya, 50603 Kuala Lumpur, Malaysia

b Mathematics, Division of Science, New York University Abu Dhabi, United Arab Emirates

c Laboratoire Paul Painlevé UMR CNRS 8524, INRIA–MODAL, Université de Lille, France

d Universiti Malaya Centre for Data Analytics, Universiti Malaya, 50603 Kuala Lumpur, Malaysia

e Center of Research for Statistical Modelling and Methodology, Faculty of Science, Universiti Malaya, 50603 Kuala Lumpur, Malaysia

*Corresponding author: dharini@um.edu.my

Keywords: Arc-length parameterisation; elastic alignment; functional data analysis; geometric morphometrics; square-root velocity function

1  Introduction

Capturing three-dimensional shape variation is crucial for disciplines spanning evolutionary biology to medical diagnostics. Over the years, a variety of shape parameterisation strategies have been developed, including landmark-based models, boundary point distributions (e.g., semi-landmarks and outlines), skeletal models, functional representations, and topological descriptors. Each has distinct advantages, such as boundary models that emphasise object contours, skeletal models that compactly encode internal structure, and topological approaches that provide scale-invariant insights. Nevertheless, landmark-based models remain dominant in geometric morphometrics due to their intuitive interpretation and compatibility with anatomical homology Drydmard16 ; zelditch2004geometric .

In classical statistical shape analysis, landmark configurations are typically aligned using generalised Procrustes analysis (GPA) to remove translation, rotation, and scale, followed by principal component analysis (PCA) to describe variation Drydmard98 ; webster2010practical . These Procrustean methods treat aligned configurations as elements of a Euclidean space, enabling standard multivariate statistical techniques. Such extrinsic approaches perform inference in a linearised setting rather than respecting the nonlinear geometry of the shape manifold patrangenaru2015nonparametric .

Intrinsic methods recognise that shape spaces are typically non-Euclidean and conduct inference directly on the manifold using geodesic distances, Fréchet means, and Riemannian metrics bhattacharya2012nonparametric ; patrangenaru2015nonparametric . These are especially relevant when the curvature of the shape space affects relationships, such as in closed curves or deformable 3D objects. However, reducing continuous curves or surfaces to sparse sets of landmarks inevitably discards fine-scale detail and may introduce bias when variation is large or the manifold is strongly curved kent1997consistency ; huckemann2010intrinsic . These limitations have prompted the development of functional data alignment-invariant models and methods that better adhere to the geometry of shape space.

Functional data analysis (FDA) presents then a promising alternative by modeling curves as smooth functions within a functional space ramsay2005functional , treating the shape as a realisation of a stochastic process rather than a set of discrete coordinates. Srivastava and Klassen (2016)SrivastavaKlassen2016 motivate the need for rigorous statistical tools to analyse functions and shapes while handling invariances such as translation, rotation, scaling, and re-parameterisation. Their geometric framework unifies registration and shape analysis, enabling tools such as functional principal component analysis (FPCA) for vector-valued functions. Wu et al. (2024)wu2024shape presents a framework that integrates geometric shape analysis with functional data methods to better capture structural features of curves, such as peaks and valleys, beyond what traditional based approaches offer. The framework developed tools for shape fitting, shape-fPCA, and shape regression, extending these ideas to manifold-valued functions and demonstrating their effectiveness on real and simulated datasets.

A critical consideration when modeling curve-like shapes is the choice of parameterisation. Arc-length parameterisation enables consistent assessment of complex-shaped signals, including hysteretic load-unload curves, by eliminating variability due to uneven sampling or velocity hartlen2022arclength . In shape analysis, the space of unparameterised curves is typically modeled as a quotient of parameterised curves under the action of the reparameterisation group, and arc-length parameterisation is often used as a canonical representative to identify each equivalence class, leading to uniformly sampled and geometry-preserving shape comparisons Tumpach2017 .

Functional data morphometrics (FDM) represents landmark trajectories as multivariate functions, allowing detailed deformation patterns and feature-based registration guo2022data ; ramsay1998curve . Multivariate functional principal component analysis (MFPCA) extends landmark trajectories to multi-dimensional trajectories, especially when temporal or arc-length ordering is essential happ2018multivariate ; yao2005functional . Parallel advances in elastic shape analysis, such as the square root velocity function (SRVF) framework, leverage the Fisher–Rao Riemannian metric to separate amplitude and phase variation, aligning curves to a Karcher mean template srivastava2011registration . This manifold-aware approach is part of a broader object-oriented data analysis perspective marron2021object that integrates geometric structure with elastic deformation metrics. Intrinsic methods offer a theoretically robust enhancement to Procrustean techniques, particularly in the context of high-dimensional or complex shape data. Katina et al. (2016)katina2021functional illustrated the application of functional data analysis to 3D facial surfaces by treating each coordinate as a function over a 2D grid and employing FPCA to encapsulate the primary variations in shape. Drawing inspiration from this functional paradigm, our study adopts a similar strategy for examining 3D landmark trajectories or curves, utilizing arc-length parameterisation and MFPCA to more effectively analyse anatomical shape variation. This paper introduces two traditional geometric morphometric methods alongside six functional data-based geometric morphometric pipelines for three-dimensional data:

GM: a classical geometric morphometrics (GM) approach was used to remove non-shape information by applying Geometric Procrustes Analysis (GPA) 2. 2.

Arc-GM: reparameterises each shape to uniform arc before applying GPA 3. 3.

FDM: a three-dimensional extension of Pillay et al. (2024)pillay2024functional which models each three-dimensional outline as a smooth multivariate functional data and applies multivariate functional principal component analysis on the raw parameterisation. 4. 4.

Arc-FDM: reparameterises each functional data (curve) to uniform arc length before spline-basis smoothing and analysis to ensure even sampling. 5. 5.

Soft-SRV-FDM: blends the identity mapping with an estimated SRVF warp, achieving a soft elastic alignment that balances global outline preservation with the ability to capture fine-scale shape variation. 6. 6.

Arc-Soft-SRV-FDM: reparameterises each shape to uniform arc before applying Soft-FDM-SVR 7. 7.

Elastic-SRV-FDM: applies full SRVF-based elastic alignment, estimating and applying the optimal warping function to remove phase variation, thereby isolating amplitude differences and capturing shape-related variation. 8. 8.

Arc-Elastic-SRV-FDM: reparameterises each shape to uniform arc before applying elastic-SVR-FDM

Figure 1 illustrates the flowchart of the eight pipelines that integrate functional data analysis and geometric morphometrics for three-dimensional morphometric data.

Of these approaches, the final seven represent innovations derived from existing methods found in the literature. Following the aforementioned geometric analyses, the results of each method are classified utilising three approaches: linear discriminant analysis (LDA), linear support vector machines (SVM), and multinomial logistic regression. These are trained on the extracted principal component scores from each method’s data. The pipelines previously mentioned were evaluated using both a simulation study and a real-world dataset consisting of 3D kangaroo skull landmarks butler20213d , illustrating significant biological diversity across various dietary categories.

Relative to traditional geometric morphometrics techniques, the functional data framework presents a cohesive and versatile paradigm that effectively accommodates both aligned and unaligned shape representations. It preserves significant phase variation through arc-parameterisation and SRVF alignment, while it also improves classification accuracy through incorporation with supervised learning models. This study elucidates the balance between alignment strategies and classification accuracy by providing both quantitative metrics and visual evaluations of reconstruction fidelity. Our findings are intended to offer practical guidance for the customisation of geometric morphometrics pipelines in accordance with specific research objectives, whether in the fields of taxonomy, evolutionary analysis, or medical morphology.

This work is organised as follows. Section 2 details the methods based on the GM and FDM frameworks. Section 3 presents a simulation study to compare FDM with its classical counterpart GM. Section 4 describes the application to 3D kangaroo landmark data. Section 5 provides the discussion for both the simulation study and empirical data. Section 6 concludes.

2  Methodology

Let

[TABLE]

$\quad k=1,\dots,n,$ be the symmetric shape landmark coordinates for $n$ specimens. Here, $\tau_{1},\dots,\tau_{N}$ are the observed landmark points, and $N$ is the number of landmarks on a dimensional domain (Figure 10) based on the crania of kangaroos Figure S1).

2.1  Generalised Procrustes Analysis

In this analytical pathway, we performed GPA. A curve-sliding matrix was provided to guide the sliding of semi-landmarks along tangent directions, minimising bending energy based on the thin-plate spline criterion. After GPA, averaged configurations were computed from the aligned coordinate array and stored. Visualisations, including grand mean shapes and reference-to-target comparisons, were used to assess alignment quality. Given the bilateral symmetry of kangaroo crania, a bilateral symmetry analysis was applied. A landmark pairing file specified corresponding bilateral landmarks. Symmetrised shape configurations were generated using object symmetry (not side or replicate symmetry), and the output was used in subsequent analyses. From the symmetrised GPA-aligned coordinates, shape variables were extracted by computing species-level means, corresponding to unique specimen identifiers. The symmetrised shape data were aggregated and reshaped into a 3D array for downstream shape-space analysis and visualisation. In parallel, centroid size, representing overall cranial scale, was averaged across replicates for each specimen and subsequently across species. This provided a size vector aligned with the shape configurations.

Generalised Procrustes analysis was conducted using the gpagen() function from the geomorph package in R adams2013geomorph , which optimally aligns configurations by minimizing the summed squared distances between corresponding landmarks. To improve alignment fidelity for curves and complex outlines, semi-landmarks were allowed to slide along tangents to the curve using thin-plate spline (TPS) bending energy minimisation, a standard method for minimising landmark placement artifacts while preserving biological signal gunz2005semilandmarks ; guo2022data . This was particularly relevant for cranial curves where homologous correspondence is geometrically approximate rather than anatomically fixed. After GPA and sliding, the aligned landmark data were analysed using classical PCA to reduce dimensionality and identify dominant modes of variation. The principal components that accounts for $95\%$ of the variation explained were retained for further analysis.

2.2  Arc-Length Parameterisation

The arc-length parameterisation reparametrises a curve so that the speed is constant (usually 1). Given a landmark curve $f(t)$ with $t\in T=[0,1]$, its arc-length is:

[TABLE]

where $\dot{f}(t)$ is the derivative of $f$ with respect to $t$ and $\|\cdot\|$ is the Euclidean norm in $\mathbb{R}^{d}$, $d\geq 1$. The total length is $L=s(1).$

Arc-length parameterisation means defining a new parameter $u\in[0,L]$ such that:

[TABLE]

Here, $t(u)$ is the inverse of $s(t)$.

Arc-length reparameterisation presents numerous benefits: it mitigates the influence of variable speed in curve comparisons and guarantees uniform distribution along the curve.

2.3  Functional Data Framework

Landmark registration offers a straightforward approach to detect and align specific data points for each observation to the corresponding mean value. This provides a better representation of the mean in terms of amplitude variation. Each observation is vector-valued, as three spatial coordinates the $x$, $y$, and $z$ coordinates are involved.

To implement functional data in an object-oriented way, the raw data is converted into functions. For example, for the $x$-coordinates, using the observation points (boundaries) $\tau_{1},\dots,\tau_{N}$ and the set of discrete raw data (landmark values), a univariate functional data sample

[TABLE]

is constructed. This procedure is also applied to the $y$-coordinates and $z$-coordinates, yielding the functional samples

[TABLE]

respectively (see for instance Figure 2).

These three univariate functional datasets form a multivariate functional dataset, with $n$ outlines, each yielding a vector of $n$ observations defined over a 3-dimensional functional domain. The MFPCA package (happ2018multivariate, ) is used to perform the conversion to functional data (ramsay2005functional, ).

2.3.1  Multivariate Functional Principal Component Analysis

After acquiring the multivariate functional data, MFPCA is performed using the univariate functional principal components. The PCA basis functions are estimated from the $\mathbf{X}$. These basis functions are then applied to $n$ observations of $\mathbf{X}$ based on the PACE (PCA through conditional expectation) approach (yao2005functional, ).

Let

[TABLE]

be a vector-valued stochastic process corresponding to the functional random variables related to the standardised landmark $x$, $y$ and $z$ coordinates, respectively. For $p\in\{x,y,z\}$, let $I_{x}$ be a compact set in $\mathbb{R}$ with finite Lebesgue measure such that

[TABLE]

where $L^{2}(I_{x})$ denotes the space of square-integrable functions on $I_{x}$. Similarly, $(I_{y},L^{2}(I_{y}))$ and $(I_{z},L^{2}(I_{z}))$ are defined.

The $P$-fold Cartesian product is defined as

[TABLE]

so that $X$ is a stochastic process indexed by $t\in I$ and taking values in the product Hilbert space

[TABLE]

The inner product on $H$ is defined as:

[TABLE]

for $f=(f_{x},f_{y},f_{z})^{\top},g=(g_{x},g_{y},g_{z})^{\top}\in H$. Then $(H,\langle\!\langle\cdot,\cdot\rangle\!\rangle)$ is a Hilbert space with induced norm $\|\cdot\|$ (happ2018multivariate, ).

Multivariate Karhunen–Loève Representation.

Assume that

[TABLE]

Let $C$ denote the $3\times 3$ matrix-valued covariance function:

[TABLE]

where the $(p,q)$th element of $C(s,t)$ is:

[TABLE]

with $s_{p}\in I_{p},t_{q}\in I_{q}$; $p,q\in\{x,y,z\}$ and $X_{x}=X$, $X_{y}=Y$ and $X_{z}=Z$.

Define $\Gamma:H\to H$ as the covariance operator of $\mathbf{X}$:

[TABLE]

By the spectral theorem for Hilbert–Schmidt operators, there exists an orthonormal basis $\{\varphi_{j}\}_{j=1}^{\infty}\subset H$ and eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\dots\geq 0$ such that:

[TABLE]

The multivariate Karhunen–Loève representation is then:

[TABLE]

where $\xi_{j}=\langle\!\langle\mathbf{X},\varphi_{j}\rangle\!\rangle$ are zero-mean uncorrelated random variables with $\operatorname{Cov}(\xi_{j},\xi_{l})=\lambda_{j}\delta_{jl}$.

The truncated expansion for $J\geq 1$ is:

[TABLE]

and for the $p$-th component:

[TABLE]

where $\varphi_{p,j}$ are the univariate FPCA basis functions of $X_{p}$.

Multivariate Functional Principal Component Analysis in Practice

In practice, the eigenvalues $\{\lambda_{j}\}$, eigenfunctions $\{\varphi_{j}\}$, and scores $\{\xi_{j}\}$ must be estimated from $n$ i.i.d. observations $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(n)}$ of $\mathbf{X}$:

For each $X_{p}$, estimate the covariance function:

[TABLE]

and obtain $\hat{\varphi}_{p,j}$ and $\hat{\xi}_{p,j}^{(i)}$ for $j=1,\dots,J_{p}$. 2. 2.

Form the score matrix $\Xi\in\mathbb{R}^{n\times J}$ with $J=\sum_{p\in\{x,y,z\}}J_{p}$. 3. 3.

Estimate $Z\in\mathbb{R}^{J\times J}$ by:

[TABLE] 4. 4.

Perform an eigen-decomposition of $\hat{Z}$ to get $\hat{\lambda}_{j}$ and $\hat{v}_{j}$. 5. 5.

Compute the estimated multivariate eigenfunctions:

[TABLE]

2.3.2  Square-Root Velocity Function

In this approach, we apply the functional data morphometric (FDM) framework to Procrustes-aligned landmark data, incorporating the square-root velocity function (SRVF) representation of landmark curves. This method treats each 3D landmark trajectory as a multivariate function, enabling elastic shape analysis that respects the non-linear geometry of shape space.

The SRVF framework allows shape comparison by separating changes in geometry (shape) from changes in parameterisation (speed along the curve). In the FDM context, SRVF transforms 3D landmark trajectories into a form where elastic distances can be computed efficiently while maintaining translation invariance. Rotation invariance is ensured via GPA, and parameterisation invariance via arc-length reparameterisation.

Let each specimen be represented by a set of $N$ ordered landmarks in $\mathbb{R}^{3}$ over a parameter domain $T$. For specimen $i$, the landmark trajectory is:

[TABLE]

with

[TABLE]

After GPA, translation and scale effects are removed, yielding aligned trajectories.

Following Srivastava et al. (srivastava2011registration, ), the SRVF of a smooth curve $f$ is defined as:

[TABLE]

The Square Root Velocity Function (SRVF) representation possesses key properties that make it valuable for shape analysis: it is naturally invariant to translations since differentiation removes any constant shifts; the $L^{2}$ distance between SRVFs corresponds to the elastic geodesic distance in shape space, enabling an efficient comparison of shapes through standard inner products. However, rotation and parameterisation invariances are not inherent and must be explicitly enforced. Rotation invariance is typically obtained by aligning SRVFs via Procrustes superimposition, while parameterisation invariance is achieved by reparameterizing all curves to a common domain using uniform arc-length parameterisation.

Once the SRVFs $\{q^{(i)}\}_{i=1}^{n}$ are obtained, they are treated as multivariate functional observations:

[TABLE]

These are then analysed using the FDM pipeline: conversion into a functional data object, smoothing, and application of multivariate functional PCA (MFPCA) for dimensionality reduction and shape variation analysis.

3  Simulation Study

We conducted a simulation study to assess the innovations driven by functional data in the morphometric analysis of three-dimensional shapes. Specifically, we compared eight pipelines: GM, arc-GM, FDM, arc-FDM, elastic-SRV-FDM, arc-elastic-SRV-FDM, soft-SRV-FDM, and arc-soft-SRV-FDM.

3.1  Simulation Design

Four distinct groups were generated based on an initial helix configuration: Group 1 maintained the original helix, Group 2 applied sinusoidal perturbations to the $x$- and $y$-coordinates, Group 3 introduced a sinusoidal perturbation to the $z$-coordinate, while Group 4 incorporated a phase shift into the $x$- and $y$-coordinates, leaving the $z$-coordinate unaltered. The simulations involved four groups with sizes $(21,\,32,\,124,\,23)$, and their proportions were carefully selected to align with those observed in kangaroo data. Each curve was defined over $t\in[0,1]$ and was sampled at $n_{\mathrm{points}}=30$ uniformly spaced intervals,

[TABLE]

The simulated coordinates are then: Group 1

[TABLE]

Group 2:

[TABLE]

Group 3:

[TABLE]

Group 4:

[TABLE]

where $\phi$ denotes a small phase shift applied to the helix.

Each specimen underwent an independent monotone time-warp,

[TABLE]

to mimic realistic phase variation in curve parameterisation. This step reflects natural variability in biological shape data, where features were locally stretched or compressed, and tests the robustness of the pipelines to misalignment. In addition, independent and identically distributed Gaussian noise was added to each coordinate with group-specific standard deviations $(\sigma_{G1},\sigma_{G2},\sigma_{G3},\sigma_{G4})=(0.05,\,0.05,\,0.10,\,0.06)$. Simulations were performed over $n_{\mathrm{reps}}=50$ independent replicates.

Pipeline Settings.

• •

Arc-length reparameterisation: 30 output points.

• •

Geometric morphometrics (GM): 30 equally-spaced pseudo landmarks aligned by GPA.

• •

FDM basis: cubic B-splines on $[0,1]$ with 10 basis functions.

• •

Multivariate FPCA (MFPCA): target $M=30$ components with safe downgrading if needed.

• •

Component selection: smallest $k$ such that cumulative variance explained $\geq 0.95$.

• •

Elastic/SRVF alignment (for elastic and Soft-FDM variants): rotation- and scale-invariant, with $\alpha_{\mathrm{soft}}=0.6$ and $\lambda_{\mathrm{soft}}=0.01$.

Evaluation.

We used a five-fold cross-validation scheme for evaluation. (Alignment, basis construction, and component extraction were recomputed separately for the training and test folds. Classifiers were fitted on scores from the training folds and applied to scores obtained from the corresponding test decompositions.) Classification was performed on the resulting principal component scores of each pipeline using:

• •

linear discriminant analysis;

• •

multinomial logistic regression;

• •

linear SVM (untuned, with centering and scaling).

3.2  Simulation results

Figure 3 shows differences in classification performance, where functional data–driven approaches achieved higher predictive accuracy and better stability between folds. Figure 4 illustrates that functional pipelines reached the 95% variance threshold with fewer components than coordinate-based methods (Table 1). The corresponding reconstruction plot (Figure 6) confirmed that the reconstructed shapes using the retained components preserved the curved structural information, thus validating the use of the 95% cutoff. Both elastic and soft-aligned pipelines outperformed the purely coordinate-based and standard functional approaches. The elastic variants generally achieved good accuracy and efficient variance capture, but their performance showed greater variability across folds, indicating reduced stability. In contrast, soft-aligned pipelines offered consistently better accuracy with reduced variability and clearer group separation in PC space (Figure 7).

4  Application to 3D Kangaroo Data

4.1  Kangaroo Cranial Landmark Dataset

To analyse cranial shape variation in extant kangaroo species, we employed two parallel workflows: one based on GPA with sliding semi-landmarks, and another based on raw (unaligned) 3D landmark configurations. This dual approach allowed us to explore the effects of shape alignment and preprocessing on morphometric outcomes.

4.1.1  Landmark Acquisition and Replicate Averaging

A dataset of 3D cranial landmarks was obtained from the publicly available Butler Kangaroo Landmark collection (butler20213d, ). Each specimen was represented by 48 anatomical landmarks digitised across 41 extant specimens. Each specimen has been assigned with dietary category labels (browser, grazer, mixed feeder, omnivore) obtained from the classification file.

Each specimen was landmarked twice. Bland–Altman analysis showed negligible differences between replicates 1 and 2. Therefore, the aligned replicates were averaged to produce a single landmark configuration per specimen for further analysis.

The kangaroo dataset (butler20213d, ) already includes a mixture of fixed landmarks and semi-landmarks digitised along cranial curves. These semi-landmarks were standardised using the established sliding procedure, which minimises bending energy relative to the Procrustes consensus. This ensured comparability with prior morphometric analyses of the dataset.

To remove non-shape variation, raw landmark configurations were first subjected to GPA. For functional pipelines, landmark curves were interpolated and reparameterised either by time (landmark index) or by cumulative arc-length. The SRVF was then applied to enable elastic shape alignment, which was performed in both strict (elastic) and regularised (soft) forms.

Eight pipelines were implemented: GM, FDM, arc-GM, arc-FDM, elastic-SRV-FDM, soft-SRV-FDM, arc-elastic-SRV-FDM, and arc-soft-SRV-FDM. Each was subjected to PCA or MFPCA depending on the representation, and the number of retained principal components was determined by a 95% cumulative variance threshold. The landmark configuration of a representative specimen (Specimen 1) was then reconstructed from the truncated components for each pipeline. In GM and Arc-GM, reconstruction added the rank-$k$ PCA approximation to the Procrustes consensus. In FDM and arc-FDM, the truncated MFPCA expansion from the functional representation was assembled across the $x$, $y$, and $z$ coordinates. For elastic and arc-Elastic pipelines, reconstruction was performed in SRVF space by adding the specimen’s truncated SRVF-PCA scores to the SRVF mean, followed by inverse SRVF mapping via time integration and resampling to the landmark grid. Soft-SRV-FDM and arc-soft-SRV-FDM used analogous truncated MFPCA predictions under soft alignment. Reconstructed shapes were Procrustes-superimposed to their corresponding originals for visual comparison, and reconstruction accuracy was quantified using mean squared error (Figure 10).

To evaluate the discriminatory capacity of each pipeline, the retained PC scores were used as predictors in three classification methods: linear discriminant analysis (LDA), multinomial logistic regression, and linear-kernel support vector machines (SVM).The supervised learning methods are trained using the retained PC scores.

All analyses were conducted in R version 4.3.2)RStudioTeam2023 . Landmark-based morphometric procedures, including GPA, were implemented using the geomorphgeomorph ; geomorphadams2021 package. Functional data representations and elastic alignment were carried out using fdasrvffdasrvf , while MFPCA was performed with MFPCAmfpca and funDatafundata packages. Classification models were fitted using caretcaret , which provided a unified framework for LDA, multinomial logistic regression, and SVM with linear kernels. Additional packages employed included Morphomorpho , MASSmass , tidyrtidyr , dplyrdplyr , ggplot2ggplot2 , and rglrgl for statistical computation, data wrangling, and visualisation.

4.2  Results

The number of principal components required to capture 95% of the variance varied substantially across pipelines (Table 3). In all cases, subsequent components contributed progressively smaller proportions of variance, with a sharp drop after the first two to three components (Figure 8). This pattern was most pronounced in FDM and Soft-FDM, where the first three components alone accounted for nearly all of the variance. Similarly, arc-soft-FDM and elastic-SRV-FDM reached the 95% threshold with only four components. By contrast, arc-GM and GM required more components, while arc-elastic-SRV-FDM was the most component-intensive, requiring 15 components to achieve the same level of variance explanation.

Classification accuracies varied across pipelines and methods. Elastic-SRV-FDM consistently outperformed other pipelines (Figure 9), achieving the highest accuracy across classifiers, particularly with LDA. In contrast, GM, arc-GM and soft-SRV-FDM pipelines yielded comparable accuracies, while arc-soft-SRV-FDM yielded the lowest classification performance.

Plots of the best PC pairs for each pipeline revealed varying degrees of dietary group separation and within-group tightness (Figure 7. Overall, elastic-SRV-FDM and soft-SRV-FDM pipelines capture the largest proportion of shape variation, that is, over 85% in the first two MFPCs. They also exhibit the most distinct groupings, especially for browsers and omnivores. While the clusters were slightly looser than those in the elastic-SRV-FDM analysis, browsers and omnivores again exhibited tighter cohesion, and the overlap between grazers and mixed feeders was reduced compared to GM and FDM. In contrast, the arc-FDM and arc-soft-SRV-FDM pipelines relied on higher-order PCs that explained less than 2% of the variance in each axis, revealed dispersed and intermixed.

5  Discussion

Functional data morphometrics is advantageous because it captures smooth variation along entire curves or surfaces. The analysis of the kangaroo cranial dataset showed that the elastic SRV-FDM pipeline performed best, followed by the two GM pipelines. This outcome is expected, as the dataset is purely landmark-based, which inherently favours GM methods. Since curvature is not represented in such datasets, the advantages of functional data approaches are limited. An important property of the elastic SRVF framework is its ability to preserve amplitude variation while aligning phase srivastava2011registration , which in this case is effectively synonymous with aligning the geometric structure of a landmark-only dataset. Consequently, the nature of this dataset does not favour soft-aligned functional data approaches, since the central idea of balancing the contributions of landmarks and continuous function shapes to reach an optimal registration guo2022data cannot be realised when only landmarks are available.

The simulation study was designed to address these shortcomings by incorporating phase variation and moderate coordinate noise, to assess how the eight pipelines perform under noisy, imbalanced folds. The soft-alignment approach of Guo et al. (2022) guo2022data has an advantage in this context, serving as a compromise that pushes landmarks closer without requiring exact overlays, while balancing contributions of landmarks and functional shapes to achieve an optimal registration.

The combination of arc-length parameterisation with SRVF alignment also provides a powerful framework for morphometric pipelines. The SRVF framework (after arc-length normalisation) simplifies the originally complicated elastic metric into a straightforward $L^{2}$ product srivastava2011registration . For the kangaroo data, arc-GM appears to perform better in terms of accuracy compared to its GPA-aligned GM counterpart, because arc-length parameterisation eliminates the influence of uneven landmark spacing or irregular sampling along the curve. When more balanced shapes are generated in the simulation study, the advantage of combining arc-length parameterisation and SRVF becomes evident, with pipelines implementing these aspects showing improved accuracy.

The soft-registration approach is likely to have an edge when images are processed from scratch, as it allows the relative contributions of functional shape and landmarks to be determined directly from the data itself guo2022data . However, its effectiveness is not evident in the kangaroo dataset, but is better illustrated in the simulated data, since the former is purely landmark-based.

For reconstruction, the number of principal components required to capture 95% of variation was used for each pipeline. The data application and simulation results in general show that pipelines combining arc-length parameterisation with SRVF alignment achieve strong reconstruction accuracy, as reflected by low mean MSE values (Table 2 and Table 4), thereby reinforcing the strength of elastic functional frameworks in handling phase variability. Notably, FDM also performs particularly well, yielding the lowest MSE and underscoring its efficiency in reconstruction tasks where smooth functional representation is advantageous. Principal component analysis and MFPCA maximise variance rather than discrimination. Consequently, the first two PCs capture the dominant variation, but it is important to note that this may not correspond to group separation but may reflect global size or noise. Discriminatory structure can reside in later PCs with smaller variance, as illustrated in the data application.

6  Conclusion

Our work compared landmark-based and functional data–driven pipelines for the analysis of three-dimensional cranial landmark data, using both simulations and real kangaroo cranial datasets to evaluate variance explained, reconstruction accuracy, and classification performance. Across both datasets, functional data approaches, and in particular arc-parameterised and SRVF-based methods, outperformed coordinate-based pipelines by preserving biologically meaningful structural information and improving dietary classification accuracy. While FDM yielded the strongest reconstruction accuracy, elastic SRVF-based methods provided superior classification accuracy, highlighting complementary strengths across pipelines. The present work was based on three-dimensional landmark coordinate data, which neglects curvature, where this limitation is addressed through the simulated data. For future studies, it would be valuable to extend the analysis to full three-dimensional images that capture both raw coordinates and curvature, thereby enabling a more robust evaluation, particularly when applying soft-aligned and elastic pipelines. Arc-length parameterisation and SRVF alignment also show great potential for advancing morphometric frameworks, and these innovations could additionally be incorporated into our two-dimensional functional data morphometrics framework pillay2024functional to enhance 2D analyses. Overall, these directions highlight the advantages of functional data analysis for complex morphological datasets, where elastic and soft alignment strategies, coupled with arc-length parameterisation, can better uncover ecologically meaningful patterns often obscured under traditional landmark-based approaches. The eight pipelines presented here provide a flexible framework for selecting methods best suited to the characteristics of the data.

Disclosure Statement

During the preparation of this manuscript, the authors used ChatGPT (OpenAI) to improve the academic tone and enhance clarity of presentation. The tool was also used to reduce potential coding errors, improve readability, and ensure tidy code. All methodological design, analyses, and scientific interpretations were conducted solely by the authors. All content generated or suggested by the tool was critically reviewed, edited, and validated by the authors, who take full responsibility for the final version of the manuscript.