Construction of 3-D Humeral Head Statistical Shape Model in CT Images

Featured Application

The research work is aimed to provide a statistical shape model of the humeral head towards the treatment of glenohumeral osteoarthritis.

Abstract

Replacing the humeral head with an artificial one via surgery is one of the options to treat glenohumeral osteoarthritis. Thus, designing the artificial humeral head is an important step to alter clinical outcomes. In order to design the artificial humeral head, the individual variety of the humeral heads should be investigated. The statistical shape model (SSM) has been attracting considerable attention to grasp 3-D shape variety; however, no method to derive the SSM of humeral heads has been studied. This paper proposes a method to construct an SSM of humeral heads based on the anatomical landmarks in shoulder computed tomography (CT) images. The proposed method consists of three steps: humeral head extraction, position and pose alignment, and finally, principle component analysis. The method was applied to 22 male subjects with leave-one-out cross validation. The proposed method obtained an average Dice coefficient of 0.92 to represent the individual shape using the constructed SSM. According to shape analysis of the humeral head, we found that the thickness of the humeral head was associated with individual characteristics of the humeral head. Therefore, it can be said that this study can provide patient-specific design of an artificial humeral head.

1. Introduction

Osteoarthritis has a negative impact on daily activities for millions of elderly people around the world [1,2]. In relation to the aging of articular cartilage or rotator cuff dysfunction, glenohumeral joint osteoarthritis prohibits smooth motion between the humeral head and glenoid [3]. Previous epidemiology studies show that 32.8% of people with over sixty years old are affected by glenohumeral joint osteoarthritis and debilitated joint function due to the age-related change [1]. The breakdown of cartilage causes the humeral head and glenoid to rub and grind against each other [3]. It induces acute pain and loss of function, i.e., mobility of the arm [2]. To compensate for the cartilage breakdown, excess cells are produced in the humeral head in the form of bone spurs. This leads to even more friction, pain, and limiting of function. Thus, one of the most popular treatments of this disease is to surgically remove the humeral head and replace it with an artificial one [3]. The design of an artificial humeral head is likely to alter the surgical outcome for glenohumeral joint osteoarthritis. Research studies showed that artificial humeral head design with some kind approximation (e.g., spherical or elliptical) did not bring an adequate outcome [4]. It also limited the range of motion (ROM) of the shoulder [5]. Some studies showed that increasing the similarity between the artificial humeral head and the actual one could increase the ROM of the shoulder [4,6]. Therefore, the focus of this research was to discuss a method that can provide patient-specific design of artificial humeral heads.

From the previous study in humeral head morphology, the humeral head was composed of different curvatures [7,8]. In addition, Humphrey et al. [9] clarified that use of elliptical heads resulted in improved replication of the normal humeral head shape in comparison with use of spherical prosthetic heads. Furthermore, the cadaveric study by Jun et al. [4] revealed that the custom, non-spherical prosthetic head more accurately replicated the head shape, rotational range of motion, and glenohumeral joint kinematics than the commercially available, spherical prosthetic heads when compared with the native humeral head.

However, humeral head size, including interindividual variation, has not been well investigated. The statistical shape model (SSM) [7] is one of the methods that expresses the interindividual variation of the shape statistically, and is basically used as prior knowledge of various algorithms. SSMs were used in number of medical applications, especially related to bones and joints. To assess the potential risk factors associated with bone fracture, an SSM was used by Gregory et al. [10] in their study. They concluded that SSMs might be useful in identifying the individuals at risk of possible femoral fracture. To aid detecting vertebral fractures, SSMs of vertebrae were used by Roberts et al. in [11]. SSM was used successfully for early identification of developing radiographic osteoarthritis of hips by Gregory et al. [12]. Apart from that, SSMs were used to automate implant design and screw placement in shoulder arthroplasty [13], surgical planning [14], medical image segmentation [7], and implant design [15,16].

