The objective of this research is to model the geometric variability of the glenoid of the scapula. The glenoid is the “socket” component of the “ball and socket” connection of the shoulder joint. The model must capture the observed variability with sufficient resolution such that it informs both operative and design decisions. Creating the model required the application of existing mathematical and statistical modeling approaches, including geometric fitting, radial basis functions (RBFs), and principal component analysis (PCA). The landmark identification process represented the glenoid in a new manner. This work was validated against existing approaches and computed tomography (CT) scans from 42 patients. Information on the range of shoulder geometries can assist with preoperative planning as well as implant design for total shoulder arthroplasty (TSA). PCA was used to quantify the variability of shape across landmarks used to represent the glenoid shape. These landmark locations could be used to generate full surface meshes of existing glenoids or new glenoid models synthesized by changing principal components (PC). The process of creation of these shoulder geometries may be useful for the study of other joints. The models created will help surgeons and engineers to understand the effects of osteoarthritis on bone geometry, as well as the range of variability present in healthy shoulders.

## Introduction

Total shoulder arthroplasty (TSA) is performed in patients who have lost mobility in their shoulder due to arthritic wear on their *glenoid*. The glenoid is the socket component of the shoulder joint. An implant is inserted in the glenoid vault to support the new ball joint implanted on the humerus. Complications, such as improper sizing and alignment of this implant, can cause poor surgical outcomes. These include failure to restore full mobility to the shoulder, losing mobility over time due to loosening of the implant, and shoulder pain. The alignment of the glenoid component is complex due to the variation of orientation among healthy patients and the wear on the original glenoid face caused by osteoarthritis. In order to predict the effects of osteoarthritis and create implants that can restore the shoulder, it is important to fully understand the variability across healthy patients.

This research was conducted in collaboration with a physician who specializes in shoulder and elbow surgery. Computed tomography (CT) scans of both healthy and osteoarthritic shoulders were provided for this research. The CT scans of healthy and osteoarthritic patients were used to create a statistical method of estimating glenoid and glenoid vault geometry.

The objective of this research is to model intact and osteoarthritic glenoids. These models must be robust to the observed variability of the glenoid, with sufficient resolution, such that this information can be used to improve surgical methods and implant design.

This paper is structured as follows: Section 2 of this paper presents the physiology of the glenoid, information on TSA, shape analyses used in anthropometry, and the need for shape analysis of the glenoid. Section 3 covers the sphere fitting analysis used to quantify the geometry of the shoulder, the meshless representation used to represent glenoid geometries, the principal component analysis (PCA) conducted, and the synthesis of new shoulder geometries. Section 4 presents the validation of the sphere fitting, the trends founds in the variability analysis, and the synthesized shoulders created through varying the principal components (PC). In Sec. 5, the implications and future work that can be conducted with the results of this study are discussed. Finally, Sec. 6 summarizes the results of this paper and the next steps for future study.

## Background

The shoulder joint consists of the humerus and scapula and is held together by a group of muscles and tendons known as the rotator cuff. The scapula (Fig. 1) has a complex shape. This research focuses on the glenoid and *vault*, the thin section behind the glenoid that supports the joint. The scapula, surrounded by the neighboring bones, can be seen in Fig. 2. The humerus and scapula are labeled. The directions used to describe orientation within the body are shown in this figure as well. *Superior* and *inferior* refer to the direction from the top of the head to the bottom of the feet, while anterior and posterior refer to direction that extends from the front to the back of the body. *Medial* refers to the direction starting at the left or right of the body toward the middle of the body, and *lateral* refers to the direction originating at the middle of the body extending out to the right or left.

The shoulder joint is not a perfect ball and socket joint, as there is some translation of the humeral head with respect to the glenoid [1,2]. Due to the complex motion of this joint, it is considered to be slightly unstable [3]. The orientation of the glenoid with respect to the scapula can be described by the glenoid *version* and glenoid *inclination*. The version represents the tilt of the glenoid around the superior–inferior axis, while inclination represents the tilt around the posterior–anterior axis.

### Osteoarthritis and Total Shoulder Arthroplasty.

