2.1. Subject data and ethics

Secondary data analysis ethical approval was granted by the University of Southampton’s Ethics Committee (ERGO II 65748A2). Subject data consisted of MRI or CT scans of 43 people with residual limbs from transtibial amputations using five different sources, collected in previously published research, and/or provided to the authors under data sharing agreements (Steer, Reference Steer2019; Bramley, Reference Bramley2020; Edgar et al., Reference Edgar, Berry, Moes, Adolphi, Bridges and Nolte2020; Mendis, Reference Mendis2021; Ding et al., Reference Ding, Henson, Sivapuratharasu, McGregor and Bull2023; Finco et al., Reference Finco, Finnerty, Ngo and Menegaz2023). The participants covered a range of amputation causes, age, and time since amputation (Table 1). Six scans were excluded due to failing to satisfy at least one of the following inclusion criteria, such that the model would describe a cohesive population and consistent anatomical features:

• • Amputation through the tibial diaphysis, that is, excluding e.g. Symes amputation.

• • The limb and scan must have the typical residual skeletal anatomy for a transtibial amputation, that is, distal femur, patella, proximal tibia, and proximal fibula should be present.

• • Standard surgical amputation technique must have been used for the amputation, without additional reconstruction, that is, excluding bone-bridging or capping or implants for prosthesis suspension.

Table 1. Details of the included individuals

Note. Four of the obtained datasets were excluded after building a preliminary SSM and their details are left out of this table.

^a For NMDID datasets (Edgar et al., Reference Edgar, Berry, Moes, Adolphi, Bridges and Nolte2020; Finco et al., Reference Finco, Finnerty, Ngo and Menegaz2023), weight and time since amputation were not provided.

^b Two participants had bilateral amputations and only their right residual limb was used in the model.

^c NMDID did not report reasons for amputation, however, of these 18 individuals three had a listed Type I diabetes diagnosis, and 10 Type II diabetes, so the vascular disease cause is likely to be an under-estimate.

Four more scans were excluded from the dataset after producing a preliminary model which identified very substantial imaging artefacts, as explained next.

2.2. Residual limb processing

2.2.1. Mesh generation and alignment

The raw MRI and CT DICOM image sets representing transverse slices were segmented to provide corresponding surface meshes of the skin and residual tibia, residual fibula, patella, and distal femur bones, in .stl format (ScanIP Version 2018.12; Synopsys, Inc., Mountain View, USA).

The MRI and CT scans were taken with the patient supine with the posterior aspect of their residual limb supported. While this gives an approximately standardized pose, the orientation still varied between limbs. By aligning the training datasets, variations in relative position and orientation can be removed such that the SSM is able to describe more detailed anatomic shape variation. Due to differences in remaining anatomy and intersubject variation in knee joint alignment, the alignment of the proximal residual tibias was prioritized.

First, each limb was aligned with respect to a global coordinate system (CS) defined by the International Society of Biomechanics convention (xz plane = transverse, yz = coronal, xy = sagittal) (Wu et al., Reference Wu2005) using the ampscan open-source shape analysis toolbox (Steer et al., Reference Steer, Stocks, Parsons, Worsley and Dickinson2020). The alignment procedure was as follows (Figure 1):

• • left-sided shapes were mirrored so that all subjects were the same side,

• • shapes were then translated so that their tibia mesh vertex centroids lay at the origin,

• • shapes were rotated to align their tibia mesh principal axes of inertia, estimated using PCA, with the global coordinate axes,

• • finally, datasets were rotated about the vertical axis so that the posterior extents of the knee condyles lay along a medial–lateral axis.

Figure 1. Coronal and sagittal views of the first 10 subjects’ training data (full and bones only) meshes before (top) and after the alignment and trimming were performed (bottom).

Second, an iterative closest point (ICP) rigid registration algorithm was used to adjust each dataset’s alignment with respect to the tibia, to maximize intersubject alignment, followed where necessary with small manual adjustments by an experienced human observer. After alignment, a trim was applied in the transverse plane at a level one patella height above the most proximal extent of the patella, to include a consistent degree of geometry above the knee relevant to supracondylar socket designs.

2.2.2. Size normalization

SSMs are often created using size-normalized training data, to separate size and shape variance. However, due to variations arising from surgery, each of the training shapes’ residual limb length represents a different proportion of the intact limb, so a simple scaling factor is not sufficient to ensure the mapping of corresponding anatomical features (Lam et al., Reference Lam, Garrison and Rozbruch2016). For example, a taller individual with a relatively high amputation may have a longer residual limb than a shorter individual with much more intact anatomy. Therefore, this study size normalized the training shapes to the full tibia length, $ T $ as the unit size. This was taken from the contralateral limb images where available, or estimated from the individual’s height, $ H $ (mm), and their age, $ A $ (years), using a regression formula (Eq. 1) (Trotter and Gleser, Reference Trotter and Gleser1952):