SSMs avoid arbitrary assumptions by trying to capture the actual patterns of variability in a set of objects of same class. It is usually constructed by applying principal component analysis (PCA) [17] to feature vectors of multiple individual shapes [18]. PCA is used in other statistical modeling as well [19]. Feature vectors can be classified into two categories: point-based or distance-based feature vectors [7]. The point-based feature vector is calculated from 2-D or 3-D coordinate values of feature points. The feature points between individuals should be associated. Thus, it requires a feature point extraction method, which limits the application. The distance-based method segments the target organ first, and then makes distance map from the extracted organ contour. The distance value of all pixels in 2-D or voxels in 3-D is utilized as the feature vector. The remaining work is to align the position and orientation of the target organ, and the alignment method is different for each organ and for each purpose. As of now, there are no studies which investigate the procedure of constructing the SSM of a humeral head.

This study aims to propose a method of constructing the SSM of a humeral head towards the design of patient-specific artificial humeral heads. The proposed method consists of three steps: segmentation of humerus, humeral head extraction, pose alignment of humeral head, and finally, construction of the SSM by applying PCA. This study particularly focuses on the shape of the humeral head, whereas our previous study [5] focused on the shape of the whole humerus. Furthermore, this experiment demonstrates the effectiveness of SSMs by studying the individual shape variability using the SSM. The rest of the paper is structured as follows: Section 2 discusses the images used in this study and Section 3 describes the proposed methods. Section 4 presents the experimental results. Finally, the paper is concluded by summarizing the experimental results and the achievements of the research, as well as future research directions.

2. Subjects and Materials

A total of 22 male subjects with age range from 18 to 79 years and mean ± standard deviation (SD) age of 39.0 ± 20.4 years old were considered to construct the humeral head SSM. The subjects had no significant disease in their shoulder joints. The local ethics committee of Nobuhara Hospital, Tatsuno, Japan approved this study and all the subjects provided written informed consent according to the approved procedure.

Axial computed tomography (CT) images of the left shoulder were acquired for each subject to cover the whole shoulder joint using a CT scanner (120 kVp, 250 mA; ECLOS, Hitachi Medical Corporation, Japan) at the Nobuhara Hospital from 21 February 2013 to 20 November 2018. The imaging technique that creates two-dimensional cross sectional images from three-dimensional body structures is known as computed tomography or CT [20]. The combination of volume and pixel is known as a voxel, which represents a value in three-dimensional space on a regular grid [21]. Voxels are extensively used in medical and scientific data, especially in medical CT scans, for visualization and analysis. The number of images for each subject varied from 554 to 734. The thickness, voxels per slice, and pixel resolution on the slice plane were 1.25 mm, 512 by 512, and 0.391–0.841 mm, respectively. Figure 1 shows the acquired axial CT images.

3. Proposed Method

SSM is an established tool in medical image analysis that represents the shape variations of a particular class of shapes. Shape variations are learned from a set of training examples and the variation of shape is represented by using the leading principle components. The process of constructing the SSM of a humeral head involves three steps, including (step 1) segmentation of the humeral head, (step 2) position and orientation alignment, and finally, (step 3) PCA shape analysis. Detailed description of each step is given below.

3.1. Segmentation of the Humeral Head Region

The humeral region is firstly segmented from CT images manually for each subject. Figure 2 shows an example of a manually segmented humerus region from shoulder CT images. To extract the humeral head from the segmented humerus region, the anatomical neck points are acquired manually. Figure 3 shows an example of the acquired point. Next, the anatomical neck points are approximated by a plane using the least squares method. Figure 4 depicts the flow of humeral head extraction from a segmented humerus.

3.2. Position and Orientation Alignment

To align the interindividual position and orientation difference of humeral head shape of each subject, a 3-D anatomical coordinate system is defined. The origin of the coordinate system O = [O_x, O_y, O_z] is the center of the gravity of the anatomical neck points. The most lateral point of the anatomical surface, P_L, is acquired from the plane approximation. An example of the most lateral point can be seen in Figure 3. The X-axis is a vector, which is directed from the origin O to the most lateral point of the anatomical surface P_L. The Z-axis, which is a vector directed from the origin toward the curvature of the humeral head, is the normal vector of the plane. Finally, the Y-axis is the cross product of Z-axis and X-axis. By determining this 3-D anatomical coordinate system, the position and orientation of each subject’s shape is aligned.

The position and orientation of shapes of all subjects are aligned in the same reference space (X-Y-Z left-handed rectangular coordinate system) by affine transformation. Figure 5 shows an example of the coordinate system.

3.3. Principal Component Analysis of Shape Matrix

