Home
Search
Collections
Journals
About
Contact us
My IOPscience
Automatic bone segmentation and bone-cartilage interface extraction for the shoulder joint from magnetic resonance images
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Med. Biol. 60 1441 (http://iopscience.iop.org/0031-9155/60/4/1441) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 137.149.200.5 This content was downloaded on 21/02/2015 at 02:31
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) 1441–1459
Physics in Medicine & Biology doi:10.1088/0031-9155/60/4/1441
Automatic bone segmentation and bone-cartilage interface extraction for the shoulder joint from magnetic resonance images Zhengyi Yang1, Jurgen Fripp2, Shekhar S Chandra1, Aleš Neubert1, Ying Xia1,2, Mark Strudwick1, Anthony Paproki1,2, Craig Engstrom3 and Stuart Crozier1 1
School of Information Technology and Electrical Engineering, The University of Queensland, Brisbane, Australia 2 The Australian e-Health Research Centre, CSIRO Digital Productivity Flagship, Brisbane, Australia 3 School of Human Movement Studies, The University of Queensland, Brisbane, Australia E-mail: [email protected] Received 13 May 2014, revised 31 October 2014 Accepted for publication 18 November 2014 Published 22 January 2015 Abstract
We present a statistical shape model approach for automated segmentation of the proximal humerus and scapula with subsequent bone-cartilage interface (BCI) extraction from 3D magnetic resonance (MR) images of the shoulder region. Manual and automated bone segmentations from shoulder MR examinations from 25 healthy subjects acquired using steadystate free precession sequences were compared with the Dice similarity coefficient (DSC). The mean DSC scores between the manual and automated segmentations of the humerus and scapula bone volumes surrounding the BCI region were 0.926 ± 0.050 and 0.837 ± 0.059, respectively. The mean DSC values obtained for BCI extraction were 0.806 ± 0.133 for the humerus and 0.795 ± 0.117 for the scapula. The current model-based approach successfully provided automated bone segmentation and BCI extraction from MR images of the shoulder. In future work, this framework appears to provide a promising avenue for automated segmentation and quantitative analysis of cartilage in the glenohumeral joint. Keywords: automatic segmentation, shoulder joint, active shape model (Some figures may appear in colour only in the online journal) 0031-9155/15/041441+19$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
1441
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
1. Introduction Magnetic resonance (MR) imaging provides clinically useful visualization of the shoulder region enabling detailed evaluation of soft-tissue and bone structures often prone to injury or pathological changes (La Rocca Vieira et al 2012, McMonagle and Vinson 2012, Vinson et al 2012). Typically, MR imaging of the muscles, cartilages and bones of the shoulder involves multiplanar examinations and, unlike computed tomography (CT), MR images are acquired without the use of ionising radiation. The segmentation of bone from 3D MR images of the shoulder is a developing area and has a number of potential applications including analyses of osteochondral defects following glenohumeral dislocation, pre-operative evaluation of glenoid geometry and version to inform total joint arthroplasty and as a precursor to segmentation of articular cartilage to quantify tissue damage due to trauma, inflammation or degenerative disease. In this paper, we present a fully automatic approach for bone segmentation in shoulder MR images focusing on the proximal humerus (head, anatomical neck, greater and lesser tuberosities) and adjacent portions of the scapula (glenoid fossa, acromion and coracoid processes). Specifically, our bone segmentation technique is based around an active shape model (ASM) approach successfully applied to 3D MR images of the knee and hip joints (Fripp et al 2007, Schmid et al 2011, Xia et al 2013, Chandra et al 2014); the present segmentation pipeline uses multiple statistical shape models (SSMs) consisting of a combined SSM to help define the spatial relationships between the humerus and scapula, and individual bone SSMs for the humerus and scapula. The segmentation process was initialized by an automated joint centre locator algorithm (Xia et al 2013, Chandra et al 2014) and incorporated a bone-cartilage interface (BCI) extraction algorithm originally proposed for the knee (Fripp et al 2010) for future development of a fully automatic cartilage segmentation pipeline for the glenohumeral joint. 1.1. Bone segmentation of the shoulder region from radiological studies
Previous applications of bone segmentation for the shoulder region from radiological studies have included morphometric evaluations of bone shape prediction in anthropometric analyses (Yang et al 2008), investigation of joint kinematics (Zhu et al 2012), joint interval measurement related to impingement diagnosis and prognosis (Hekimoğlu et al 2013) and imageguided shoulder surgery (Chaoui et al 2011). A spectrum of segmentation methods have been employed across various studies using modalities such as fluoroscopic imaging (Zhu et al 2012), CT (Chaoui et al 2011) and MR imaging. In CT images, the high contrast between bone and other tissues is well suited for automated bone segmentation based on thresholding (Kawasaki et al 2008). However, pathological changes and non-uniform density of the bone may complicate the segmentation process, often causing Hounsfield value thresholding to fail or lead to poor results (Krekel et al 2010). In CT images of severely damaged arthritic shoulder joints, Krekel et al used a semi-automatic method for combining surface curvature and iterative refinement via a Hough transform for bone segmentation (Krekel et al 2010). Comparison of the bone volumes for the proximal humerus and scapula derived from the manual and semi-automated measures found a median Dice Similarity Coefficient (DSC) (Dice 1945) of 98.9% and median Hausdorff distances of 3.06 mm and 3.92 mm, respectively. However, the disadvantages of such CT analyses include the use of ionizing radiation to acquire images of the shoulder (upper thoracic) region and, in comparison with MR imaging, the poorer visualization of important soft-tissues associated with the glenohumeral joint. In previous MR imaging studies of the shoulder, bone segmentation has been primarily performed using manual approaches (Massimini et al 2011, Izadpanah et al 2012). Manual 1442
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
segmentations are both time- and expertise-intensive and inherently analyst dependant given the subjective interpretations involved in visually distinguishing between bones and surrounding tissues in clinical MR images. In order to overcome a number of the above issues related to manual segmentation approaches (e.g. analyst time, training, experience, intra- and inter-rater reliability), a number of semi- and fully-automated approaches have been developed for bone segmentation of the shoulder region based on 2D and 3D MR examinations. As per manual bone segmentations from MR images, such automated approaches have to perform segmentations in the presence of variable tissue boundary strengths (contrast) over the examination/s and various signal inhomogeneities related to imaging sequence/system factors. Essentially, these more automated approaches can be classified in relation to the (non)use of shape priors for driving the bone segmentations. 1.1.1. Bone segmentation without shape priors. In the absence of knowledge about an object’s shape, segmentation is generally based on either pixel intensities (e.g. Otsu thresholding and k-means classification), edges defined by intensity gradient (e.g. edge tracking), or object connectivity inferred from regional intensity similarity (e.g. watershed and region-growing) (Pham et al 2000). Several of these algorithms have been applied to bone segmentation in the shoulder region. Graichen et al (2000) constructed 3D models of the humerus, clavicle and scapula for motion analysis using a region-growing based semi-automatic segmentation method followed by shape-based interpolation from 2D MR images acquired on a 0.2 T scanner; qualitative or quantitative validation on the segmentation accuracy was not reported. Branzan-Albu et al segmented the humerus and scapula from 2D MR images using the isolabel method originally developed for tumour segmentation (Branzan Albu et al 2002) to obtain smooth 3D surface reconstructions of bone elements of the shoulder; segmentation accuracy was not reported. Nhat Tan et al (2007) extended this work by proposing a semi-automatic method using combined region-based and gradient-based supervised segmentation which was applied to contour delineation of the humerus from 2D T1-weighted images acquired at 1.5 T. An area overlap ratio of 96.32% between the semi-automated measures and manual segmentations by a radiologist was achieved on 60 slices. Pérez et al (2009) implemented a 4-class (bone, muscle, cartilage and fat) classification method combining statistical and geometrical model approaches based on an ‘active contours without edge’ framework to segment 2D coronal shoulder images. The results on two 2D images were qualitatively evaluated and the proposed method was found to outperform k-means and Bayesian density mixture approaches. In a motion analysis study of shoulder joint rotation via image segmentation, Koishi et al obtained contours of the humerus and the scapula initially with automatic segmentation followed by manual modification using 3D Fast Low Angle SHot (FLASH) MR images acquired at 1.5 T, details of the segmentation method were not reported (Koishi et al 2011). 1.1.2. Bone segmentation with shape priors. Bone segmentation in clinical MR images
using approaches based solely on image intensities is challenging. The use of shape priors significantly improves segmentation accuracy, particularly when image intensity information is unreliable (Schmid and Magnenat-Thalmann 2008). Shape priors can be represented using SSMs and Statistical Appearance Models (SAM) as reviewed in (Heimann and Meinzer 2009, Hufnagel 2011). In recent years, a number of shape model based methods have been developed for segmenting the bones (and cartilages) in joints such as the knee (Williams et al 2003, Fripp et al 2007) or hip (Xia et al 2013). Articulated SSMs have been proposed for determination of the relative motion between bony structures in the knee (Bindernagel et al 2011) and the hip (Kainmueller et al 2009). In bone segmentation of the hip joint, Schmid et al developed a robust SSM based method to deal with the difficulties in building shape 1443
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
models from small field of view images (Schmid et al 2011). Kainmueller et al have also implemented a deformation model allowing each vertex to move in any direction within a surrounding sphere to handle highly curved surface meshes as relevant to ball-and-socket joints such as the hip and shoulder (Kainmueller et al 2013). More recently, Chandra et al have developed focused shape models, where shape priors were spatially up-weighted to regions of interest, to improve the accuracy of bone segmentation and BCI extraction in the hip joint (Chandra et al 2014). There has been limited work on the use of shape model approaches for bone segmentation in the shoulder region. In a motion analysis study, Liu et al (2008) employed the live wire method for interactive bone segmentation to generate a rigid deterministic 2D shape model of bony elements at one position of the joint. Subsequent images were segmented as a search for the same bone by minimizing an energy function that utilizes both boundary and regionbased information. This work was based around a robust SSM approach (Schmid et al 2011) and involved an initialization step tailored for ball-and-socket joints whereby the concept of a kinematics-inspired initialization was proposed as highly suitable for the shoulder. To the best of the authors’ knowledge, no fully automatic bone segmentation approach using 3D shape priors has been reported for shoulder MR images. 1.1.3. Bone segmentation based on active shape models. The application of ASMs (Coo-
tes et al 1995) has been used successfully for bone segmentation from MR images of the knee (Fripp et al 2007) and hip joints (Schmid et al 2011, Xia et al 2013). In ASM-based approaches, shape priors of the object to be segmented are represented by an SSM and the object boundary is found by iteratively deforming the SSM to fit in the image. Once an SSM has been built, the initialization of ASMs plays an important role in its convergence at the object boundary. A coarse-to-fine strategy for bone segmentation in the knee was developed to use atlas selection for ASM initialization (Park et al 2012). In this paper, we present an ASM-based automatic bone segmentation method using hybrid SSMs, consisting of a combined SSM incorporating the spatial relationship between the humerus and scapula, and individual SSMs for each of these bones. The initialization was performed using an automated joint centre locator algorithm, similar to that in previous studies (Xia et al 2013, Chandra et al 2014). 1.2. Extraction of the bone-cartilage interface
The BCI is defined as the region on the external bone surface that is covered with articular cartilage tissue. In this study, we adopted the BCI extraction algorithm originally proposed for the knee (Fripp et al 2010) to quantitatively evaluate its performance with the aim of developing a fully automatic shoulder cartilage segmentation pipeline. This approach is based on our previous work where the extraction of the BCI after bone segmentation in the knee provides a robust reference for automatic segmentation and morphological analysis of articular cartilage (Fripp et al 2010). The overarching motivation behind developing a system for automatic cartilage morphometric analysis stems from contemporary work on new therapeutic measures including cartilage repair and pharmaceutical interventions. A major focus of research for automatic cartilage segmentation has centred on MR imaging given its’ capacity to visualize cartilage morphology, to provide measures on parameters such as volume, thickness, surface area and curvature and investigate the biochemical composition of cartilage (Link et al 2007). Currently, the majority of the research centres on the knee (Folkesson et al 2007, Brem et al 2009, Yin et al 2010, Swamy and Holi 2012, Tamez-Pena et al 2012, Marques et al 2013) and hip joints (Cheng et al 2011). Studies on shoulder 1444
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
cartilage morphometry have used interactive measurements (Hekimoğlu et al 2013) and manual segmentations (Massimini et al 2011). A 2D semi-automatic glenohumeral cartilage segmentation method based on B-Spline Snakes with interactive initialization was developed (Stammberger et al 1999a, 1999b). 2. Materials and methods 2.1. MR image acquisition
Twenty-five healthy subjects (12 females and 13 males, 33.6 ± 9.1 year-old, 62–95 kg) with no history of shoulder injury underwent 3 T MR imaging (Siemens Trio) of the right shoulder. Images were acquired using a dedicated shoulder coil (4 element phased array coil) with a true fast imaging with steady state precession (TrueFISP) sequence. All subjects were imaged in the supine position with the right hand placed palm down under the ipsilateral hip, both for immobilisation and reproducibility. A standardized imaging plane was aligned through the humeral head and parallel to the spine of the scapula. The imaging parameters used were: TR/TE = 9.23/3.99 ms, iPAT = 2, field of view (FOV) 140 × 140 mm2. The image dimensions were 350 × 350 × 160 with a voxel size of 0.4 × 0.4 × 0.47 mm3. The in-plane resolution used in subsequent analysis was 0.2188 mm isotropic achieved by sinc interpolation. The medical research ethics committee of the University of Queensland approved the study and informed consent was obtained from all participants. The image datasets were anonymized. 2.2. Manual segmentation of the MR images
For each MR image, the contours of the humerus and scapula were manually segmented on each slice (160 coronal slices) by ZY, with expert guidance from CE using ITK-SNAP (Yushkevich et al 2006). To acquire BCI priors for developing an automatic BCI extraction method, the cartilages covering the humeral head and glenoid were manually segmented by MS, an experienced radiographer, on all subjects using the coronal slices as the primary basis for analyses. To evaluate intra- and inter-observer variability, cartilage segmentation was repeated by MS and CE on 8 randomly selected subjects. 2.3. Segmentation pipeline
A central issue with segmentation of the bones and cartilages of the shoulder joint from clinical MR images (TrueFISP in this study) relates to weak or missing boundaries between adjacent tissues. As demonstrated in figure 1(a), bone must be distinguished from adjacent fat, tendons and background, while the cartilage needs to be separated from the synovial fluid, labrum, ligaments and muscles associated with the glenohumeral joint. In these situations of multi-tissue interfaces, segmentation with shape priors, such as those represented as SSMs, has proven advantages over methods based solely on image intensity information. Following the pipeline of bone segmentation and BCI extraction in the knee developed previously (Fripp et al 2007), the fully automated method for the shoulder joint consists of four stages: (i) preprocessing, (ii) initialization, (iii) ASM-based bone segmentation and (iv) BCI extraction (figure 2). (i) In preprocessing, the first step is bias field correction for intensity inhomogeneity caused by B0 non-uniformity and B1 non-linearity. The ITK (Ibanez et al 2003) implementation of the N4 algorithm based on B-spline fitting was employed in this study (Tustison et al 2010) 1445
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
Figure 1. Example slices of MR images overlaid with manual bone segmentations
of the proximal humerus (red) and scapula (blue) on the left panel, and humeral head cartilage (cyan) and glenoid cartilage (green) on the right panel. Axial and coronal views are shown in the top and middle rows, respectively. Volume or surface rendered pictures of the segmented objects are shown in the bottom row. (a) Examples of individual axial and coronal images and a volume rendered MR image block encompassing the glenohumeral joint region. The cartilage has a relatively bright signal intensity, while the adjacent bone, tendons and adipose tissue have dark signal intensities; (b) example slices with superimposition of manually segmented humerus and scapula with 3D rendered volume of these bones; (c) cartilage segmentation by CE with resulting volume rendering demonstrating the larger humeral and smaller glenoid cartilage plates; (c) cartilage segmentation by MS (first measurement occasion); (d) cartilage segmentation by MS (retest measurement occasion). The arrows in (c)–(e) illustrate several instances of inter- and intra-observer differences in manual segmentations.
with the following parameters: spline order 3, node distance 150, sigmoid alpha 0, sigmoid beta 0.5, maximum number of iterations at each resolution 100 × 80 × 60, convergence threshold 0.0001 and shrink factor 3. The second step is to remove high frequency noise by smoothing. We tested median filtering, mean filtering and anisotropic diffusion; median filtering was chosen due to the balance between de-noising outcomes and computational efficiency. (ii) The initial position of the SSM is critical to the segmentation accuracy and convergence speed. The initialization involves propagating the known bone surfaces of an atlas image to the patient space via affine registration as implemented in ITK. Intensity-based affine registration was used with normalized mutual information as similarity metric and tri-linear interpolation as interpolator. (iii) The ASM-based segmentation approach involved three steps: coarse segmentation using a combined SSM, fine segmentation using individual SSMs of the proximal humerus and scapula and final segmentation using bone surface relaxation, where the surface is free-form deformed without shape priors. This hybrid SSM-based segmentation scheme is similar to that used for segmentation of the articulating bones in the knee (Fripp et al 2007) and hip joints (Xia et al 2013). Details are provided below in subsection 2.3.2 after a brief review on SSM construction in subsection 2.3.1. (iv) The BCI extraction (see 2.3.3.) followed the above processes. 1446
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
Preprocessing Raw Image Bias Correction
Median Filtering
Preprocessed Image
Initialization T
Initialized Surface
Affine Registration Atlas Image
Transformation Atlas Surface ASM-based Segmentation Surfaces
Labels
Coarse Segmentation using Combined Model
Fine Segmentation using Separate Models
Final Segmentation after Model Relaxation
BCI Extraction
Initial BCI Extraction
Tissue Estimation
Until Converge
BCI growing
Figure 2. Flowchart of the ASM-based bone segmentation and BCI extraction pipeline.
2.3.1. Statistical shape model construction. The shapes of the manually segmented bones were represented using a point distribution model (PDM) (Cootes et al 1995). The correspondence of each point across the training set of triangular meshes was established using a template surface optimization algorithm (Chandra et al 2014). Denote the training shapes as 1447
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
Figure 3. Statistical shape models of the humerus (left) and scapula (right). The first (top) and second (bottom) modes of variation are illustrated by shapes generated at ±2 standard deviations of each mode, and by directional vectors describing vertex displacement associated with each mode. The vectors are superimposed on the mean shape and colour-coded by the magnitude of normalized relative variation.
xi, i = 1, 2, ..., N , where N is the number of shapes. The mean shape x and covariance matrix C were computed as: x =
1 N
N
∑ xi, C = i=1
N
1 ∑ (xi − x) (xi − x)T . N − 1 i=1
(1)
Principal component analysis (PCA) was then applied to C yielding the modes of variation as the eigenvectors vi and their associated variances as the eigenvalues wi. A shape instance y can be approximated by: y = x + Vb,
(2)
where, V is the matrix with the first t modes of shape variations as columns and b is a weighting vector controlling the influence of each mode on the reconstructed shape (Cootes et al 1995). As illustrated in figure 3, individual bone SSMs for the humerus and scapula, along with a ‘combined’ SSM, constructed to address the relative spatial relationship between the humerus and scapula, were built in this study. 2.3.2. ASM-based segmentation. The automated bone segmentation was performed in a
three-level multi-resolution Gaussian pyramid (Fripp et al 2007). The ASM segmentation was performed by iterating the following two steps until the stopping criteria were satisfied (i.e. surface change was less than a threshold or the maximum number of iterations was reached): • Deformation. The surface was deformed by moving each vertex independently along its normal towards the mean gradient of the 1D intensity profile determined at that vertex along its normal. • Shape constraints. The pose and shape parameters, as encoded in b, of the deformed surface were estimated using SSM to regularize anatomically unlikely shapes. 1448
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
After initialization of the model parameters, a coarse-to-fine approach was taken to segment the proximal humerus and the scapula: • Coarse segmentation using the combined SSM which contained the spatial relationship between the two bones. This step was used for improving the atlas-based affine initialization. • Fine segmentation using the individual SSM of the humerus and the scapula. • Final segmentation involving bone surface relaxation (Fripp et al 2007) to fine-tune the bone boundaries by allowing the surface to free-form deform for a few iterations, not regularized by shape constraints but volume preserved smoothing. 2.3.3. BCI extraction. From the manually segmented cartilages, prior knowledge about the
BCI can be acquired by computing the probability p of each point on the bone surface being covered by cartilage. This probability was estimated as the frequency of cartilage tissue appearing above a given point on the bone surface in all of the training images. After bone segmentation, the BCI is automatically detected and extracted in a fashion of seeded region growing. The highly possible points on BCI (e.g. pi > 0.95, where i is the index to the points) were selected as seed BCI. Tissue intensity above the seed BCI were sampled along the normal to the surface and modelled as a Gaussian distribution. The distribution was used to grow the BCI in the following way. Every point with at least one neighbour on the current BCI was examined for the presence of cartilage tissue in the extracted intensity profile along its normal. If more than half of the profile is cartilage tissue, according to the estimated Gaussian distribution, this point was included as BCI. Estimation of tissue intensity distribution and BCI growing were repeated until convergence. 2.4. Validation method
To validate the developed automatic bone segmentation method for the shoulder, a leaveone-out strategy was employed. To segment the individual bones (humerus, scapula) on a given subject’s image, the manual segmentations on the remaining 24 datasets were used to build an SSM for segmentation. There were 25 runs of SSM construction and ASMbased segmentation. In each run, the segmentation accuracy was quantified using DSC values, mean absolute surface distance (MASD) and relative volume difference (RVD) measures (Xia et al 2013) with the manual segmentation as the reference. These quantities were computed before and after bone segmentation relaxation. Calculation of the DSC is defined as: DSC =
2TP , (2TP + FP + FN )
(3)
where, TP, FP, and FN are true positives, false positives and false negatives, respectively. They can be calculated using set operations on the automatic segmentation (A) and the reference (M). The DSC value indicates the relative true positive ratio with a larger value indicating better performance (1 is the maximum value). The MASD d is the average distance between points on two surfaces and defined as: d (A, M ) =
davg (SA, SM ) + davg (SM , SA) 2
1449
,
(4)
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
where, davg (SA, SM ) is the average directed surface distance from all points on surface SA to surface SM . The RVD is the volume difference in percentage of automatic segmentation A with respect to the reference volume M: RVD (A, M ) =
A − M × 100% . M
(5)
To evaluate any improvement in terms of volume overlap achieved by the relaxation step, paired t-tests were employed to examine for differences between the DSC values before and after relaxation. The segmentation accuracy of regions covered by the humeral and glenoid cartilages is of particular interest, given that the performance of the follow-up step of BCI extraction will be substantially affected by this step. To investigate the segmentation error distribution over the bone surfaces, with a particular focus on BCI regions, Hausdorff distance maps (Xia et al 2013) were generated for both the humerus and scapula. The accuracy of BCI extraction was evaluated using a surface equivalent to DSC values (Fripp et al 2007), measuring the agreement between the automatically and manually extracted BCI. To reduce the sampling error, the reference BCI was determined from the resampled manual cartilage segmentation using shape-based interpolation with a slice thickness of 0.125 mm. The TP was defined as the surface area where the point-to-point or point-to-surface distance between the extracted BCI to the reference BCI was less than 0.5 mm. This distance was calculated in both directions and the average was used. The FP region was defined as the surface area on the extracted BCI that was more than 0.5 mm away from any point on the reference BCI, while the FN region was defined as the surface area on the reference BCI that was more than 0.5 mm away from any points on the extracted BCI. Eight subjects were randomly selected and manual cartilage segmentations were performed by two expert analysts (MS and CE). After appropriate consensus to set rating parameters, the MR examinations were segmented by these experts in an independent fashion. Taking the first segmentation from expert MS as M, we simply replace A with BCI determined from a repeated manual segmentation performed by the same expert, or the other expert CE, for estimating intra- and inter-rater reliability, respectively. 3. Results 3.1. Statistical shape modelling
The SSMs built from the manual bone segmentations are shown in figure 3. Using the shrinkwrap surface extraction algorithm (Fripp et al 2007), segmented surfaces of the proximal humerus and the scapula were extracted from the 25 cases with each of the individual bone surfaces having 8,192 points and 16 380 triangles. The humeral, scapular and combined SSMs had different compactness (i.e. the ability to represent shape variation using as few modes as possible) and are plotted in figure 4 as the reconstruction precision against the number of modes. Here, reconstruction precision was defined as the cumulative relative shape variance in percentage terms. For example, the variation of shape accounted for by the first 5 modes was 72% for the humerus, 51% for the scapula and 78% for the combined humerus-scapula SSM. In this work, we used SSMs which accounted for greater than 95% of shape variations which equated to the first 16, 20 and 15 modes of the humeral, scapula and combined models, respectively. 1450
Z Yang et al
1
1
0.9
0.9
0.8
0.8
0.8
0.7
0.7
0.7
0.6 0.5 0.4
Precision
1 0.9
Precision
Precision
Phys. Med. Biol. 60 (2015) 1441
0.6 0.5 0.4
0.6 0.5 0.4
0.3
0.3
0.3
0.2
0.2
0.2
0.1
0.1
0 1
5
9
13
17
21
25
0 1
0.1 5
9
13
17
21
25
0 1
5
9
13
Mode
Mode
Mode
(a)
(b)
(c)
17
21
25
Figure 4. The compactness of each SSM constructed. (a) Humerus; (b) scapula; (c) combined.
3.2. Bone segmentation
An example (case #3) of the intermediate and final results from the automated segmentation of the humerus and scapula is given in figure 5. After the initialization and coarse segmentation stage, the placement of the combined ASM using the combined model is anatomically realistic, as shown in figure 5(b), although in the provided example, the joint space in the glenohumeral cartilage region appears overly narrow. This is most likely due to the spatial relationship encoded in the combined model not adequately expressing the variation present in this particular case. Following this step in the ASM process, the subsequent fitting of the individual humerus and scapula SSMs universally resulted in improved segmentation of these bones in relation to their corresponding anatomical profiles in the MR images for the case demonstrated in figure 5(c). Specifically, in this stage of the pipeline, the fitting of the individual humerus model (table 1) resulted in DSC values in the range of 0.684–0.947 with a mean of 0.909. In the more anatomically complex scapula, the range of DSC values was 0.549–0.825 with a mean of 0.723. The final relaxation step provided a localized fine-tuning of the bone boundaries to modify instances of likely under- or over-segmentation of the humerus and scapula, as demonstrated in figure 5(d). At the end of this pipeline, there were 23/25 and 5/25 cases having DSC values above 0.85 for the humerus (0.926 ± 0.050) and scapula (0.806 ± 0.054), respectively. In line with our future aim of developing a robust automated pipeline for cartilage segmentation in the glenohumeral joint, we calculated the DSC values for the scapula in an identified region of interest (ROI) centred around the glenoid fossa which encloses the BCI region of this bone, which reduced the effects of the large anatomical variance of the scapula outside this area, resulting in 15/25 cases having DSC values above 0.85 with an overall mean of 0.837. A boxplot showing the DSC values for different segmentation stages is given in figure 6. The improvements in terms of DSC values between the manual and automated segmentations of bone volume (table 1) achieved by relaxation were statistically significant for both the humerus (t = 12.104, p
Search
Collections
Journals
About
Contact us
My IOPscience
Automatic bone segmentation and bone-cartilage interface extraction for the shoulder joint from magnetic resonance images
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Med. Biol. 60 1441 (http://iopscience.iop.org/0031-9155/60/4/1441) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 137.149.200.5 This content was downloaded on 21/02/2015 at 02:31
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) 1441–1459
Physics in Medicine & Biology doi:10.1088/0031-9155/60/4/1441
Automatic bone segmentation and bone-cartilage interface extraction for the shoulder joint from magnetic resonance images Zhengyi Yang1, Jurgen Fripp2, Shekhar S Chandra1, Aleš Neubert1, Ying Xia1,2, Mark Strudwick1, Anthony Paproki1,2, Craig Engstrom3 and Stuart Crozier1 1
School of Information Technology and Electrical Engineering, The University of Queensland, Brisbane, Australia 2 The Australian e-Health Research Centre, CSIRO Digital Productivity Flagship, Brisbane, Australia 3 School of Human Movement Studies, The University of Queensland, Brisbane, Australia E-mail: [email protected] Received 13 May 2014, revised 31 October 2014 Accepted for publication 18 November 2014 Published 22 January 2015 Abstract
We present a statistical shape model approach for automated segmentation of the proximal humerus and scapula with subsequent bone-cartilage interface (BCI) extraction from 3D magnetic resonance (MR) images of the shoulder region. Manual and automated bone segmentations from shoulder MR examinations from 25 healthy subjects acquired using steadystate free precession sequences were compared with the Dice similarity coefficient (DSC). The mean DSC scores between the manual and automated segmentations of the humerus and scapula bone volumes surrounding the BCI region were 0.926 ± 0.050 and 0.837 ± 0.059, respectively. The mean DSC values obtained for BCI extraction were 0.806 ± 0.133 for the humerus and 0.795 ± 0.117 for the scapula. The current model-based approach successfully provided automated bone segmentation and BCI extraction from MR images of the shoulder. In future work, this framework appears to provide a promising avenue for automated segmentation and quantitative analysis of cartilage in the glenohumeral joint. Keywords: automatic segmentation, shoulder joint, active shape model (Some figures may appear in colour only in the online journal) 0031-9155/15/041441+19$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
1441
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
1. Introduction Magnetic resonance (MR) imaging provides clinically useful visualization of the shoulder region enabling detailed evaluation of soft-tissue and bone structures often prone to injury or pathological changes (La Rocca Vieira et al 2012, McMonagle and Vinson 2012, Vinson et al 2012). Typically, MR imaging of the muscles, cartilages and bones of the shoulder involves multiplanar examinations and, unlike computed tomography (CT), MR images are acquired without the use of ionising radiation. The segmentation of bone from 3D MR images of the shoulder is a developing area and has a number of potential applications including analyses of osteochondral defects following glenohumeral dislocation, pre-operative evaluation of glenoid geometry and version to inform total joint arthroplasty and as a precursor to segmentation of articular cartilage to quantify tissue damage due to trauma, inflammation or degenerative disease. In this paper, we present a fully automatic approach for bone segmentation in shoulder MR images focusing on the proximal humerus (head, anatomical neck, greater and lesser tuberosities) and adjacent portions of the scapula (glenoid fossa, acromion and coracoid processes). Specifically, our bone segmentation technique is based around an active shape model (ASM) approach successfully applied to 3D MR images of the knee and hip joints (Fripp et al 2007, Schmid et al 2011, Xia et al 2013, Chandra et al 2014); the present segmentation pipeline uses multiple statistical shape models (SSMs) consisting of a combined SSM to help define the spatial relationships between the humerus and scapula, and individual bone SSMs for the humerus and scapula. The segmentation process was initialized by an automated joint centre locator algorithm (Xia et al 2013, Chandra et al 2014) and incorporated a bone-cartilage interface (BCI) extraction algorithm originally proposed for the knee (Fripp et al 2010) for future development of a fully automatic cartilage segmentation pipeline for the glenohumeral joint. 1.1. Bone segmentation of the shoulder region from radiological studies
Previous applications of bone segmentation for the shoulder region from radiological studies have included morphometric evaluations of bone shape prediction in anthropometric analyses (Yang et al 2008), investigation of joint kinematics (Zhu et al 2012), joint interval measurement related to impingement diagnosis and prognosis (Hekimoğlu et al 2013) and imageguided shoulder surgery (Chaoui et al 2011). A spectrum of segmentation methods have been employed across various studies using modalities such as fluoroscopic imaging (Zhu et al 2012), CT (Chaoui et al 2011) and MR imaging. In CT images, the high contrast between bone and other tissues is well suited for automated bone segmentation based on thresholding (Kawasaki et al 2008). However, pathological changes and non-uniform density of the bone may complicate the segmentation process, often causing Hounsfield value thresholding to fail or lead to poor results (Krekel et al 2010). In CT images of severely damaged arthritic shoulder joints, Krekel et al used a semi-automatic method for combining surface curvature and iterative refinement via a Hough transform for bone segmentation (Krekel et al 2010). Comparison of the bone volumes for the proximal humerus and scapula derived from the manual and semi-automated measures found a median Dice Similarity Coefficient (DSC) (Dice 1945) of 98.9% and median Hausdorff distances of 3.06 mm and 3.92 mm, respectively. However, the disadvantages of such CT analyses include the use of ionizing radiation to acquire images of the shoulder (upper thoracic) region and, in comparison with MR imaging, the poorer visualization of important soft-tissues associated with the glenohumeral joint. In previous MR imaging studies of the shoulder, bone segmentation has been primarily performed using manual approaches (Massimini et al 2011, Izadpanah et al 2012). Manual 1442
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
segmentations are both time- and expertise-intensive and inherently analyst dependant given the subjective interpretations involved in visually distinguishing between bones and surrounding tissues in clinical MR images. In order to overcome a number of the above issues related to manual segmentation approaches (e.g. analyst time, training, experience, intra- and inter-rater reliability), a number of semi- and fully-automated approaches have been developed for bone segmentation of the shoulder region based on 2D and 3D MR examinations. As per manual bone segmentations from MR images, such automated approaches have to perform segmentations in the presence of variable tissue boundary strengths (contrast) over the examination/s and various signal inhomogeneities related to imaging sequence/system factors. Essentially, these more automated approaches can be classified in relation to the (non)use of shape priors for driving the bone segmentations. 1.1.1. Bone segmentation without shape priors. In the absence of knowledge about an object’s shape, segmentation is generally based on either pixel intensities (e.g. Otsu thresholding and k-means classification), edges defined by intensity gradient (e.g. edge tracking), or object connectivity inferred from regional intensity similarity (e.g. watershed and region-growing) (Pham et al 2000). Several of these algorithms have been applied to bone segmentation in the shoulder region. Graichen et al (2000) constructed 3D models of the humerus, clavicle and scapula for motion analysis using a region-growing based semi-automatic segmentation method followed by shape-based interpolation from 2D MR images acquired on a 0.2 T scanner; qualitative or quantitative validation on the segmentation accuracy was not reported. Branzan-Albu et al segmented the humerus and scapula from 2D MR images using the isolabel method originally developed for tumour segmentation (Branzan Albu et al 2002) to obtain smooth 3D surface reconstructions of bone elements of the shoulder; segmentation accuracy was not reported. Nhat Tan et al (2007) extended this work by proposing a semi-automatic method using combined region-based and gradient-based supervised segmentation which was applied to contour delineation of the humerus from 2D T1-weighted images acquired at 1.5 T. An area overlap ratio of 96.32% between the semi-automated measures and manual segmentations by a radiologist was achieved on 60 slices. Pérez et al (2009) implemented a 4-class (bone, muscle, cartilage and fat) classification method combining statistical and geometrical model approaches based on an ‘active contours without edge’ framework to segment 2D coronal shoulder images. The results on two 2D images were qualitatively evaluated and the proposed method was found to outperform k-means and Bayesian density mixture approaches. In a motion analysis study of shoulder joint rotation via image segmentation, Koishi et al obtained contours of the humerus and the scapula initially with automatic segmentation followed by manual modification using 3D Fast Low Angle SHot (FLASH) MR images acquired at 1.5 T, details of the segmentation method were not reported (Koishi et al 2011). 1.1.2. Bone segmentation with shape priors. Bone segmentation in clinical MR images
using approaches based solely on image intensities is challenging. The use of shape priors significantly improves segmentation accuracy, particularly when image intensity information is unreliable (Schmid and Magnenat-Thalmann 2008). Shape priors can be represented using SSMs and Statistical Appearance Models (SAM) as reviewed in (Heimann and Meinzer 2009, Hufnagel 2011). In recent years, a number of shape model based methods have been developed for segmenting the bones (and cartilages) in joints such as the knee (Williams et al 2003, Fripp et al 2007) or hip (Xia et al 2013). Articulated SSMs have been proposed for determination of the relative motion between bony structures in the knee (Bindernagel et al 2011) and the hip (Kainmueller et al 2009). In bone segmentation of the hip joint, Schmid et al developed a robust SSM based method to deal with the difficulties in building shape 1443
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
models from small field of view images (Schmid et al 2011). Kainmueller et al have also implemented a deformation model allowing each vertex to move in any direction within a surrounding sphere to handle highly curved surface meshes as relevant to ball-and-socket joints such as the hip and shoulder (Kainmueller et al 2013). More recently, Chandra et al have developed focused shape models, where shape priors were spatially up-weighted to regions of interest, to improve the accuracy of bone segmentation and BCI extraction in the hip joint (Chandra et al 2014). There has been limited work on the use of shape model approaches for bone segmentation in the shoulder region. In a motion analysis study, Liu et al (2008) employed the live wire method for interactive bone segmentation to generate a rigid deterministic 2D shape model of bony elements at one position of the joint. Subsequent images were segmented as a search for the same bone by minimizing an energy function that utilizes both boundary and regionbased information. This work was based around a robust SSM approach (Schmid et al 2011) and involved an initialization step tailored for ball-and-socket joints whereby the concept of a kinematics-inspired initialization was proposed as highly suitable for the shoulder. To the best of the authors’ knowledge, no fully automatic bone segmentation approach using 3D shape priors has been reported for shoulder MR images. 1.1.3. Bone segmentation based on active shape models. The application of ASMs (Coo-
tes et al 1995) has been used successfully for bone segmentation from MR images of the knee (Fripp et al 2007) and hip joints (Schmid et al 2011, Xia et al 2013). In ASM-based approaches, shape priors of the object to be segmented are represented by an SSM and the object boundary is found by iteratively deforming the SSM to fit in the image. Once an SSM has been built, the initialization of ASMs plays an important role in its convergence at the object boundary. A coarse-to-fine strategy for bone segmentation in the knee was developed to use atlas selection for ASM initialization (Park et al 2012). In this paper, we present an ASM-based automatic bone segmentation method using hybrid SSMs, consisting of a combined SSM incorporating the spatial relationship between the humerus and scapula, and individual SSMs for each of these bones. The initialization was performed using an automated joint centre locator algorithm, similar to that in previous studies (Xia et al 2013, Chandra et al 2014). 1.2. Extraction of the bone-cartilage interface
The BCI is defined as the region on the external bone surface that is covered with articular cartilage tissue. In this study, we adopted the BCI extraction algorithm originally proposed for the knee (Fripp et al 2010) to quantitatively evaluate its performance with the aim of developing a fully automatic shoulder cartilage segmentation pipeline. This approach is based on our previous work where the extraction of the BCI after bone segmentation in the knee provides a robust reference for automatic segmentation and morphological analysis of articular cartilage (Fripp et al 2010). The overarching motivation behind developing a system for automatic cartilage morphometric analysis stems from contemporary work on new therapeutic measures including cartilage repair and pharmaceutical interventions. A major focus of research for automatic cartilage segmentation has centred on MR imaging given its’ capacity to visualize cartilage morphology, to provide measures on parameters such as volume, thickness, surface area and curvature and investigate the biochemical composition of cartilage (Link et al 2007). Currently, the majority of the research centres on the knee (Folkesson et al 2007, Brem et al 2009, Yin et al 2010, Swamy and Holi 2012, Tamez-Pena et al 2012, Marques et al 2013) and hip joints (Cheng et al 2011). Studies on shoulder 1444
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
cartilage morphometry have used interactive measurements (Hekimoğlu et al 2013) and manual segmentations (Massimini et al 2011). A 2D semi-automatic glenohumeral cartilage segmentation method based on B-Spline Snakes with interactive initialization was developed (Stammberger et al 1999a, 1999b). 2. Materials and methods 2.1. MR image acquisition
Twenty-five healthy subjects (12 females and 13 males, 33.6 ± 9.1 year-old, 62–95 kg) with no history of shoulder injury underwent 3 T MR imaging (Siemens Trio) of the right shoulder. Images were acquired using a dedicated shoulder coil (4 element phased array coil) with a true fast imaging with steady state precession (TrueFISP) sequence. All subjects were imaged in the supine position with the right hand placed palm down under the ipsilateral hip, both for immobilisation and reproducibility. A standardized imaging plane was aligned through the humeral head and parallel to the spine of the scapula. The imaging parameters used were: TR/TE = 9.23/3.99 ms, iPAT = 2, field of view (FOV) 140 × 140 mm2. The image dimensions were 350 × 350 × 160 with a voxel size of 0.4 × 0.4 × 0.47 mm3. The in-plane resolution used in subsequent analysis was 0.2188 mm isotropic achieved by sinc interpolation. The medical research ethics committee of the University of Queensland approved the study and informed consent was obtained from all participants. The image datasets were anonymized. 2.2. Manual segmentation of the MR images
For each MR image, the contours of the humerus and scapula were manually segmented on each slice (160 coronal slices) by ZY, with expert guidance from CE using ITK-SNAP (Yushkevich et al 2006). To acquire BCI priors for developing an automatic BCI extraction method, the cartilages covering the humeral head and glenoid were manually segmented by MS, an experienced radiographer, on all subjects using the coronal slices as the primary basis for analyses. To evaluate intra- and inter-observer variability, cartilage segmentation was repeated by MS and CE on 8 randomly selected subjects. 2.3. Segmentation pipeline
A central issue with segmentation of the bones and cartilages of the shoulder joint from clinical MR images (TrueFISP in this study) relates to weak or missing boundaries between adjacent tissues. As demonstrated in figure 1(a), bone must be distinguished from adjacent fat, tendons and background, while the cartilage needs to be separated from the synovial fluid, labrum, ligaments and muscles associated with the glenohumeral joint. In these situations of multi-tissue interfaces, segmentation with shape priors, such as those represented as SSMs, has proven advantages over methods based solely on image intensity information. Following the pipeline of bone segmentation and BCI extraction in the knee developed previously (Fripp et al 2007), the fully automated method for the shoulder joint consists of four stages: (i) preprocessing, (ii) initialization, (iii) ASM-based bone segmentation and (iv) BCI extraction (figure 2). (i) In preprocessing, the first step is bias field correction for intensity inhomogeneity caused by B0 non-uniformity and B1 non-linearity. The ITK (Ibanez et al 2003) implementation of the N4 algorithm based on B-spline fitting was employed in this study (Tustison et al 2010) 1445
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
Figure 1. Example slices of MR images overlaid with manual bone segmentations
of the proximal humerus (red) and scapula (blue) on the left panel, and humeral head cartilage (cyan) and glenoid cartilage (green) on the right panel. Axial and coronal views are shown in the top and middle rows, respectively. Volume or surface rendered pictures of the segmented objects are shown in the bottom row. (a) Examples of individual axial and coronal images and a volume rendered MR image block encompassing the glenohumeral joint region. The cartilage has a relatively bright signal intensity, while the adjacent bone, tendons and adipose tissue have dark signal intensities; (b) example slices with superimposition of manually segmented humerus and scapula with 3D rendered volume of these bones; (c) cartilage segmentation by CE with resulting volume rendering demonstrating the larger humeral and smaller glenoid cartilage plates; (c) cartilage segmentation by MS (first measurement occasion); (d) cartilage segmentation by MS (retest measurement occasion). The arrows in (c)–(e) illustrate several instances of inter- and intra-observer differences in manual segmentations.
with the following parameters: spline order 3, node distance 150, sigmoid alpha 0, sigmoid beta 0.5, maximum number of iterations at each resolution 100 × 80 × 60, convergence threshold 0.0001 and shrink factor 3. The second step is to remove high frequency noise by smoothing. We tested median filtering, mean filtering and anisotropic diffusion; median filtering was chosen due to the balance between de-noising outcomes and computational efficiency. (ii) The initial position of the SSM is critical to the segmentation accuracy and convergence speed. The initialization involves propagating the known bone surfaces of an atlas image to the patient space via affine registration as implemented in ITK. Intensity-based affine registration was used with normalized mutual information as similarity metric and tri-linear interpolation as interpolator. (iii) The ASM-based segmentation approach involved three steps: coarse segmentation using a combined SSM, fine segmentation using individual SSMs of the proximal humerus and scapula and final segmentation using bone surface relaxation, where the surface is free-form deformed without shape priors. This hybrid SSM-based segmentation scheme is similar to that used for segmentation of the articulating bones in the knee (Fripp et al 2007) and hip joints (Xia et al 2013). Details are provided below in subsection 2.3.2 after a brief review on SSM construction in subsection 2.3.1. (iv) The BCI extraction (see 2.3.3.) followed the above processes. 1446
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
Preprocessing Raw Image Bias Correction
Median Filtering
Preprocessed Image
Initialization T
Initialized Surface
Affine Registration Atlas Image
Transformation Atlas Surface ASM-based Segmentation Surfaces
Labels
Coarse Segmentation using Combined Model
Fine Segmentation using Separate Models
Final Segmentation after Model Relaxation
BCI Extraction
Initial BCI Extraction
Tissue Estimation
Until Converge
BCI growing
Figure 2. Flowchart of the ASM-based bone segmentation and BCI extraction pipeline.
2.3.1. Statistical shape model construction. The shapes of the manually segmented bones were represented using a point distribution model (PDM) (Cootes et al 1995). The correspondence of each point across the training set of triangular meshes was established using a template surface optimization algorithm (Chandra et al 2014). Denote the training shapes as 1447
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
Figure 3. Statistical shape models of the humerus (left) and scapula (right). The first (top) and second (bottom) modes of variation are illustrated by shapes generated at ±2 standard deviations of each mode, and by directional vectors describing vertex displacement associated with each mode. The vectors are superimposed on the mean shape and colour-coded by the magnitude of normalized relative variation.
xi, i = 1, 2, ..., N , where N is the number of shapes. The mean shape x and covariance matrix C were computed as: x =
1 N
N
∑ xi, C = i=1
N
1 ∑ (xi − x) (xi − x)T . N − 1 i=1
(1)
Principal component analysis (PCA) was then applied to C yielding the modes of variation as the eigenvectors vi and their associated variances as the eigenvalues wi. A shape instance y can be approximated by: y = x + Vb,
(2)
where, V is the matrix with the first t modes of shape variations as columns and b is a weighting vector controlling the influence of each mode on the reconstructed shape (Cootes et al 1995). As illustrated in figure 3, individual bone SSMs for the humerus and scapula, along with a ‘combined’ SSM, constructed to address the relative spatial relationship between the humerus and scapula, were built in this study. 2.3.2. ASM-based segmentation. The automated bone segmentation was performed in a
three-level multi-resolution Gaussian pyramid (Fripp et al 2007). The ASM segmentation was performed by iterating the following two steps until the stopping criteria were satisfied (i.e. surface change was less than a threshold or the maximum number of iterations was reached): • Deformation. The surface was deformed by moving each vertex independently along its normal towards the mean gradient of the 1D intensity profile determined at that vertex along its normal. • Shape constraints. The pose and shape parameters, as encoded in b, of the deformed surface were estimated using SSM to regularize anatomically unlikely shapes. 1448
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
After initialization of the model parameters, a coarse-to-fine approach was taken to segment the proximal humerus and the scapula: • Coarse segmentation using the combined SSM which contained the spatial relationship between the two bones. This step was used for improving the atlas-based affine initialization. • Fine segmentation using the individual SSM of the humerus and the scapula. • Final segmentation involving bone surface relaxation (Fripp et al 2007) to fine-tune the bone boundaries by allowing the surface to free-form deform for a few iterations, not regularized by shape constraints but volume preserved smoothing. 2.3.3. BCI extraction. From the manually segmented cartilages, prior knowledge about the
BCI can be acquired by computing the probability p of each point on the bone surface being covered by cartilage. This probability was estimated as the frequency of cartilage tissue appearing above a given point on the bone surface in all of the training images. After bone segmentation, the BCI is automatically detected and extracted in a fashion of seeded region growing. The highly possible points on BCI (e.g. pi > 0.95, where i is the index to the points) were selected as seed BCI. Tissue intensity above the seed BCI were sampled along the normal to the surface and modelled as a Gaussian distribution. The distribution was used to grow the BCI in the following way. Every point with at least one neighbour on the current BCI was examined for the presence of cartilage tissue in the extracted intensity profile along its normal. If more than half of the profile is cartilage tissue, according to the estimated Gaussian distribution, this point was included as BCI. Estimation of tissue intensity distribution and BCI growing were repeated until convergence. 2.4. Validation method
To validate the developed automatic bone segmentation method for the shoulder, a leaveone-out strategy was employed. To segment the individual bones (humerus, scapula) on a given subject’s image, the manual segmentations on the remaining 24 datasets were used to build an SSM for segmentation. There were 25 runs of SSM construction and ASMbased segmentation. In each run, the segmentation accuracy was quantified using DSC values, mean absolute surface distance (MASD) and relative volume difference (RVD) measures (Xia et al 2013) with the manual segmentation as the reference. These quantities were computed before and after bone segmentation relaxation. Calculation of the DSC is defined as: DSC =
2TP , (2TP + FP + FN )
(3)
where, TP, FP, and FN are true positives, false positives and false negatives, respectively. They can be calculated using set operations on the automatic segmentation (A) and the reference (M). The DSC value indicates the relative true positive ratio with a larger value indicating better performance (1 is the maximum value). The MASD d is the average distance between points on two surfaces and defined as: d (A, M ) =
davg (SA, SM ) + davg (SM , SA) 2
1449
,
(4)
Z Yang et al
Phys. Med. Biol. 60 (2015) 1441
where, davg (SA, SM ) is the average directed surface distance from all points on surface SA to surface SM . The RVD is the volume difference in percentage of automatic segmentation A with respect to the reference volume M: RVD (A, M ) =
A − M × 100% . M
(5)
To evaluate any improvement in terms of volume overlap achieved by the relaxation step, paired t-tests were employed to examine for differences between the DSC values before and after relaxation. The segmentation accuracy of regions covered by the humeral and glenoid cartilages is of particular interest, given that the performance of the follow-up step of BCI extraction will be substantially affected by this step. To investigate the segmentation error distribution over the bone surfaces, with a particular focus on BCI regions, Hausdorff distance maps (Xia et al 2013) were generated for both the humerus and scapula. The accuracy of BCI extraction was evaluated using a surface equivalent to DSC values (Fripp et al 2007), measuring the agreement between the automatically and manually extracted BCI. To reduce the sampling error, the reference BCI was determined from the resampled manual cartilage segmentation using shape-based interpolation with a slice thickness of 0.125 mm. The TP was defined as the surface area where the point-to-point or point-to-surface distance between the extracted BCI to the reference BCI was less than 0.5 mm. This distance was calculated in both directions and the average was used. The FP region was defined as the surface area on the extracted BCI that was more than 0.5 mm away from any point on the reference BCI, while the FN region was defined as the surface area on the reference BCI that was more than 0.5 mm away from any points on the extracted BCI. Eight subjects were randomly selected and manual cartilage segmentations were performed by two expert analysts (MS and CE). After appropriate consensus to set rating parameters, the MR examinations were segmented by these experts in an independent fashion. Taking the first segmentation from expert MS as M, we simply replace A with BCI determined from a repeated manual segmentation performed by the same expert, or the other expert CE, for estimating intra- and inter-rater reliability, respectively. 3. Results 3.1. Statistical shape modelling
The SSMs built from the manual bone segmentations are shown in figure 3. Using the shrinkwrap surface extraction algorithm (Fripp et al 2007), segmented surfaces of the proximal humerus and the scapula were extracted from the 25 cases with each of the individual bone surfaces having 8,192 points and 16 380 triangles. The humeral, scapular and combined SSMs had different compactness (i.e. the ability to represent shape variation using as few modes as possible) and are plotted in figure 4 as the reconstruction precision against the number of modes. Here, reconstruction precision was defined as the cumulative relative shape variance in percentage terms. For example, the variation of shape accounted for by the first 5 modes was 72% for the humerus, 51% for the scapula and 78% for the combined humerus-scapula SSM. In this work, we used SSMs which accounted for greater than 95% of shape variations which equated to the first 16, 20 and 15 modes of the humeral, scapula and combined models, respectively. 1450
Z Yang et al
1
1
0.9
0.9
0.8
0.8
0.8
0.7
0.7
0.7
0.6 0.5 0.4
Precision
1 0.9
Precision
Precision
Phys. Med. Biol. 60 (2015) 1441
0.6 0.5 0.4
0.6 0.5 0.4
0.3
0.3
0.3
0.2
0.2
0.2
0.1
0.1
0 1
5
9
13
17
21
25
0 1
0.1 5
9
13
17
21
25
0 1
5
9
13
Mode
Mode
Mode
(a)
(b)
(c)
17
21
25
Figure 4. The compactness of each SSM constructed. (a) Humerus; (b) scapula; (c) combined.
3.2. Bone segmentation
An example (case #3) of the intermediate and final results from the automated segmentation of the humerus and scapula is given in figure 5. After the initialization and coarse segmentation stage, the placement of the combined ASM using the combined model is anatomically realistic, as shown in figure 5(b), although in the provided example, the joint space in the glenohumeral cartilage region appears overly narrow. This is most likely due to the spatial relationship encoded in the combined model not adequately expressing the variation present in this particular case. Following this step in the ASM process, the subsequent fitting of the individual humerus and scapula SSMs universally resulted in improved segmentation of these bones in relation to their corresponding anatomical profiles in the MR images for the case demonstrated in figure 5(c). Specifically, in this stage of the pipeline, the fitting of the individual humerus model (table 1) resulted in DSC values in the range of 0.684–0.947 with a mean of 0.909. In the more anatomically complex scapula, the range of DSC values was 0.549–0.825 with a mean of 0.723. The final relaxation step provided a localized fine-tuning of the bone boundaries to modify instances of likely under- or over-segmentation of the humerus and scapula, as demonstrated in figure 5(d). At the end of this pipeline, there were 23/25 and 5/25 cases having DSC values above 0.85 for the humerus (0.926 ± 0.050) and scapula (0.806 ± 0.054), respectively. In line with our future aim of developing a robust automated pipeline for cartilage segmentation in the glenohumeral joint, we calculated the DSC values for the scapula in an identified region of interest (ROI) centred around the glenoid fossa which encloses the BCI region of this bone, which reduced the effects of the large anatomical variance of the scapula outside this area, resulting in 15/25 cases having DSC values above 0.85 with an overall mean of 0.837. A boxplot showing the DSC values for different segmentation stages is given in figure 6. The improvements in terms of DSC values between the manual and automated segmentations of bone volume (table 1) achieved by relaxation were statistically significant for both the humerus (t = 12.104, p