Osteoarthritis is the loss of cartilage at a joint. Glenohumeral osteoarthritis (GHOA) may result in the deformation of the glenoid, resulting in a loss of mobility and pain at the joint. Posterior wear induced on the glenoids in cadavers resulted in increased translation of the humeral head [4]. TSA, a type of shoulder implant, is used as a way to regain mobility of the shoulder and lessen the pain. During TSA, a front incision is made on the patient's shoulder. The humerus is dislocated, and the top of the humeral head is sawed off. A metal prosthesis is cemented into the humeral shaft, and a hemispherical implant is screwed onto the shaft.

One of the complications of TSA is the glenoid component of the implant loosening over time. A long term study on 113 patients with the Neer prosthesis found implant survival was 94% after 10 years and 87% after 15 years [5]. However, 44% displayed evidence of glenoid loosening, which was associated with pain [6]. By preventing the glenoid from loosening, it may be possible to improve surgical outcomes.

The modern glenoid component is made from polyethylene and is cemented into the glenoid. In an effort to prevent the loosening that can occur with a cemented glenoid, metal-backed implants were created to encourage bone growth into the implant. A polyethylene liner is placed between the glenoid and the humerus implants to prevent metal on metal contact. The polyethylene liner was found to wear faster than the all-polyethylene implant [7,8]. While the metal back itself does not loosen, the liner can disassociate from the implant [9].

### Glenoid Studies.

The version (superior–inferior tilt) of the glenoid is an important factor in TSA, and many studies have worked on identifying trends across the version of healthy and osteoarthritic shoulders. The forces on the glenoid are directly affected by the version of the glenoid component [10]. CT scans of scapulae provide accurate three-dimensional (3D) models that can be used for bone measurement [11]. In a two-dimensional (2D) scan, variables such as patient posture can influence the results [12]. The traditional method of measuring version at the rim in 2D view changes at different glenoid heights, making the exact angle of version difficult to estimate. Lewis and Armstrong [13] used sphere fitting as a method of predicting glenoid version and found that there was a range of version and inclination values across intact glenoids.

To compensate for the retroverted surface observed in glenoids with osteoarthritis, a posterior step is suggested to turn the implant back to a healthy angle [14]. This suggestion is complicated due to the range of version angles that can be considered to be healthy [13]. Preliminary work is being done on using 3D planning tools to place glenoid implants on cadavers [15]. This leads to the possibility of 3D printed implants specific to the patient's level of osteoarthritis. Simulations are conducted to improve upon glenoid implant placement in 3D models, but fail to take into account the natural range of healthy version angles [16]. There is, therefore, a need to isolate the changes in shape due to osteoarthritis from the natural variability that occurs across intact shoulders.

### Studying Variability.

To design for a population with a large amount of variability, it is important to have an accurate representation of the population's anthropometry. However, due to the time and cost involved in obtaining data, it can be difficult to obtain enough information to represent a population. Data can be synthesized based on the available information obtained in a sample to represent the variability in the population.

Simple methods such as proportionality constants (scaling to stature) fail to take into account variability in the shape of the individual. Another approach is to formulate linear regression models relating the shape and landmarks to easily obtainable *predictors* such as stature and body mass index, or glenoid height and width. These models can be based on the detailed information from one study and may then be extrapolated to a desired target population for which the predictors are known [17]. Approaches using PCA have been shown to be very effective as well [18,19]. Shape analysis across varying populations has been completed in several fields, including the study of total knee arthroplasty [20,21].

Radial basis functions (RBFs) are functions whose values are solely based on the distance from a point in space to the function's center. A combination of radial basis functions can be used to interpret multivariate scattered data by approximating a mapping function. These functions can depend on many variables and are defined by many data. Unlike various other methods, the data do not need to be in a grid and can be scattered in their domain.

Radial basis functions are used to represent 3D surfaces through a system of linear equations [22]. These functions can be used to create surfaces to represent point clouds and can be used to fill holes in meshes [23]. By using them, and optimizing for speed, it is possible to cut down on computation time when compared with mesh-based methods [24]. Instead of using RBFs to represent a surface, they can also be used to deform known surfaces [25]. For example, Bennink et al. [26] created a 3D warping model to track the changes of brain shape during neurosurgery.