A distance transform [22], which assigns positive distances from the surface to the outer of the humerus, is applied to the aligned volume of the humeral region. The 6 neighborhood Manhattan distance is used as the distance function. The Manhattan distance d(A, B) between two points $A=\left(x_0,y_0,z_0\right)$ and $B=\left(x_1,y_1,z_1\right)$ is calculated using the equation below:

$$d\left(A,B\right)=\left|x_0-x_1\right|+\left|y_0-y_1\right|+\left|z_0-z_1\right|$$

(1)

The distance on the humerus bone surface voxel is zero. The obtained 3-D distance image of subject i is converted to a 1-D vector to create a feature vector matrix ${\mathit{F}}_{\mathit{i}}$ ($1\times N_v$). The conversion is applied to all training data, and N_s vectors are obtained. N_v is the number of voxels, and N_s is the number of subjects. Equation (1) shows a feature vector matrix F, in which the number of rows equals to the number of voxels in the image N_v (256 × 256 × 256), and the number of columns equals to the number of subjects N_s.

$$\mathit{F}={\left({\mathit{F}}_{\mathbf{1}},\cdots ,{\mathit{F}}_{\mathit{i}},\cdots ,{\mathit{F}}_{{\mathit{N}}_{\mathit{s}}}\right)}^T$$

$$\mathit{F}=\left(\begin{array}{ccccc}f{\left(1\right)}_{1,1,1}& \cdots & f{\left(i\right)}_{1,1,1}& \cdots & f{\left(N_s\right)}_{1,1,1}\\ f{\left(1\right)}_{2,1,1}& \cdots & f{\left(i\right)}_{2,1,1}& \cdots & f{\left(N_s\right)}_{2,1,1}\\ ⋮& \cdots & ⋮& \cdots & ⋮\\ f{\left(1\right)}_{256,256,256}& \cdots & f{\left(i\right)}_{256,256,256}& \cdots & f{\left(N_s\right)}_{256,256,256}\end{array}\right)$$

(2)

In the next step, the average shape matrix ($\overline{\mathit{F}}$) of all subjects is calculated using the following equation

$$\overline{F}=\frac{1}{N_s}{\displaystyle \sum }_{i=1}^{N_s}{F}_i$$

(3)

where F_i is the feature vector matrix of ith subject. Covariance matrix C of feature vector matrix F is calculated by the equation below

$$C=\frac{1}{N_s-1}{\displaystyle \sum }_{i=1}^{N_s}\left(F_i-\overline{F}\right){\left(F_i-\overline{F}\right)}^T$$

(4)

The eigenvalue decomposition of the covariance matrix (C) is used to calculate the eigenvalue $\left(\mathit{\lambda }\right)$ and the eigenvector (V). Covariance matrix C can be expressed in terms of the eigenvalue $\left(\mathit{\lambda }\right)$ and the eigenvector (V) using the equation below:

$$CV=\lambda V$$

(5)

where the eigenvalue $\mathit{\lambda }$ and the eigenvector V can be expressed by the equations below:

$$\lambda =\left[\begin{array}{cccc}{\lambda }_1& 0& \dots & 0\\ 0& {\lambda }_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\lambda }_{N_E}\end{array}\right]$$

(6)

$$V=\left[\begin{array}{cccc}{e}_{11}& e_{12}& \dots & e_{1N_E}\\ e_{21}& e_{22}& \dots & e_{2N_E}\\ ⋮& ⋮& \ddots & ⋮\\ e_{N_{v}1}& e_{N_{v}2}& \dots & e_{N_v{N}_E}\end{array}\right]$$

(7)

where ${\mathit{\lambda }}_1$, ${\mathit{\lambda }}_2$, …, ${\mathit{\lambda }}_{N_E}$ refer to the eigenvalues of corresponding ranks, and N_E represents the maximum number of dimensions of the reduced matrix that is determined by the cumulative contribution rate. e_j = [e_1j,e_2j,…, e_Nvj]^T refers to the eigenvector of jth rank, where each column is arranged by the rate of contribution. The principle component (PC) score (θ_i) of the feature vector of the ith subject is calculated using the following equation

$${\theta }_i=\left(F_i-\overline{F}\right)V$$

(8)

