1. Introduction
Pes planus of the foot, an idiopathic [
1] condition commonly seen in children which may continue to adulthood, is associated with a low or absent medial longitudinal arch and some degree of hindfoot and forefoot malalignment [
2]. The true incidence of flatfoot is unknown; however, many research studies have tried to understand this condition [
3,
4,
5,
6,
7,
8], with conflicting outcomes contributing to a general disagreement even on its explicit definition and the differences between a “normal” foot and a flatfoot [
9]. For this reason, flatfoot is often considered as a normal variant of the foot shape rather than a condition [
9]. There are mainly three types of flatfeet [
10]: the flexible flatfoot (FFF), which is usually asymptomatic and most common, the flexible flatfoot with short Achilles tendon (FFF-STA), which affects one in four people with FFF, causing pain and often requiring surgery, and the rigid flatfoot, which is always symptomatic and caused by neuromuscular disorders [
11]. The term “flexible” means that the foot appears flat when it is weight-bearing (when the patient stands up), while it returns to normal when the foot is no longer weight-bearing (when the patient sits down with their feet suspended) [
12].
Corrective surgery in FFF is recommended in adolescents only for symptomatic cases to prevent disability and arthritis developing in adulthood [
13] when other conservative methods such as orthotics have failed to relieve pain [
9,
14]. However, surgery brings its own challenges. Given the complex foot structure and the not-fully understood foot biomechanics, there is no universal agreement on how to treat flatfoot successfully, and no optimal surgical methodology has yet been established [
15]. Surgical planning needs to consider many factors that are not well defined to remodel the foot accordingly, relying on educated decisions from foot surgeons that are based on empirical evidence and experience. Patient outcomes can be varied [
16,
17,
18,
19], and therefore, there is some disagreement about optimal surgery for individual patients undergoing this procedure.
This study aims to develop parametric bone models to improve our understanding of the population variance of the calcaneus and medial cuneiform, which are the two bones that are typically reshaped to achieve surgical correction in adolescents with FFF-STA. The most common surgical approach is the Evans and Cotton osteotomies [
20] to convert a symptomatic FFF-STA to an asymptomatic FFF. This approach requires a reshaping of the calcaneus and medial cuneiform and lengthening of the calf muscles to release the gastrocnemius fascia [
21]. A lateral incision is made over the calcaneus to perform the lateral column lengthening osteotomy (Evans osteotomy) using a tri-cortical bone graft harvested from the iliac crest of the patient [
20]. Another incision is made over the medial cuneiform to perform an opening wedge medial cuneiform osteotomy (Cotton osteotomy) [
22,
23].
There are no studies published that have looked at the variation in the anatomical shape of these two specific bones, the calcaneus and the medial cuneiform, in skeletally mature adolescents with and without flatfoot deformity. Bone variation may be an important factor in the outcomes of flatfoot corrective surgery, and therefore, a study that can quantify and analyze this variation is an essential step in understanding its potential role in patient outcomes. The present paper outlines the methodology to achieve these goals within the context of this study. Principal Component Analysis (PCA) is a dimensionality reduction method commonly used in the medical imaging community to simplify the representation of shape variation. This approach will be applied to a dataset of calcanei and medial cuneiform bones to help improve understanding of variation in the bone morphology of these two bones across the population. The authors note that normal foot and flexible flatfoot demonstrate similar bone morphologies, where the key differences in the flatfoot condition lie in bone positioning under weight-bearing, soft tissue structures and pain experienced by the patient. Thus, using general population data will include both groups and will be analogous to analyzing a sub-set of flexible flatfeet.
2. Materials and Methods
The methodology applied in this study is based on the parametrization paper and MATLAB (v. 17, The MathWorks, Inc., Natick, MA, USA) script written by Pascoletti et al. (2021) [
24] which was used to analyze the morphological differences in the mandible and the proximal femur. The same general methodology is applied here on the calcaneus and the medial cuneiform of the foot.
2.1. Dataset
Initially, five computed tomography (CT) scans of feet from three subjects were used for a preliminary analysis [
25]. These were acquired from cadavers, assuming that they did not suffer from any foot-related condition and without specific regard for sex, age or weight. To expand the sample size, the New Mexico Decedent Image Database [
26] was used to collect the lower limb CT scans of 10 men and 10 women (40 calcanei and 40 medial cuneiform bones from right and left feet) that have been selected from a database of more than 15,000 subjects who died between 2010 and 2017. The dataset included relevant medical history such as cause of death, any chronic conditions, sex, age and weight. The subjects have been selected by looking at their medical history and by making sure that they did not have any previous pathologies or surgeries that could affect their foot bones’ shapes and that the cause of death did not influence the condition of the foot bones. After further scrutiny, the CT scans of three subjects were not included in the analysis due to prior lower limb surgeries which could have affected the bone morphology of the foot bones and were replaced by three other subjects without foot- and ankle-related issues.
Table 1 outlines the subjects that were finally analyzed, which totaled 23 subjects and 45 CT scans of the calcanei and medial cuneiform bones. The sample size used in this study is comparable to the number of bones analyzed in similar bone parametrization studies [
24,
27,
28].
2.2. Preparation of the Surface Meshes
The CT-scanned images were segmented using the open source software 3D Slicer (v. 5.0.3) with a threshold range of 250 HU (grayscale Hounsfield units) to maximum to create a surface mesh. The surface mesh was re-meshed using an isotropic explicit re-meshing tool in MeshLab (v. 2021.10, open source) to make a smoother homogeneous mesh. Amira (v. 2021.1, Thermo Fisher Scientific, Waltham, MA, USA) was used to place 18 corresponding landmarks on the calcanei and 6 corresponding landmarks on the medial cuneiform bones. Therefore, each bone mesh had the same landmark identification number located on the same anatomical position. The number of landmarks and their locations were defined by the co-author (BB), who is a specialist foot and ankle surgeon (
Figure 1,
Table 2 and
Table 3).
2.3. Creation of Iso-Topological Meshes
Transformations were carried out on each individual shape to scale, rotate and locate the models to one individual bone, which was used as a benchmark, so that they were considered as shapes ’only’, which is defined by D.G. Kendall [
29] as follows:
…all the geometrical information that remains when location, scale and rotational effects are filtered out from an object.
In this way, all shapes were comparable, and information related to their location, scale and rotation was eliminated. Generalized Procrustes Analysis (GPA) was applied, using the same procedure implemented by Pascoletti [
24], to the current dataset in order to find the average shape [
30].
In order to perform Principal Component Analysis (PCA) on shapes, it is necessary to have corresponding landmarks in every shape in the training database, to simplify the correspondence problem related to 3D shapes. Iso-topological meshes were created implementing a mesh morphing technique, making the PCA results extensive, where the (x, y, z) coordinates of all nodes of the surface meshes were considered as variables for the PCA (5310 nodes for the calcaneus and 1886 nodes for the medial cuneiform). Hence, all the surface meshes are required to have the same number of nodes, which are connected in the same way and where nodes are identically listed by position order. This resulted in all shapes in the set having identical meshes, with the coordinates of the nodes being the only difference between meshes, based on shape differences.
This can be achieved by many approaches, but here the Radial Basis Function (RBF) method was applied, which consisted of selecting one shape as a ’target mesh’ and using every other shape in the dataset as a ’standard mesh’ that morphed onto the ’target mesh’. The MATLAB script written by Pascoletti [
24] ran as follows:
The ‘standard mesh’—the mesh that will be morphed—was superimposed onto the ‘target mesh’, eliminating the location and rotational effects, by applying some transformations of rotation and translation, calculated from the landmarks.
The mesh was then morphed through the RBF method, and to reproduce a similar shape, all nodes in the standard mesh looked for the closest center of gravity of a triangle in the target mesh and were then projected perpendicularly.
A Laplacian smoothing was applied to the standard mesh to improve mesh quality, and then step number 2 was repeated. These two steps were then applied iteratively until the standard mesh fit the target mesh well, without seeing any big gaps that made the surface mesh look irregular.
Having all the shapes in the database aligned, an average shape was calculated using the Procrustes mean (also called the Frechét mean), which is found using Equation (
1) [
31]:
2.4. Principal Component Analysis (PCA)
With 45 iso-topological meshes of the calcanei and medial cuneiform bones completed, the PCA was applied to each database. This led to eigenvectors () that were the principal components (also called modes) explaining the shape variation and to corresponding eigenvalues () that were scalars corresponding to the importance (quantified by the percentage of the variance in the data explained) of the corresponding principal components. Eigenvectors were sorted starting with those with larger eigenvalues, which explain a larger percentage of shape variation.
Using the MATLAB script written by Pascoletti et al. [
24,
27], the number of principal components needed to be able to describe 90% of the shape variation of the database was found. Once the principal components were calculated and it was decided how many were needed, the average shape of the bone was then compared with the shape described by each PC, considering the principal components one by one in a way that a qualitative analysis could be made, linking a shape variation to each principal component. For example, in order to observe what variation PC1 described, Equation (
2) was applied, where
was the average shape,
was the mode and
was the corresponding weight of the first PC.
The script also provided a printout of the coordinates of the landmarks of the analyzed bones once they were all aligned in the same coordinate system. This step was important because the landmarks used were picked by the foot and ankle expert as anatomical landmarks that corresponded to exact points in the bones, so knowing their position allowed for the assessment of how those relevant points varied across the population. These were analyzed using MiniTab (v. 2021, MiniTab, LLC, State College, PA, USA), now applying PCA only on the landmarks’ coordinates.
3. Results
Figure 2 and
Figure 3 show two graphs of cumulative explained variance for the calcaneus and medial cuneiform datasets, respectively, against the number of principal components that are needed to explain said variance, produced by applying Principal Component Analysis (PCA) on all nodes of the surface meshes of the two datasets. To explain 90% of the total variance of the calcaneus dataset, 16 modes are sufficient, while for the medial cuneiform, 90% of the cumulative variance was explained by 12 modes, as shown in
Table 4 and
Table 5.
The first PC is always the most important eigenvector, explaining the highest variability, with the following PCs describing increasingly less variability within the dataset. For the calcaneus dataset, the first two PCs explain 16.9% and 12.3%, respectively, while for the medial cuneiform, they describe 35.2% and 17.1% of the variance, respectively. All the principal components are represented in
Figure 4 and
Figure 5 for the two bones, respectively, from the first one in the top left corner, to the last one in the bottom right corner. Each is shown with a colormap where the PC shape is superimposed by the average shape of the bone and the locations showing substantial deviation are highlighted in red. This qualitative analysis of the results of the PCA was made to better understand what each principal component expresses dimensionally. In the colormap of the first PC of the calcaneus dataset (
Figure 6), the area with most difference is at the top of the posterior surface of the bone, where the Achilles attachment is located (around landmark 18). Similarly, in the colormap of the second PC, the area with the most difference is the posterior talar articular surface (between landmarks 17 and 7), and then with the least difference there is the fibular trochlea (around landmark 11) and the sustentaculum tali (around landmark 14) [
32,
33], which are indicated in more detail by the red and yellow areas in
Figure 6.
In the colormap of the first PC for the medial cuneiform (
Figure 5), the area highlighted in red is in the medial part of the bone (around landmark 5), while in the colormap for the second PC, the areas in red are placed in the superior part of the bone (around landmarks 1 and 2) and in the inferior part close to the neighboring navicular bone (around landmark 4).
A quantitative analysis was also carried out, applying PCA only on the coordinates of the landmarks, where each landmark became three variables in the analysis, given by its x, y and z coordinates (54 variables for the calcaneus and 18 for the medial cuneiform).
Figure 7 and
Figure 8 show the score plots of each subject for the calcaneus and medial cuneiform, respectively. For the calcaneus dataset, there is a clear outlier, which is the black point in the lower left corner, which represents one of the first five CT scans with no additional data. It also shows that the female datapoints are clustered together, showing different variance in PC2 compared to male datapoints, also clustered together. For the medial cuneiform dataset, there were no clear outliers and the datapoints of the two sexes are scattered through the plot, showing no evident clusters or trends according to sex or body mass index (BMI) score.
Figure 9 and
Figure 10 show the loading plots for the two datasets, respectively, where the variables with the largest effect on the first two components are shown according to the length of the lines (magnitude). These are essentially a graphical summary of the correlation structure of the data. Numbers denote the landmark point, and letters denote the Cartesian coordinate (x, y, z) of that landmark. For the calcaneus, point 18x (landmark 18 in the x direction) has a large negative loading on the first component, indicating that the variation in 18x contributed the most to the 20.1% variation in the dataset. Similarly, other landmarks that had large effects include 11y and 11z with large positive loadings and 11x, 17z and 17y with negative loadings on the second component, so these variables contribute the most to the 12.3% shape variation. For the medial cuneiform, the loading plot shows that 5y and 6x have large positive loadings and 1y, 4y, 5x and 6y have negative loadings on the first component (accounting for much of the 31.0% variance), while 2x has a large positive loading and 1x, 2y, 4x and 4z have large negative loadings on the second component, contributing to the 17.0% variability in the shape. A large positive or negative loading (magnitude and direction) gives information on how a variable correlates with another, subject to the caveat this is always approximate, given by the percentage variance explained by a given PC. For example, two arrows going in the same direction are strongly positively correlated between each other (both positives and negatives), meaning that if one variable increases, the other increases, as well. For example, in
Table 6, 11y and 11z have almost identical angles to the x-axis (83.8° and 84.2°, respectively); therefore, these two variables are positively correlated. However, two arrows going in opposite directions (almost 180° between them) are strongly negatively correlated, meaning that if one increases, the other decreases. For example, in
Table 6, points 2z and 18x are almost completely opposite in direction, meaning that they are negatively correlated. If two arrows are perpendicular relative to each other, they are uncorrelated, such as points 1y and 4z in
Table 7. The magnitude and direction of the arrow explain the importance of that variable to either the first or second component. For example, 18x (
Table 6) is the most influential landmark for the first component of the calcaneus because of its great magnitude and its angle being almost horizontal on the axis. In comparison, 11x (
Table 6) is the most influential landmark for the second component, as shown in the plot, based on the vertical position and magnitude.
4. Discussion
Flatfoot surgery is complex and has many variable factors that influence surgical outcomes. Therefore, this study investigates the extent of shape variation in the population as one of these key factors to be understood. This research study examines the bone morphology differences in the calcaneus and medial cuneiform, two bones that are re-shaped during flatfoot surgery to convert a painful flexible flatfoot with short Achilles tendon (FFF-STA) to an asymptomatic flexible flatfoot (FFF). With the help of 18 landmarks for the calcanei and 6 landmarks for the medial cuneiform bones in 45 foot scans, the meshes were altered to be iso-topological meshes, all with the same rotation, location and size. The graphs from the Principal Component Analysis (PCA) show the cumulative explained variance against the principal components. Ninety percent of the total variability can be explained with 16 PCs and 12 PCs for the calcaneus and medial cuneiform datasets, respectively (
Figure 2 and
Figure 3). This can be explained by the geometry of the bone itself. The calcaneus is the largest bone of the foot and has a complex shape that allows it to transmit the weight of the body to the back of the foot (heel) and to the front of the foot through the cuboid bone [
34]. Depending on the forces exerted on the bone, it will change its shape and composition, increasing or lowering its amount of cortical and cancellous bone. The medial cuneiform is a wedge-shaped bone located in the midfoot, acting like the keystone of the medial longitudinal arch of the foot. It has a more regular shape and it is smaller compared to the calcaneus. For these reasons, it seems reasonable to determine that four additional PCs are needed to explain a bigger and hence more complex bone [
35] for the calcaneus.
The qualitative colormaps (
Figure 4 and
Figure 5) and the quantitative loading plots (
Figure 9 and
Figure 10) present the same data in a way that effectively allows us to observe the anatomical location (colormaps) and quantify the direction and deviation (loading plots) of key shape variations. For example, in the colormap of PC1 for the calcaneus, the area highlighted in red is identified on the loading plot as landmark 18 in the x direction. The same happens for PC2. The location of the osteotomy in the calcaneus is lateral, and although PC2 also describes the area around the fibular trochlea, which is affected by the Evans osteotomy, PC2 mainly explains the area of the posterior talar articular surface, and overall, it describes only 12.3% of the total variance of the whole calcaneus database. The following PCs explain even less than 12.3% each, which is not a considerable contribution to the shape variation. For the medial cuneiform, the most notable observations in
Figure 5 are the areas highlighted in red in the colormap for PC1 and PC2, which are confirmed by the loading plot (
Figure 10). The location of the Cottons osteotomy is in the superior part of the medial cuneiform, and although PC2 describes only the vertices of the superior part of the bone, it is possible that if we had placed a landmark in the middle of these two points, it would have given us more information about this location and hence inform us about the population variability of this osteotomy site.
This research study informed us that there is not as considerable a difference in bone shape in the calcaneus and medial cuneiform as we thought there would be. PC2 for both datasets could potentially describe the osteotomy sites for both bones. Further studies focusing only on the osteotomy sites in both bones would possibly reveal more information regarding their implications in the surgical planning and outcomes. Another potential study with this dataset would be to investigate how many landmarks are sufficient to identify the shape of the calcaneus, particularly using landmarks that are only visible from the sinus-tarsi approach. The ability to characterize the shape and orientation of a bone from a small surgical cut would potentially allow the possibility of using robotic-guided surgery in the foot and ankle.
Learning about shape patterns on bones could have implications in surgery, as the surgeon will have a better understanding of where shape variation is most critical, and therefore they may influence where a surgeon may cut the bone. For example, this could help planned surgeries to avoid cuts in areas of high shape variation, which could reduce post-surgery complications and uncertainty and improve surgical reliability.
What is most interesting about these shape patterns in the loading plots is that the direction and magnitude of each landmark variable give information about the interdependence between landmarks. This could show coupled or clustered shape changes where one landmark proportionately or inversely proportionally changes another. We have noted from the results that the osteotomy site does not appear to show high variability in the population; this may be an advantage for surgeons as the surgical cuts can be more consistently planned, knowing which bone landmarks showing low population variability could be used as reference points. A research study that focuses on this simplification of the dataset could be useful to reduce the number of landmarks needed to fully identify a bone shape. Due to the complexity of the foot bone shapes, multiple bones and joints in the foot, there are still logistical challenges in implementing robotic surgery in the foot and ankle. However, given the complexity in these surgeries and the many factors affecting surgical outcomes, this study may be a step towards quantifying and understanding bone shape variabilities of the foot. This work could lead to further studies characterizing interdependent anatomical patterns of foot bones, which could be a game-changer for surgical robotics in the foot and ankle to improve outcomes and surgical consistency.
While the score plot for the medial cuneiform (
Figure 8) dataset gives no additional information, the score plot for the calcaneus (
Figure 7) shows female subjects and male subjects each clustered together, indicating more similarities among them. An anthropology study could further this finding by studying the bone morphology differences of the calcaneus with a gender study and how this could affect the osteotomy location for lateral column lengthening osteotomy or similar foot surgeries.
An important limitation that needs to be mentioned is that the foot and ankle specialist did not position each landmark on the bones meshes himself. The positions of these landmarks were taught by the surgeon (BB) and applied by the main author (YC). Although this means the interpretation of where the landmarks were was decided by a non-clinical researcher, the surgeon identified the number of landmarks and where they should be marked on 2D images of the bones with clinical reasoning as to why he was placing a landmark in that position. This was then implemented by a single person (YC), eliminating any inter-observer variation. Furthermore, analyzing a dataset of only 45 bones to represent a population is another limitation of this research paper. Although Pascoletti et al. [
24,
27] used 40 and 50 bones, respectively, in their parametrization of the mandible, this is not an extensive dataset in order to reach confident results to be applicable in a surgery setting. Another study that analyzes a bigger pool of subjects’ bones may achieve results that strengthen the results presented here. In order to perform a successful flatfoot surgery, the surgeon needs to consider other patient variabilities that go beyond the shape differences of the calcaneus and medial cuneiform only and would need to include a dataset of clinically diagnosed painful flexible flatfoot with short Achilles tendon (FFF-STA). Despite these limitations, this project was focusing on the bones that are directly cut and reshaped in flatfoot surgery, and therefore, those two bones’ shapes are the most critical variables to look at as a starting point. A further study could look at specific surgical approaches and how critical the position and cutting approach is, considering bone shape variability. Understanding where bone variability is found greatest and least will help with surgical planning in the future.