The combination of RBFs and PCA has been used to study changes in skull shape of infants and torso shape of adults to create better impact safety predictions [27]. Radial basis functions are useful to this endeavor, since they can be used to visualize predicted changes in shape or fit a grid-based data function for computational measures that require organized data.

## Methodology

It is important that the glenoid models are of sufficient resolution to be accurate. This involves processing the provided CT scans in order to produce meshes for each shoulder. The orientation of the glenoids relative to their scapulae was found by sphere fitting the original meshes. In order to study trends within the geometry, it was crucial to have a consistent set of landmarks placed for each glenoid. These landmarks were used to create new meshes through the process of RBFs, which could be sphere fit in order to validate the sufficient accuracy of the landmarks. The variation between landmark locations was identified using PCA. By changing the PC scores, it was possible to synthesize a new set of shoulder geometries. A visualization of this process is shown in Fig. 3.

### Quantifying Geometry.

The orientation of the glenoid with respect to the scapula can be described through its version and inclination. Predicted version of a glenoid changed at various heights when measuring in a 2D view. Lewis and Armstrong [13] introduced sphere fitting as a method to identify version consistently across the whole surface of the face. When the concave section of the glenoid face is isolated, a sphere can be fit to the resulting curved surface. The position of the sphere center determined the version and inclination of the glenoid face, along with the curvature of the glenoid face. A visualization of a sphere fit to the face of the glenoid is shown in Fig. 4. Two methods of sphere fitting were conducted: using a sphere with a fixed radius of 30 mm to determine the version and inclination alone, and using a sphere with a variable radius to determine version, inclination, and curvature of the glenoid face. The 30 mm standard was adapted from the Lewis and Armstrong study.

### Meshless Representation.

To restore full mobility to the shoulder, the hypothesis is that the implant should have the orientation of a healthy glenoid. The sphere fitting and vault fitting tactics are useful to identify trends in glenoid orientation. However, they do not represent all of the variability present in the shape of the glenoid. In order to fully model the effects of GHOA, information beyond univariate measures are needed.

Comparing the 3D meshes directly can be complicated as the meshes do not have consistent representation; the vertex index for one glenoid mesh will not represent that same point on another glenoid. In order to quantify the variability of shape and orientation, the shoulders need to be represented with a consistent and repeatable method. Landmarks were placed on the meshes representing healthy and osteoarthritic shoulders, such that their locations relative to the glenoid would be constant across the shoulders. RBFs allow for a mesh to be morphed according to changes in key landmark locations [22].

#### Automatic Landmark Setup.

Landmark identification needs to be consistent and reproducible across shoulders, with minimal error. These landmarks must be placed such that they could be used to identify and reproduce the glenoid geometry and must be unaffected by user error. Identifying landmarks on the glenoid is difficult due to the lack of bony protrusions. In order to limit the amount of error in the placement of the landmarks, three main regions for landmark locations were identified: the face, the rim, and the vault of the glenoid. These landmarks were projected onto the surface meshes of the glenoid using code created in matlab. The number of landmark locations was determined by limiting the amount of error caused when recreating shoulder meshes using only the landmark locations.

#### Face Landmarks.

*K*) was set to determine the number of landmarks (

*N*

_{LM}) placed on the face. The value of

*K*was adjusted and found through the sphere-fit validation. A total of 2

*K*+ 1 rows of landmarks were placed at set heights across the inferior to superior length of the glenoid, where from the inferior-most row the number of landmarks increased by two until the center row, which had 2

*K*+1 landmarks. The number of landmarks in rows superior to the center row decreased by two such that the following pattern was obtained:

These landmarks were used to recreate the geometry of the glenoid face using radial basis functions. If the geometry of a given face was well represented through RBFs, then the RBFs could also be used to represent new landmark locations of synthesized glenoid faces.

#### Rim Landmarks.

The rim of the glenoid face is the protrusion around the concave face of the glenoid. The glenoid was oriented to the plane that fit the glenoid face, and the vertices were expressed within a cylindrical coordinate system, using radius (*r*), angle (*θ*), and height (*z*). For each value between a set range, the vertices with the maximum height were collected. A curve was fit to these points. At constant *θ* values of this curve, the length of the radius *r* was found, and the resulting point was projected onto the face.