(1) $$ {\displaystyle \begin{array}{l}{T}_{male}=\frac{H-786.2+0.6\alpha }{2.52},\\ {}{T}_{female}=\frac{H-615.3+0.6\alpha }{2.90}\end{array}} $$

where α is defined in

(2) $$ \alpha =\left\{\begin{array}{c}0\; if\hskip0.24em A\le 30\\ {}\left(A-30\right) if\;A>30.\end{array}\right. $$

Following size normalization a visual check of alignment was made

2.2.3. Registration

Mesh registration is necessary for point-to-point comparison between each of the SSM’s training shapes so that they are represented using a corresponding set of vertices, allowing their shape variance to be extracted. A “baseline” shape was chosen and each of its mesh bodies mapped onto the surface of the corresponding “target” data subject bodies using Amberg and Romdhani’s (Amberg et al., Reference Amberg, Romdhani and Vetter2007) nonrigid ICP algorithm in the trimesh package (Dawson-Haggerty, 2019). This algorithm was selected as it is less sensitive to outliers and missing data than conventional elastic matching methods. This deforms the baseline toward the target incrementally, permitting accurate registration of anatomy present in both datasets, and appropriately distributed mesh vertices in cases where the anatomy in the baseline mesh is not present on the target dataset. Figure 2 depicts an example of the key steps applied to one of the skin meshes, including over scaling the baseline to the target to ensure registration of the model’s open proximal boundary.

Figure 2. The main steps in registration (original, aligned, scaled, and registered meshes, from left to right).

2.3. Statistical analysis

Principal component analysis (PCA) is considered the standard approach used for creating statistical shape models (Stegmann and Gomez, Reference Stegmann and Delgado Gomez2002). PCA objectively identifies patterns of shape variance (modes) in a set of training data. This analysis enables the dimensionality of the population’s characteristics to be reduced by selecting a limited number of important modes of variation, that is, principal components (PCs) (Audenaert et al., Reference Audenaert, Pattyn, Steenackers, De Roeck, Vandermeulen and Claes2019). PCA describes how each training dataset compares to the population mean shape, described by a “score” in each mode of variation. Two statistical limb shape models were produced by PCA, representing the skin and bone surfaces (“full model,” 𝜇) and the skin surface only (“skin model,” $ {\mu}^s $ ), using the scikit-learn toolbox (Abraham et al., Reference Abraham2014). The procedure was an extension to method which was previously reported for analyzing the residual limb’s external shape (Dickinson et al., Reference Dickinson, Diment, Morris, Pearson, Hannett and Steer2021).

First, the mesh vertices $ {\boldsymbol{x}}_i $ for each $ i $ of the $ n $ training individuals were each represented as a column vector, Eq. (3), where $ m $ represents the number of vertices in the baseline mesh and $ x $ , $ y $ , and $ z $ ; are the vertex coordinates,

(3) $$ {\boldsymbol{x}}_i={\left[{x}_1,{y}_1,{z}_1,\dots, {x}_m,{y}_m,{z}_m\right]}^T. $$

A column vector representing the mean limb shape was then calculated using, Eq. (4),

(4) $$ \overline{\boldsymbol{x}}=\frac{1}{n}\sum \limits_{i=1}^n{\boldsymbol{x}}_i, $$

where $ n $ represents the number of limbs in the training dataset, giving Eq. (5),

(5) $$ \overline{\boldsymbol{x}}={[{\overline{x}}_1,{\overline{y}}_1,{\overline{z}}_1,\dots, {\overline{x}}_m,{\overline{y}}_m,{\overline{z}}_m]}^T; $$

and the meshes were mean centered according to Eq. (6),

(6) $$ {\boldsymbol{x}}_{i{\hskip2pt }_{mc}}={\boldsymbol{x}}_i-\overline{\boldsymbol{x}}. $$

PCA by singular value decomposition (SVD) was then performed. This enabled each limb shape to be described by Eq. (7), where $ \boldsymbol{A}j $ represents the eigenvectors corresponding to the each of the c modes of population shape variation,

(7) $$ {\boldsymbol{x}}_i=\overline{\boldsymbol{x}}+\sum \limits_{j=1}^c\boldsymbol{A}j\boldsymbol{\lambda} {j}_i. $$

and $ \boldsymbol{\lambda} j $ is a vector of weighting coefficients or “mode scores,” associated with the eigenvectors to describe each training dataset shape’s deviation from mean, with $ \boldsymbol{\lambda} j $ being represented by Eq. (8),

(8) $$ \boldsymbol{\lambda} j=\left[{\lambda j}_1,\dots, {\lambda j}_n\right]. $$

The approximated 95% range of population variation was plotted by generating synthetic shapes ( $ \eta $ ) with weighting coefficients defined using the training dataset’s 2.5th and 97.5th percentile mode scores for each mode $ j $ in turn using Eqs. (9 and 10), allowing visual interpretation of the shape variance contained within the mode,

(9) $$ {{\boldsymbol{\eta}}_j}_{0.025}=\boldsymbol{A}{\boldsymbol{\lambda} j}_{.025}+\overline{\boldsymbol{x}}, $$

(10) $$ {{\boldsymbol{\eta}}_j}_{0.975}=\boldsymbol{A}{\boldsymbol{\lambda} j}_{.975}+\overline{\boldsymbol{x}}. $$

During model generation, all training shapes’ mode scores were inspected. Four shapes were observed to have very distorted soft tissues associated with positioning in the MRI or CT scanner, identified by outlying influence in both extremes of leave-one-out cross-validation tests, so they were set aside leaving a training dataset of 33.

2.4. Statistical shape model performance and validity testing

2.4.1. Compactness: how many modes of variability does the model require?

As PCA is a dimensionality reduction technique, the compactness describes the ability of the SSM to represent the population variation with a minimal subset of modes (Pearson, Reference Pearson1901). Compactness was assessed by reconstructing each training individual ( $ i $ ) with progressively increasing number of modes ( $ j $ ) and calculating the root mean square error (RMSE or $ {\varepsilon}_j $ ) of vertex surface deviation between the reconstruction ( $ {{\overset{\sim }{\boldsymbol{x}}}_i}_j $ ) and the original data ( $ {\boldsymbol{x}}_i $ ), such that Eq. (11) represents each reconstruction

(11) $$ {\displaystyle \begin{array}{c}{\tilde{\boldsymbol{x}}}_{i_0}=\overline{\boldsymbol{x}}\\ {}{\tilde{x}}_{i_1}=\overline{\boldsymbol{x}}+{A}_1{\lambda}_{1_i}\\ {}{\tilde{\boldsymbol{x}}}_{i_2}=\overline{\boldsymbol{x}}+{A}_1{\lambda}_{1_i}+{A}_2{\lambda}\\ {}\vdots \\ {}{\tilde{\boldsymbol{x}}}_{i_c}=\overline{\boldsymbol{x}}+{A}_1{\lambda}_{1_i}+{A}_2{\lambda}_{2_i}\cdots + Ac\lambda {c}_i\end{array}} $$

and Eq. (12) represents the error

(12) $$ {\displaystyle \begin{array}{c}{\varepsilon}_{n_0}={\boldsymbol{x}}_i-{\tilde{\boldsymbol{x}}}_{i_0}\\ {}{\varepsilon}_{n_1}={\boldsymbol{x}}_i-{\tilde{\boldsymbol{x}}}_{i_1}\\ {}\vdots \\ {}{\varepsilon}_{n_c}={\boldsymbol{x}}_i-{\tilde{\boldsymbol{x}}}_{i_c}\end{array}} $$

2.4.3. Generality by recreation: how accurately can the model describe left-out shapes?

Generality was further assessed by evaluating how well the SSM can describe a new data set not used in its creation by using each LOO-SSM model from the cross-validation test to recreate the left-out shape. The left-out shape’s mode scores $ \lambda $ ( $ m\times 1 $ column vector) were estimated by solving the least squares matrix problem. Conventionally this takes the form $ \boldsymbol{A}\boldsymbol{x}=\boldsymbol{b} $ , though we are prioritizing nomenclature convention used in SSM studies, in Eq. (13),

(13) $$ \boldsymbol{A}\boldsymbol{\lambda } =\boldsymbol{x}, $$

where $ \boldsymbol{A} $ is the $ m\times c $ matrix representing the principal components ( $ Aj $ ) and $ \boldsymbol{x} $ is the $ m\times 1 $ column vector representing the left-out shape’s vertex coordinates, see Eq. (2). Note $ c $ , the number of modes (PCs), will be one less in the LOO-SSMs than in the full SSM because the number of modes that can be determined is equal to $ n-1 $ . As $ \boldsymbol{A} $ is not a square matrix it is not invertible, therefore, the Moore–Penrose pseudoinverse was computed to solve Eq. (13) (Generalized Inverses, 2003), using NumPy’s linear algebra submodule.

2.5. Example use case: Shape prediction from partial data

A common use of SSMs is prediction of a shape given only partial data. An application of the present model might be to predict the internal geometry from the external shape. Two methods were compared.

First, a similar pseudoinverse method to that used earlier for model recreation was tested. For this purpose, a consistent partial shape, the skin, and consistent missing data, the bones, were used. This means that point-to-point correspondence between the partial shape and the SSM was known. This method is also known as Gappy POD (proper orthogonal decomposition) (Bui-Thanh et al., Reference Bui-Thanh, Damodaran and Willcox2004). In place of Eq. (13), Eq. (14) provides,

(14) $$ {\boldsymbol{A}}^s\boldsymbol{\lambda} ={\boldsymbol{x}}^s, $$

where $ {\boldsymbol{A}}^s $ is the subset of $ \boldsymbol{A} $ that only contains the first $ s $ rows that act on the skin nodes, and $ {\boldsymbol{x}}^s $ is the SSM skin mesh registered to the partial shape. Employing the earlier method, mode scores can be estimated and then used with the complete principal components to predict a complete shape.

An alternative method employing linear regression was also tested (Figure 3). Training dataset mode scores calculated from each leave-one-out skin only model were collected alongside the corresponding mode scores from each leave-one-out full model, and a linear regression between each pair of mode scores was formulated using scikit-learn (Abraham et al., Reference Abraham2014), such that the full mode scores were the dependent variable. The recreation method, provided in “Generality by recreation,” was used upon the skin-only leave-one-out model to estimate mode scores and thus recreate each left out shape’s skin. The corresponding leave-one-out regression was then fitted to the skin-only mode score estimates to predict mode score values for the full model’s description of the left-out shape. Using these predicted scores, a predicted shape was generated using Eq. (8) and plotted. This was then aligned with the actual left-out shape and the RMSE calculated. In addition to mean centering and scaling, the PCA used scikit-learn’s “whitened” components scores to minimize linear correlation across features.

Figure 3. Flowchart of steps in linear regression prediction method.

The prediction accuracy was calculated with both methods as the RMSE between the actual and predicted shapes for (i) the bones, (ii) the skin, and (iii) the whole limb. The RMSE was calculated with the normalized data and expressed in mm by rescaling back to the actual size. Differences in accuracy between the prediction methods were assessed for statistical significance with a nonparametric Wilcoxon signed rank test, for paired data.

3.1. Modes of variation

The participants’ known or estimated full tibia lengths indicated that the training shapes represented a median 36% of the intact anatomy. The proportionally shortest and longest limbs preserved 18% and 60% of the intact tibia length, respectively (Appendix 1). After size normalization, the average locations of vertices in the aligned, registered training shape meshes were calculated to produce the population mean shape (Figure 4). The PCA calculation provides the proportion of population variance attributed to each mode of shape variation (Figure 5) and indicates that 69% of the population variance was contained within the first 2 modes and 95% within 10 modes, implying that these are most important for describing or classifying gross shape.

Figure 4. Mean of residual limb training shapes.

Figure 5. Individual and cumulative variance represented by each mode of the full statistical shape model.

The shapes of these nine primary independent modes of limb variation were plotted to permit inspection of the variance they represent (Figure 6). The variations predominantly manifested as amputation height in Mode 1 which encompasses 52.2% of the population variance, and slenderness/soft tissue bulk in Mode 2, including 16.7% of the variance. As discussed later, these external shape variations are consistent with previous SSMs. However, this model provides first quantitative insights into how the internal and external anatomic shape variation is related. For example, Mode 3 (8.2% of variance) represents the relative scale of the bones within the soft tissue, and along with modes 5 (3.1%) and 6 (2.6%) includes varying distal soft tissue coverage over the distal tibia bony prominence, a notable site of particular soft tissue vulnerability and relevance in socket design. These modes illustrate how a relatively shorter residual bone allows the posterior soft tissue more freedom to deform posteriorly and medially. The model also shows variance which might be associated with scanning; Mode 4, 7.2% of the population variance, describes knee flexion, and Modes 5 (3.1) and 8 (1.3%) describe posterior and medial variation in soft tissue distribution which might be associated with soft tissue deformation under the body’s self-weight and contact with the scanner bed, respectively, in the calf and thigh.

Figure 6. SSM mode shapes as described by 2.5th (blue) to 97.5th percentile (red) estimated variance range from in the training dataset. These permit the principal modes of residual limb shape variance to be inspected.

3.2. Validity tests

Compactness testing (Figure 7) demonstrated that by employing the initial 10 modes which encompass 95% of the variance, all training shapes could be reconstructed with an RMSE below 5 mm. In context, the training shapes were a median of 14 mm RMSE (range 4–30 mm) from the mean shape, indicated by their reconstruction without using any SSM modes. In cross-validation testing (Table 2), the mean and 2.5th and 97.5th percentile extreme shapes in each mode were generated for each of the LOO-SSMs and compared to the corresponding shapes obtained from the full SSM. Leaving out training shapes had little influence on the mean, at 0.14–0.88 mm RMSE surface deviation (median 0.42 mm). The RMSEs associated with leaving out individual shapes were also small for the extreme shapes in Modes 1 and 2 (medians of 0.36–0.70 mm, and interquartile ranges [IQR] of 0.66–1.06 mm, respectively), however, given the limited training dataset size, individual shapes could disproportionately affect single modes. This is evidenced by two outliers in Mode 1, represented by training datasets 10 and 20, which were considerably longer than the rest (Table A1); leaving each of them out influenced the long-length extreme shape by 8.29 and 8.97 mm, respectively, but did not substantially influence the short length extreme.

Figure 7. Compactness: range of RMSE between each subject’s actual shape and its reconstruction from the SSM using a reduced number of modes. Normalized data rescaled back to actual size for expression in mm.

Table 2. Variation of the leave-one-out statistical shape model (LOO-SSM) for mean and first two mode shapes compared to the full SSM

Note. Normalized data rescaled back to actual size for expression in mm. RMSE: root mean squared error; IQR: interquartile range.

3.3. Prediction of bones from residual limb surface

The model’s ability to predict each residual limb’s bony anatomy from its external surface was evaluated and compared to the model recreation (Figure 8). Four shapes are presented to illustrate relatively good and poor examples of the model’s predictive capabilities. Participant IDs 1 and 2 (Figure 8, rows 1 and 2) are both near the population mean in Mode 1. Both methods predicted the soft tissue shape well, and displayed some inaccuracies, particularly with respect to predicting individual bony anatomy characteristics. For example, ID1 has a relatively short residual fibula, and ID2 has a quite lateralized fibula and femur. The PI method showed higher error in predicting soft tissue coverage over notable bony prominences such as the fibula head and distal tibia, which is potentially problematic for its use in socket design or analysis. The linear regression (LR) method produced more accurate bone shapes and slightly lower error in all the locations listed earlier. Subject ID8 is included as an example which was predicted poorly by both methods (Figure 8, row 3) as it features unique bony characteristics which were far from the population mean and not indicated by the external soft tissue shape: the very short fibula and the relatively short tibia compared to the full limb length, represented as high tissue thickness at the limb’s distal end. The other main weakness in prediction is demonstrated by Subject 10 (Figure 8, row 4) who has a considerably longer residuum (60% of intact tibia) compared to the rest of the training dataset. This is difficult to predict because it features anatomy not present in the remainder of the training data. Relative to the recreation, the PI method performs particularly well especially in the prediction of distal soft tissue coverage, however, both methods are limited by the model’s inability to describe the longer limb shape and very slender distal bony anatomy.

Figure 8. Four example subjects’ original (red) compared to recreated (blue) shapes, and shapes predicted using pseudoinverse (PI) and linear regression (LR) methods alongside error measurements (mm). Normalized data rescaled to actual size for expression in mm; % in row headers refers to proportion of intact tibia remaining.

Quantifying the overall errors for predicting the full shapes (Figure 9), both the linear regression and pseduoinverse methods give a similar spread of values (LR: median 4.97 mm, IQR 1.94 mm; PI: median 5.42 mm, IQR 1.50 mm), but the LR method outperforms the PI method overall (p = .012, effect size = −0.4). This trend became stronger considering the bone prediction alone, approaching the RMSE values for recreation (LR: median 6.66 mm, IQR 2.60 mm; PI: median 8.58 mm, IQR 2.84 mm; p < .001, effect size = −0.8). Since recreation shows the model’s accuracy in describing the dataset given full information, it can be considered as the lower achievable bound for prediction error values (all: median 3.40 mm, IQR 0.99 mm; bones: median 3.46 mm, IQR 1.21 mm). Similarly, looking at the distal tissue thickness (Figure 10), LR predictions fall closer to the actual values than the PI predictions on average, but this did not reach significance (p = 0.272, effect size = 0.2).

Figure 9. RMSE of the recreated and predicted shapes compared to the actual shape. Normalized data rescaled back to actual size for expression in mm.

Figure 10. Range of distal tissue thickness error (measured as most distal point of the tibia to most distal point of the skin compared to the original shape) for the recreated and predicted shapes. Normalized data rescaled back to actual size for expression in mm.