Subject shape ${\tilde{\mathit{F}}}_{\mathit{i}}$ can be reconstructed with the mean shape, eigenvector (V), and PC score ${\mathit{\theta }}_{\mathit{i}}$ as shown in the equation below:

$${\tilde{F}}_i=\overline{F}+{\theta }_i{V}^T$$

(9)

By assigning a PC score (θ_i), an arbitrary shape can be synthesized within the range of individual difference.

In this study, they are synthesized by changing the coefficients along each eigenvector to the mean shape (middle) using the following equation:

$$SSM=\overline{\mathit{F}}+a{\sigma }_1{\mathit{e}}_1+b{\sigma }_2{\mathit{e}}_2$$

(10)

where $\overline{\mathit{F}}$ is the mean shape; ${\sigma }_1$ and ${\sigma }_2$ are the eigenvalues of the first and the second PC, respectively; and a and b are parameters. e₁ and e₂ are the eigenvectors. The reconstructed signed distance map is then converted to a binary image.

3.4. Shape Parameter Definition Using Statistical Shape Model

To quantitatively evaluate the interindividual variations of humeral head shape, the following measurements are performed. Figure 5 shows the anatomical direction of the humeral head. At first, the thickness of the humeral head is measured, which is in the superior-inferior (SI) direction as shown in Figure 5a. In the projection of the Z-axis direction, which is the SI direction, shown in Figure 5b, the vertices that are parallel to the Y-axis and have the longest line segment connecting the two boundary points are defined as P₁ and P₂, respectively. The vertices that are parallel to the X-axis and have the longest line segment connecting the two boundary points are defined as P₃ and P₄, respectively. The point of intersection of line segments P₁–P₂ and P₃–P₄ is point O. Next, the distances between two points (i.e., P₁–P₂, P₃–P₄, P₁–O, P₂–O, P₃–O, and P₄–O) are calculated.

4. Experimental Results

4.1. Construction of Humeral Head SSM and Its Performance Evaluation

This study involved a total of 22 subjects, as described earlier in the materials section. The signed distance was used to create the feature vector matrix. Eigenvector, mean vector, and PC scores were calculated using PCA. Figure 6 shows shape variations of the humeral heads synthesized by using the SSM. The horizontal direction of Figure 6 corresponds to the first PC, whereas the vertical direction corresponds to the second PC. Figure 7 shows the cumulative contribution rate (CCR) corresponding to the number of PCs. It shows that the SSM can represent around 95% of shape variation by using the first five PCs.

An SSM can express the interindividual variation of the shape by the mean (i.e., average shape), the eigenvectors, and the PC scores (i.e., shape parameter). It can also reconstruct the original individual shape [5]. The shape of one subject was chosen as the test data, whereas the shapes of all other subjects were used to build an SSM. The shape parameter of the test data was calculated by PCA. The test data were reconstructed by using the SSM and the calculated parameter by Equation (8). Therefore, verification of the generalization ability [23] of the SSM was evaluated by comparing the test data shapes with the reconstructed one. The Sørensen–Dice coefficient (DSC) was used as an evaluation metric. Figure 8 shows a schematic diagram of the overall process.

The leave-one-out cross validation (LOOCV) procedure [23] was used to validate the generalization ability of the SSM. As the total subjects were 22, LOOCV became a 22-fold cross validation process in this case. For every left-out subject (test subject), a new SSM was computed from the remaining training subjects. The parameters of the left-out subject were computed to reconstruct the subject. In this experiment, a total of 22 different evaluations were available. The average DSC was 0.918. The selected number of PCs was six because the cumulative contribution ratio of the PCs was 95% or more.

4.2. Humeral Head Shape Analysis Using Statistical Shape Model

To evaluate the interindividual humeral head shape quantitatively, the measurements were done for 25 statistical shape models constructed by changing the standard deviation along the first PC from $-3{\mathsf{\sigma }}_1$ to $+3{\mathsf{\sigma }}_1$ at intervals of $0.25{\mathsf{\sigma }}_1$, where ${\mathsf{\sigma }}_1$ is the Eigen value of the first PC. The measurement was performed on the shape of each subject.

Table 1. The average and the variance of the measured values between different points on the humeral head (in mm).