#### Vault Landmarks.

The vault landmarks were placed by projecting points onto the neck for the posterior, anterior, and inferior sides of the vault. The points were set at various depths that were dependent on the medial and lateral-most points of the extrusion of spine of the scapula. The vault landmarks were useful in determining the implant depth as well as the erosion to the posterior side of the glenoid. In order to place the points on the anterior side of the vault, a minimum and maximum depth compared with the rim was set, and landmarks were placed at these depths. The distance was set as 25% of the distance from the glenoid center to the point at which the medial border meets the scapular spine. This value was selected to provide sufficient information on the vault angle. The depth was set relative to the size of the total scapula, due to prediction that the base of the acromion would not change with GHOA. The distance from the rim of the glenoid to the lateral base of the spine and acromion is noticeably shorter in osteoarthritic shoulders. Therefore, the depth of the neck landmarks on the posterior side of the glenoid was determined by using the lateral protrusion of the scapular spine. The inferior-most point of the glenoid vault was selected at several different depths to match the depths of the posterior points.

#### Mesh Creation.

*baseline mesh*can be transformed to represent these new locations by creating a displacement field using its matching landmark. Using the baseline landmark locations ($p1$…$pN$) and the differences ($F1$…$FN$) between the landmark locations, radial basis functions can be applied in order to create a displacement field. Any point $p$ on the baseline mesh will have a displacement $F$. The 3D displacement vector can be defined as follows:

*σ*($p\u2212pi$) or

*σ*($r$), is a function whose value is based on

*r*, the distance a given point is from $pi$. In this research, the RBFs were set to a thin plate spline as written in the following equation:

**B**can be created for the RBF between each set landmark and the new landmark location

**Q**

_{rbf}is created

*F*,

_{x}*F*, and

_{y}*F*can each be represented by the following matrices:

_{z}By setting **F** to the known displacement values to displace the points from the original glenoid to the *comparison glenoid*, it is possible to find the weights for the RBFs ($ax,\u2009ay,\u2009az$) and the polynomial constants ($cx,\u2009cy,\u2009cz$). Inputting the original vertices through Eq. (3) with the known weights results in the displacements that can be used to transform the original glenoid to be represented using the vertices of the new glenoid.

#### Landmark Verification.

In order to validate that the landmarks alone could be used to produce accurate meshes for predicted shoulders, RBFs were used with the same base glenoid across the landmark locations for each shoulder. If the error of sphere fitting of the generated meshes averaged less than 5%, the landmark placements were deemed sufficient.

A sphere-fit calculation was conducted to see how accurately the transformed glenoid represented the comparison glenoid mesh. A flow chart of the validation process is shown in Fig. 5. Once a suitable glenoid was identified, its mesh was set as the baseline to represent the new shoulders.

### Quantifying Variability.

Principal component analysis was used in order to quantify the variability of the landmarks placed to represent the shape of the glenoids. PCA reorients data into a reference frame where each coordinate is set such that there is no correlation between components. These primary sources of variability were found and used in mathematical models to generate glenoid geometries. In a multidimensional problem, such as the landmark locations placed on the glenoid, transforming the landmarks into the uncorrelated space allows new locations to be generated. Generating a new landmark location directly in the original reference frame would be nearly impossible, due to the correlations that exist between landmarks. It is possible to generate new landmarks in the PC space; a linear combination of the PCs with variable PC scores that exist within the range of identified data will result in feasible new geometries.

A PCA transforms data from a correlated space into a reference frame where the scores for each principal component are uncorrelated. This allows new data to be synthesized easily. In order to create new members of a multivariable population, PCA can be used in the endeavor of creating *eigenfaces*. An eigenface is an eigenvector that represents a matrix of data. If the relevant data concerning a member of a population are represented in matrix form, the data must be reorganized into vertices. The data are organized such that each column represents a member of the total population, and each row represents a different point of information. For this research, the population consisted of the processed glenoids, and each glenoid was represented by its landmark locations. The PCA was conducted with these vertices which represent the matrix data to identify the eigenfaces.