Measured Points	Proposed Method (SSM)	Original Data
The thickness	17.2 ± 5.7	17.6 ± 2.3
P1–P2	43.2 ± 20.5	44.0 ± 8.1
P3–P4	41.4 ± 24.2	41.9 ± 7.4
P1–O	20.1 ± 20.3	19.9 ± 15.8
P2–O	23.1 ± 3.9	24.2 ± 10.3
P3–O	18.7 ± 12.8	18.4 ± 10.3
P4–O	22.6 ± 5.5	23.4 ± 10.0

shows the average and variance of the measured values between different points. The results showed that the first PC had a large variance in the distance between P₁–P₂ and P₃–P₄, i.e., the interindividual variation in the anterior-posterior to lateral-medial (AP-LM) direction was large. In particular, it was found that the variance of the distance P₁–O was large. The results showed that the variance in the distance P₁–O was large in both of the measurements, i.e., measurement from the SSM and measurement from the subject’s shape. Thus, it confirmed that the SSM by the proposed method can grasp the tendency of the subject population. Table 1 also shows the measured values for original data. Figure 9 and

Table 2. The correlation coefficient between thickness and each measured values of different points on the humeral head.

	P1–O	P2–O	P3–O	P4–O
The proposed method (SSM)	0.87	0.24	0.89	0.67
Original data	0.40	0.06	0.26	0.20

show the correlation coefficients between the thicknesses and the distances between different points, i.e., P₁–O, P₂–O, P₃–O, and P₄–O, considering both the measurements. These results showed that the thicknesses and distances of P₁–O and P₃–O had a strong positive correlation (>0.8) when using the proposed method, by the SSM. In other words, it could be considered that the prosthesis should be deformed by considering the same tendency at each size, because the head is not just a scale-up, but rather the extension of anterior and lateral is dominant.

5. Discussions

In glenohumeral osteoarthritis treatment, one of the most preferred ways is to replace the humeral head surgically with an artificial one. Popular consensus among experts is that accurate replication would bring better functional outcome [4,24,25,26]. Although it is well documented that the humeral head is ovoid [4,24,25,26,27], the use of spherical artificial humeral heads during the reconstructive surgery of the shoulder is general practice. Even small alterations (4–5 mm) of size and position of the articular surface could cause adverse effects on the biomechanics of the glenohumeral during arthroplasty surgery on the glenohumeral [4,28,29]. A potential concern arises from this, as the spherical artificial humeral head often exceeds a 4 mm mismatch, when compared with the anatomical measurement of a normal humeral head [30]. Thus, this study proposed a method to construct a humeral head SSM for the design of a patient-specific artificial humeral head. The generalization ability of SSMs was verified using the LOOCV technique, and the results indicated that its generalization ability was quite strong. The study also investigated the interindividual shape of the humeral head quantitatively using the SSM. The investigation demonstrated a strong correlation between individual characteristics of the humeral head and the thickness of the humeral head, which the conventional method failed to identify by simply analyzing the subject data. In line with the hypothesis, this newly obtained knowledge might be useful in designing patient-specific artificial humeral heads. On the contrary, the reliability of the study was impacted by the absence of a sufficient number of subjects. On top of that, the methodological choices were constrained by manual segmentation, which was a very tedious and inefficient process, especially for the persons involved in it. Thus, to improve the reliability and acceptability of the proposed method, the number of subjects should be increased. In addition, the manual segmentation process should be replaced by automatic segmentation.

6. Conclusions

One of the most important criteria for a successful glenohumeral osteoarthritis treatment is that the shape of the artificial humeral head be similar to the actual shape of the humeral head from the viewpoint of the ROM. This study proposed a method to construct a humeral head SSM from CT images. The SSM expressed the interindividual shape variation of the humeral head, and the generalization ability of the SSM was validated by calculating the Dice coefficient, which was 0.918. Humeral head shape analysis using the SSM was conducted, and the anatomical features of the humeral heads were obtained in that process. The quantitative analysis of the humeral head SSM discovered the fact that the thickness of the humeral head had a strong positive correlation with the anterior and lateral expansion of the humeral head. The encouraging segmentation results and the information obtained by quantitative humeral head SSM shape analysis could be considered as an important tool towards successful design of artificial humeral heads. Future works should be directed to increase the reliability of the model by increasing the number of subjects, optimize the extraction method of the humeral head, and elucidate the differences of interindividual variations of humeral head shape in other principal components.