The eigenface PCA was performed in the manner described in Refs. [28] and [29]. Each of the (*n*) three-dimensional vertices was stacked in the same column, and each separate example of the population (*k*) is listed in a separate row of height 3*n*. The *n* eigenfaces can be calculated by finding the eigenvalues of the *n* × *n* correlation matrix. However, only *k* − 1 eigenfaces will have nonzero corresponding eigenvalues.

Landmarks on the glenoid were placed in order to represent the shape of the glenoid, in an order that was repeated across the glenoids. In order to quantify the variability of the shape changes that occur with GHOA, PCA was conducted using these landmarks. Each set of landmarks was scaled and rotated using Procrustes analysis, such that only the differences in shape were considered.

### Synthesizing New Glenoids.

*σ*) of that PC score and setting all other PC scores to zero. Two standard deviations identify the maximum set point to provide an extreme but realistic view of the variation

The new stacked landmark data *s _{i}* created with Eq. (9) can be reorganized into an $n\xd73$ matrix of landmark data. To create new meshes to represent the changes in landmark data, the baseline mesh is centered and scaled in the same manner as its landmarks. This mesh was then morphed to represent the new landmark locations using RBFs. Finally, a sphere fit is conducted on the new mesh to see the effects of changing the PC.

## Results

This study was conducted using CT scans of intact and osteoarthritic scapulae provided by Dr. April Armstrong at Penn State Hershey Medical. Using this procedure, 42 scapulae (34 intact and eight osteoarthritic) were processed. In this research, 113 landmarks were projected on the isolated glenoid faces, 20 on the rim, and 174 on the vault. The validation of these landmarks is displayed in Table 1. The results of the PCA on the landmark locations are listed in Table 2. A visualization of the average glenoid is shown in Fig. 6, and the transformed glenoids are shown in Fig. 7. The sphere-fit results across synthetic shoulders are tabulated in Table 3 and shown visually in Fig. 8.

PC | Variability (%) | Total variability (%) |
---|---|---|

1 | 30.3 | 30.3 |

2 | 24.5 | 54.8 |

3 | 7.2 | 62.0 |

4 | 5.9 | 68.9 |

5 | 4.9 | 72.8 |

6 | 3.7 | 76.6 |

7 | 3.5 | 80.1 |

8 | 2.9 | 83.0 |

9 | 2.1 | 85.1 |

10 | 1.9 | 87.1 |

PC | Variability (%) | Total variability (%) |
---|---|---|

1 | 30.3 | 30.3 |

2 | 24.5 | 54.8 |

3 | 7.2 | 62.0 |

4 | 5.9 | 68.9 |

5 | 4.9 | 72.8 |

6 | 3.7 | 76.6 |

7 | 3.5 | 80.1 |

8 | 2.9 | 83.0 |

9 | 2.1 | 85.1 |

10 | 1.9 | 87.1 |

Adjustable | Sphere | |||
---|---|---|---|---|

($deg$) | Inc. ($deg$) | Radius (mm) | ||

Varying PC 1 and 2 | Min | −24.5 | −4.6 | 31.4 |

Max | 6.6 | 2.6 | 33.9 | |

Varying PC 3 and 4 | Min | −16.3 | −3.1 | 30.1 |

Max | −1.8 | 1.0 | 36.4 |

Adjustable | Sphere | |||
---|---|---|---|---|

($deg$) | Inc. ($deg$) | Radius (mm) | ||

Varying PC 1 and 2 | Min | −24.5 | −4.6 | 31.4 |

Max | 6.6 | 2.6 | 33.9 | |

Varying PC 3 and 4 | Min | −16.3 | −3.1 | 30.1 |

Max | −1.8 | 1.0 | 36.4 |

### Validation of Landmark Locations.

The results of using the radial basis function transformed glenoids to predict sphere fit are listed in Table 3. Each glenoid was set as the baseline glenoid, and its mesh was displaced to the locations of the remaining comparison using RBFs. A sphere was then fit to the new surface. The combined average error for any given glenoid being used as the baseline mesh was 1.1 mm, which is just over a 3% error. Out of all the glenoids, those of patient 9 and patient 72 had the least amount of error when used as the baseline mesh. Patient 72's glenoid had an average error for sphere radius of under 1 mm, which was less than 2% of the total value, and the version angles were all within 1 deg, with the mean being within 0.22 deg of the original predicted value. This was determined to be satisfactory, and patient 72 was set as the baseline mesh.

### Variability Analysis.

To identify the total variability of the glenoid shape, PCA was conducted using the landmarks placed on the glenoids' faces and vaults. The effects of altering scores for the top four PCs were explored in order to show the primary variation within the data. The sources of variation were tested to determine if there were any trends that correlated with GHOA.

By using linear combinations of the first ten PCs with the mean of the original data, 87% of the variability of the landmark locations can be represented. The variability expressed by the top ten PCs is displayed in Table 1. The mean values for the landmarks were transformed into a mesh using RBF and are displayed in Fig. 6.

### Synthesized Shoulders.

PC 1 and PC 2 represented over 60% of the variation of the data and had different mean values between intact and GHOA glenoids. The values of PC 1 and PC 2 were systematically varied in a parametric study. The other PC values were kept at the mean value of zero. By varying these two values, a set of shoulders could be created to represent potential variations in nature. This process was repeated with changing values of PC 3 and PC 4 scores, with the PC 1 and 2 scores set to represent the mean glenoid. The PCs value ranged from a minimum value at two standard deviations below the mean to a maximum value of two standard deviations above the mean. The summary of the effects of changing the PC scores is listed in Table 2. A visualization of changing the first two PC scores is displayed in Fig. 8. The effect of changing the scores of the first two PCs on the sphere fitting results is displayed in Fig. 7.

## Discussion

This project was conducted to identify the variability between glenoids for the purpose of improving TSA outcomes. For this research, CT scans of intact and osteoarthritic scapulae were processed and studied in order to quantify the variability present in glenoids. The landmark locations were cast to represent the topology as consistently as possible across the glenoid geometries. If the locations were exactly consistent across topologies, the sphere-fit validation testing would have been without any error. However, since the error was small, the landmarks represent the topology of the glenoids satisfactorily.

Using this methodology of landmark representation, it may be possible to predict the effects of osteoarthritis in unhealthy shoulders and restore glenoids to their original shape. By increasing the number of osteoarthritic shoulders in the study, it is possible to identify the PC scores that correlate with osteoarthritis. Preliminary differences across the osteoarthritic versus intact populations can be seen in the top two PC scores, but not the other top ten PC scores. This indicates that the first two PCs may represent the shape change that is caused by osteoarthritis, while the other PCs represent the variability present across intact shoulders. If this trend continues across multiple shoulders, then it may be possible to restore osteoarthritic shoulders to their intact geometry by adjusting the scores of the first two PCs to that of the average intact shoulder. The predicted intact geometry would be used as a reference for how an implant should be oriented.

With enough evidence that this model is useful, these data can be used to design implants for glenoids in the future. Predicting the original healthy angle could improve prosthetic design. It may be possible to use these results in collaboration with current step prosthetic designs. At this point, it is possible to create a group of models with varying geometry. The next step would be to adjust this method to include sizing changes.

These osteoarthritic models can be useful in the endeavor to design appropriate implant families. By having a prediction of the healthy geometry, the angle of restoration and the size of the glenoid can be used to select implant step levels. Implant sizing would need to reflect the variability in sizes of the glenoid face, and the different restoration angles needed to support the humeral head. This method of landmark setting, shape morphing, and population generation may be useful for implant design for any bone, as well as other design for human variability endeavors.

## Conclusion

The purpose of this research was to create a model that can generate intact and osteoarthritic shoulders. Due to the lack of distinct features on the glenoid, it can be difficult to represent the complexity of its shape. A set of landmark locations were selected to represent the topology of the glenoid. These landmark locations were transformed into a PC frame of reference to find the variability that spans across glenoid geometries in intact shoulders. By manipulating the PC scores, a suite of shoulder geometries was created within the observed variability from existing shoulders. These geometries have the potential to be used in the endeavor of designing implants for total shoulder arthroplasty.

## Acknowledgment

Thank you to Dr. April Armstrong from Hershey for providing the CT scans and her advice.

Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessary reflect the views of the National Science Foundation.

## Funding Data

National Science Foundation (Award No. 0729386).