Medical Image Analysis 18 (2014) 567–578

Contents lists available at ScienceDirect

Medical Image Analysis journal homepage: www.elsevier.com/locate/media

Focused shape models for hip joint segmentation in 3D magnetic resonance images Shekhar S. Chandra a,1, Ying Xia a,⇑, Craig Engstrom b, Stuart Crozier c, Raphael Schwarz d, Jurgen Fripp a a

Australian e-Health Research Centre, CSIRO, Australia School of Human Movement Studies, The University of Queensland, Australia c School of Information Technology and Electrical Engineering, The University of Queensland, Australia d Siemens Healthcare, Erlangen, Germany b

a r t i c l e

i n f o

Article history: Received 31 July 2013 Received in revised form 29 January 2014 Accepted 5 February 2014 Available online 18 February 2014 Keywords: Hip joint WPCA Bone segmentation Shape models MRI

a b s t r a c t Deformable models incorporating shape priors have proved to be a successful approach in segmenting anatomical regions and specific structures in medical images. This paper introduces weighted shape priors for deformable models in the context of 3D magnetic resonance (MR) image segmentation of the bony elements of the human hip joint. The fully automated approach allows the focusing of the shape model energy to a priori selected anatomical structures or regions of clinical interest by preferentially ordering the shape representation (or eigen-modes) within this type of model to the highly weighted areas. This focused shape model improves accuracy of the shape constraints in those regions compared to standard approaches. The proposed method achieved femoral head and acetabular bone segmentation mean absolute surface distance errors of 0:55  0:18 mm and 0:75  0:20 mm respectively in 35 3D unilateral MR datasets from 25 subjects acquired at 3T with different limited field of views for individual bony components of the hip joint. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Active Shape Models (ASMs) developed by Cootes et al. (1995) (also referred to as deformable models) have been successfully applied to segmenting three dimensional (3D) MR images of musculoskeletal structures (Neubert et al., 2012; Schmid et al., 2011; Fripp et al., 2010) as well as for other structures such as the prostate (Martin et al., 2010; Chandra et al., 2012), liver (Heimann et al., 2006) and the heart (Zheng et al., 2008; Ecabert et al., 2008). In this approach, a triangulated surface is deformed to fit an object of interest while simultaneously utilising other image features and/ or priors, such as the object shape and appearance. For example, Fripp et al. (2010) used deformable models to (fully) automatically segment the articulating bone elements of the knee joint (femur, tibia and patella) from 3D T2-weighted water-excited Double-Echo Steady State (weDESS) MR images utilising shape priors for each bone. The bone segmentations were then used to accurately localise the bone-cartilage interfaces (BCIs) from which individual joint cartilage plates can be segmented and cartilage thickness ⇑ Corresponding author. Tel.: +61 732533660. E-mail addresses: [email protected] (S.S. Chandra), [email protected] (Y. Xia). 1 Principal corresponding author. http://dx.doi.org/10.1016/j.media.2014.02.002 1361-8415/Ó 2014 Elsevier B.V. All rights reserved.

determined for early Osteoarthritis (OA) assessment (Altman et al., 2004). Accurate 3D segmentations of the osseous structures of the hip joint are important as current OA assessments utilise plain radiographs that are inconsistent in the early disease stages as they rely on bony features (abnormalities) that typically occur in later pathological stages of OA (Pollard et al., 2008). Most of the previous bone segmentation work has been performed using computerised tomography (CT) (Lamecker et al., 2004; Seim et al., 2008; Audenaert et al., 2011; Masjedi et al., 2012), plain radiographs (Ding et al., 2007; Dong et al., 2007; Zheng et al., 2010; Baka et al., 2011) or ultrasound images (Barratt et al., 2008). MR imaging enables noninvasive, 3D assessment of the joint structure including biochemical changes with no ionizing radiation and excellent soft tissue contrast (Burstein et al., 2000). This enables the simultaneous visualisation of both the bone and cartilage tissue as opposed to other imaging modalities. Segmentation of these osseocartilaginous structures of the hip joint is however more challenging in 3D MR images due to the presence of partial bone coverage resulting from the large joint geometry, relatively thin cartilage plates that are present (Hodler et al., 1992) and the partial voluming errors caused by its spherical structure (Naish et al., 2006). Schmid et al. (2011) developed a robust shape model (RSM) approach incorporating the robust Principal Component Analysis

568

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

(PCA) of Skocˇaj et al. (2007) to segment pelvic and femoral bone elements, visualised as partial structures, in small field of view (FoV) MR images. Whilst they reported encouraging results for bone segmentation with an average distance error of 1:12  0:46 mm, further improvements in these bone segmentations, especially for both the femoral and acetabular BCI regions, would likely be of benefit for subsequent analyses of the relatively thin (1–2 mm thick) hip cartilages (Shepherd and Seedhom, 1999). This improved accuracy would also be particularly important in subsequent segmentation of the separated cartilage plates. Other techniques for segmenting bones, such as utilising the image information in the scans directly, include the works of Dalvi et al. (2007), Zoroofi et al. (2004), Nguyen et al. (2007) and Bourgeat et al. (2007). In this work, we introduce a weighted shape learning approach (via the spatially weighted PCA of Thomaz et al. (2010)) for regionspecific accuracy and compactness applied to the human hip joint (see Fig. 1). Predetermined weights may be set to each corresponding point of an anatomically selected region in the shape model allowing the focusing of the shape constraint energy on areas deemed most important from a research or clinical perspective. This focused shape model (FSM) orders the shape representation to the up-weighted regions, in contrast to other approaches, where the representation is ordered from the largest variations to the smallest, regardless of the regions that might be of more importance clinically. In the FSM, the variation in the highly weighted areas become the dominant representations within the model and the resulting eigen-modes are ordered by importance to the problem at hand. The predicted outcomes from our work with the FSM approach are that it provides: 1. A more compact shape representation which requires substantially fewer modes to achieve equivalent accuracy to the RSM proposed by Schmid et al. (2011) for small FoV MR images when focusing on the hip joint.

2. Lower reconstruction errors and higher accuracy when using the same precision of shape representation for the hip joint when compared to the RSM. The paper is structured as follows. In the next section, the background information for this work is covered. Section 3 introduces the proposed FSM and Sections 4 and 5 present the results and discussion for its application to the hip joint respectively. The experiments conducted compare the FSM to the work of Schmid et al. (2011) on small FoV MR image bone segmentation with the eventual goal of utilising these as the basis for segmenting the femoral head and acetabular cartilages of the hip joint in future work. 2. Background In medical imaging, it is important to keep the shape of the object of interest diagnostically interpretable so that it can be useful to a clinician. Typically at the deformation stage of the ASM, a surface S (having vertices v i with 0 6 i < N, where N is the total number of vertices) is deformed in a direction that is normal to the surface at each vertex, limited only by an allowed maximum length of displacement m and independently of its neighbouring vertices. Neighbouring vertices can be accounted for by either B-spline fields (Martin et al., 2010) or surface smoothing algorithms (Chandra et al., 2012), such as those of Taubin et al. (1996). However, even with these neighbourhood constraints, the surface can still diverge from the realistic (expert visualised) anatomic representation of the object. 2.1. Shape models in image segmentation Cootes et al. (1995) demonstrated that a Statistical Shape Model (SSM) can be built from a set of landmark points (also called a point distribution model), so long as the points have correspondence, i.e. each point on the surfaces represent the same part of the anatomy

Fig. 1. A focused approach to shape models applied to the bones of the hip joint (femur bottom and pelvis top), where the shape representation (shown as scalar and vector fields of the primary mode (unit-scale) variations from the mean shape by 3r) can be ordered by a priori selected anatomical structures in contrast to a robust model. These structures are produced by manually contouring regions of a template surface (femoral head, femoral bone-cartilage interface, fovea, acetabulum, acetabular bone-cartilage interface and fossa) and Gaussian weighting their distance maps at the label interfaces to reduce discontinuities within the model. The rest of the surface is retained in the segmentation process by down-weighting their contribution in the model to still allow their use in pose adjustments during the segmentation process.

569

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

and are same in number over the entire surface. The most popular approach in building these shape models is to utilise the PCA (Pearson, 1901), which requires a rectangular (3N  s) data matrix X ¼ xj to determine the eigen-decomposition (see Algorithm 1). Thus, meaningful variations are obtained for each row in X that represents the same regions of the surfaces S j .

2.2. Robust shape models The RSM approach utilises the robust PCA of Skocˇaj et al. (2007) that can handle missing or corrupted data (hereon referred to as missing data) using binary weights wij to represent missing or known data for each point in each training surface. In this

Algorithm 1. PCA Algorithm (Pearson, 1901, adapted) 1: X 2: l b 3: X

Sj mean shape from X. . Each column consists of variational points .Inner product is to reduce matrix size

Xl

bTX b. 4: Compute inner product C ¼ ð1=sÞ X b to get 5: Compute the singular value decomposition (SVD) of C and re-project onto X estimate of U. 6: return U and ki

The PCA produces an eigen-decomposition of the covariance b T X, b where X b represents the training shapes with matrix C ¼ ð1=sÞ X the mean shape l removed. The eigen-decomposition results in a set of eigenvectors U ¼ ½u1 ; . . . ; uN  ¼ ui and eigenvalues ki that effectively represent a set of modes that model how each vertex varies from the mean vertex given the training surfaces. One can then project a shape S ¼ ½v 1 ; . . . ; v N  ¼ v i as

bi ¼ U T ðv i  li Þ;

ð1Þ

into the variation (PCA) model to acquire the variation coefficients ^i bi . The shape S can then be reconstructed as a new shape b S ¼v within the constraints of the training set as

v^ i ¼

n X

bji uji þ li  v i ;

ð2Þ

j

where n is the number of modes, so that n 6 s. In general, only the most significant modes are used, up to a precision of representation P : 0 < P 6 1, for computational efficient and stable reconstructions. For a more detailed review of ASMs see Heimann and Meinzer (2009). Most ASM based works have utilised the PCA for shape learning until Schmid et al. (2011) introduced robust shape models (RSMs).

. Re-projection undoes the matrix inner product . The shape model

robust PCA, the variations are modelled from the mean of the known data as

1

li ¼ PN j

wij

N X wij xij :

ð3Þ

j

Any detected outliers are replaced by iteratively reconstructed values from a model of the reliable or known data using an Expectation Maximisation (EM) algorithm (see Algorithm 2). The final model then consists of known data and reconstructed data. The EM steps of the iterative reconstruction algorithm are used to recover missing data during the projection and reconstruction stage. For a more detailed discussion of the method see Algorithm 4 of Skocˇaj et al. (2007). Schmid et al. (2011) were motivated to use the robust PCA as Skocˇaj et al. (2007) had shown that it outperformed the standard PCA in the presence of missing data. Schmid et al. (2011) then applied the RSM approach to the special case of small FoV MR images, so that the parts of the surfaces outside the image were ignored during the training and ASM fit. The iterative reconstruction algorithm of Skocˇaj et al. (2007) was used to train the RSM with partial surfaces by replacing unknown parts with reconstructed ones. Then, at each iteration of the ASM fit, the vertices of the surface were identified as being outside or inside the image using binary valued weights. To apply the shape constraints, the surface was initially projected as Eq. (1) and then

Algorithm 2. Iterative Reconstruction PCA Algorithm (Skocˇaj et al., 2007) 1: X 2: W 3: l b 4: X

S j with known data. binary weights representing missing and known data. weighted mean from known data from Eq. (3).

. Can accommodate partial surfaces in training

Xl

b via Algorithm 1. 5: Apply PCA on known data in X b X 6: Y 7: while not converged do 8: E-step: Project known data with the pseudo-inverse of weighted U. b with current model. 9: M-step 1: Reconstruct X b with missing data replaced by reconstructed values. 10: M-step 2: Y X 11: M-step 3: Re-estimate U by PCA on Y. 12: end while 13: return U and ki

. Initialise the algorithm with the known model

. Element-wise product with W . Recovers missing data here

. The robust shape model

570

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

the (binary) weighted Procrustes alignment algorithm used to ensure the accurate reconstruction of the known parts of the projected surface (those inside the image). This improves accuracy and convergence of the segmentation in the presence of a small FoV. Being able to process small FoV is important for the hip joint as the cartilage is relatively thin compared to other joints such as the knee (Hodler et al., 1992), requiring MR images to have higher resolutions. We show that the work of Schmid et al. (2011) on the hip joint can be improved to incorporate weighting of the shape model regions during both the learning stage and in the application of the shape constraints. The weighting not only allows one to handle small FoV as well, but also enables the focusing of the shape constraint energy into the specific areas deemed most important during the segmentation (e.g. clinically) for more accurate segmentations. 3. Focused shape models In a number of applications it is necessary to get a more accurate segmentation of a specific region (s) of an object. For example, our eventual goal is to segment the hip cartilages for which some methods rely on an accurate bone segmentation near the BCI as priors (Fripp et al., 2010). By utilising the spatially weighted PCA, the BCI area of the hip joint can be up-weighted to account for important areas, all the while utilising the other areas to ensure proper pose alignment and more accurate mapping to learned shape subspace that the regions of interest (ROI) may be dependent on, such as its adjacent areas (see Fig. 1). 3.1. Weighted shape learning The spatially weighted PCA of Thomaz et al. (2010)2 is simpler to implement than the robust PCA as can be seen from comparing Algorithms 2 and 3.

bi ¼ V y wi ðv i  li Þ;

ð5Þ

where V represents the row-weighted eigenvectors U and y represents the pseudo-inverse of V, which effectively solves the over-determined set of linear equations in a least squares sense (Skocˇaj et al., 2007). To reduce computation complexity in the projection stage, only the first n 6 s need to be used. Lastly, the weights (4) are also used to improve accuracy during the optimisation of the shape model pose at the projection stage by using only the upweighted areas during pose alignment. This is achieved by cropping the surfaces to the same up-weighted regions and computing the registration (Gower, 1975) and applying the resulting transform to the entire surface. The spatially weighted PCA of Thomaz et al. (2010) has important distinctions to the temporally and spatially weighted PCA proposed by Skocˇaj et al. (2007). The temporally weighted PCA weights each training shape (or column of the training matrix), which is suitable for handling unreliable or low quality training shapes as a whole. The spatially weighted PCA of Skocˇaj et al. (2007) is a general scheme solved using an EM algorithm for weighting arbitrary points in any of the training shapes and not just each corresponding point (or row of the training matrix) as with the spatially weighted PCA of Thomaz et al. (2010), so that w is a matrix and not a vector. In the spatially weighted PCA of Thomaz et al. (2010), entire regions maybe weighted and their variations from the mean shape scaled as desired. This weighting of a corresponding point for all training surfaces such as (4) is not possible with the spatially weighted PCA of Skocˇaj et al. (2007) as the weightings are lost from the eigenvectors in the process of removing the weighted mean from the data matrix3. The authors also found that the eigenvalues of this spatially weighted PCA was difficult to accurately initialise for reconstructing bone surfaces. Skocˇaj et al. (2007) suggested using random values, which in the authors experience, led to convergence towards very sub-optimal local minima.

Algorithm 3. Spatially weighted PCA algorithm (Thomaz et al., 2010) 1: X 2: wi 3: l b 4: X

Sj . square root of the weights representing ROIs. weighted mean from Eq. (3).

. Weights represent important regions

Xl . Weight each corresponding point

b with wi to get Y. 5: Multiply each row of X 6: Apply PCA on Y via Algorithm 1. 7: return U and ki

. The focused shape model

The weights

wi ¼

pffiffiffiffi  pffiffiffiffi pffiffiffiffi w1 ; w2 ; . . . ; wN ;

ð4Þ

are determined from the surface regions that were contoured (see Fig. 1) and are applied to all surfaces via the point-to-point correspondences present. Smooth transitions were created at the contour interfaces by using signed distance maps computed from the region boundaries on the template surface. This was done in order to avoid discontinuities in the shape constraints, such as steps caused by dramatic weight changes. The resulting weights incorporating these smooth transitions were created by superimposing each region (largest in spatial extent to smallest) after Gaussian weighting the positive values of the signed distance maps that represent the distance from the region boundaries. The vertices of the surface outside the image were handled by providing weights during the projection stage as 2

See also the related work of Thomaz and Giraldi (2010).

The last consideration when using the weighted PCA is that it cannot accommodate missing data in the same way as the robust PCA. This means that zero-valued weights cannot be used, but can be handled by using very small values (e.g. 1  106 ) instead. These weight values can be used to denote the vertices of the surface outside the image in both the training stage or the projection stage using Eq. (5). 3.2. Method For this work, a template surface was manually contoured to mark all the ROIs for the hip joint by S. Chandra under expert guidance from C. Engstrom. The regions included the up-weighting of the overall likely BCI region, the acetabulum (and surrounding 3 Eqs. (7) and (8) of the work of Skocˇaj et al. (2007) separates the weights from the eigenvectors and so no weights (such as (4)) are directly encoded in the resulting eigenvectors leading to no weighted dimensionality reduction normally afforded to us by the PCA.

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

sacrum), the femoral head and down-weighting the highly variable but areas identified as less important, such as the iliac crest and the non-cartilagenous fovea and fossa areas that contain ligaments that connect the pelvis to the femur. A total of 56 full FoV manual contours of CT scans from Calvary Mater Newcastle Hospital (Dowling et al., 2012) were utilised in obtaining training surfaces with correspondences via a simple template surface optimisation technique developed for this work (see Appendix and Algorithm 5). The same 56 training surfaces were utilised in building both the RSM and FSM. The weights illustrated in Fig. 1 were then used to build two FSM of the hip joint. The first FSM was built directly from the weights in Fig. 1 taking into account each sub-region of the joint and BCI (hereon referred to as the ‘‘Focused BCI’’ model), while down-weighting the remaining parts of the bones to 0.1. The second FSM was designed for the entire joint area with equal weighting for all constituents of the joint (hereon referred to as the ‘‘Focused Joint’’ model) and downweighting the remaining parts of the bones to 0.1. The weights in the latter model were designed to encompass the immediately adjacent areas such as the sacrum and labrum (rim of the acetabulum) and the full femoral head to better fit the BCI areas of the joint. The primary mode of the Focused BCI model for each bone is also shown in Fig. 1. Once the FSM were obtained, the bone segmentation pipeline described in Algorithm 4 was applied to 3D MR images.

571

Fig. 2. The atlas image (as coronal and sagittal views) and the template mesh used for initialising the deformable model utilised in this work. The green box shows the region used to mask the segmentations during the validation of the approach. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

registrations were also quality controlled by checking for exceptions during the registration process thrown by the metric when the FoV between the images caused problems. For these exceptions, the femoral head locator was used as the initialisation for the surface.

Algorithm 4. Hip segmentation pipeline 1: input case image, atlas image, atlas surfaces (combined bone, pelvis, femur) and their shape models. 2: Initialise atlas via a Hough transform (Duda and Hart, 1972) based femoral head locator with the radius as the scale (Nishii et al., 2004). 3: Compute the rigid registration between the atlas image and the case image (Ibanez et al., 2005ITK Toolkit). 4: Propagate the atlas surfaces representing each bone to the new case using the transforms computed. 5: for each of the combined, pelvis and femoral bone models do . Combined model for small bone pose adjustments after initialisation 6: for two image pyramid levels do 7: for each iteration up to maximum iterations i do 8: Freely deform the surface (i.e. without shape constraints) toward the image gradient along the vertex normals. . Up to m 9: 10: 11: 12: 13: 14: 15: 16: 17: 18:

Identify regions of the surface inside the image and weight appropriately. . Robust: 0, Focused: 1  106 Use the surface smoothing algorithm of Taubin et al. (1996) to smooth the deformations. Use the relevant shape model (either robust or focused) to reconstruct the surface within the shape constraints. end for end for if fitting combined model, initialise the remaining bones with the combined bone surface. .Improves the individual bone poses end for Allow the surface to freely deform along the vertex normals for a few iterations to get closer fit. Smooth surface as before to get final bone surface (s). return Surfaces that represent the bone segmentations.

In this pipeline, the template surface was initialised into the new scans using a femoral head locator based on detecting a sphere via the Hough transform (Duda and Hart, 1972) as first shown by Nishii et al. (2004). This initialisation was further refined by using a Mattes Normalised Mutual Information (Mattes et al., 2003) rigid registration (Ibanez et al., 2005, ITK Toolkit) of the new scan (with scaling provided by the femoral head detector) to the template surface’s corresponding atlas image (see Fig. 2). The atlas image, which was constructed from 10 weDESS images having a similar FoV using the non-rigid diffeomorphic demons algorithm of Vercauteren et al. (2009) after affine initialisations, was resampled into the new patient space using the femoral head locator transform parameters before each registration. The

Once the surface was initialised, it was fit to the image using gradient information as described in Algorithm 4 that incorporates a hierarchical object fitting with combined and individual bones (Yokota et al., 2009). The combined bone fit provides an initial rotational estimate of the pose while the individual bone refinements correct the flexion pose of the bones. The resulting surface was voxelised to produce a binary image and validated with manual delineations of 25 healthy volunteers. 3.3. Validation The FSM and the bone segmentation pipeline in Algorithm 4 were validated against 35 3D MR image datasets of 25

572

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

Table 1 The acquisition parameters for each MR sequence. Parameter

weDESS

truFISP

Plane Spacing (mm) Slice Thickness (mm) Repetition Time (ms) Echo Time (ms) Bandwidth (Hz) Flip Angle (degrees)

0.67 0.7 15.46 5.16 202 25

0.23 0.49 10.65 4.46 140 30

healthy volunteers acquired using a 3T Siemens Trio Scanner. The volunteers had activity history levels ranging from normal everyday activities to high level participation in sports such as water polo and rugby. A large 4-channel body matrix was used during the MR acquisition with Generalized Autocalibrating Partially Parallel Acquisition (GRAPPA) (2) enabled. The medical research ethics committee of the University of Queensland approved the study and informed written consent was obtained from all participants involved the research. The MR examinations consisted of 25 bilateral weDESS and 10 unilateral T2-weighted True Fast Imaging with Steady-state Precession (truFISP) sequences. A summary of the main acquisition parameters used in this work is provided in Table 1. The weDESS sequences were manually segmented by Y. Xia using ITK-SNAP (Yushkevich et al., 2006) under expert guidance from C. Engstrom with each case taking several hours to complete. These manuals were rigidly co-registered to their corresponding truFISP images using ITK (Ibanez et al., 2005). A simple automatic cropping was performed on the bilateral weDESS images by bisecting the x-axis to produce the right side of the hip in order to approximate unilateral scans. Both MR sequences had partial FoV, each with different coverage for the bones (see Fig. 3). All images were pre-processed by correcting for bias fields via N4 (Tustison et al., 2010), statistically normalising (zero mean and unit variance) and median smoothing them to reduce the effects of noise and image artefacts. The segmentation performance was analysed using mean absolute surface distance (MASD) near the BCIs of the hip joint and Dice’s similarity coefficient (DSC) (Dice, 1945) scores within the masked joint area highlighted in Fig. 2 by propagating the joint mask region from atlas space via the initialisation transform. 4. Results To assess the performance of the FSMs, a comparison was made to an in-house implementation of the RSM proposed by Schmid

Fig. 4. The reconstruction errors (in terms of MASD) of the robust and focused shape models for the (a) pelvic and (b) femoral bones of the hip joint up to 95% shape representation precision. The focused model of the joint area retains an accuracy improvement of about 0.2 mm over the RSM (which is quite large when compared to the hip cartilages which can be of the order of 1 mm in thickness) regardless of the number of modes utilised in both models.

et al. (2011) (with one minor improvement) in three experiments with respect to the hip joint (noting that the eventual aim is to segment the hip joint cartilages in future work). Firstly, the generalisability (i.e. the reconstruction errors of the training surfaces using a leave-one-out approach) and relative compactness (i.e. the dimensional reduction of the shape representation) of the models were tested and compared within the BCIs of the hip joint. Secondly, the FSMs were applied to hip joint segmentation and the results validated against manual analyses as described in Section 3.3. Thirdly, the performance of a region-based RSM with respect to the FSMs were examined to evaluate down-weighting areas in the FSMs with regards to segmentation accuracy for the selected anatomical regions. The minor improvement to the proposed method of Schmid et al. (2011) involved using the RSM in the projection stage as well as the weighted Procrustes alignment to improve finding the closest shape representation in the shape model. 4.1. Generalisability

Fig. 3. Examples of the MR images and their FoV shown as coronal and sagittal slices within a 3D rendering of the propagated atlas surface utilised for this work. (a) shows an example of a bilateral weDESS image (0:67  0:67  0:6 mm, cropped to one side) and (b) shows an example of a unilateral truFISP image (0:23  0:23  0:49 mm). Black regions around the truFISP scan are not part of the actual volume.

The Focused Joint model weighted within the entire joint region was found to be always more accurate than the RSM for both the pelvic and femoral bony components with an improved accuracy of around 0.2 mm on average regardless of the number of modes used within the RSM (see Fig. 4). When using 90% shape representation for each bone with the Focused Joint model and using the same number of modes for the RSM, the improved accuracy was found to be 0:24  0:12 mm; this is salient given that hip joint cartilage thickness is of the order of 1 mm. The Focused BCI model

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

573

Fig. 7. The femoral validation MASD errors (in mm) for 25 healthy subjects in the joint area. The statistics were computed from the automatic segmentations of 25 weDESS and 10 truFISP unilateral MR examinations using Algorithm 4.

Fig. 5. The compactness of the robust and focused shape models for the (a) pelvic and (b) femoral bones of the hip joint. The focused model compacts better due to the most important information being preferentially up-weighted a priori.

ordering of the primary eigen-modes to best suit the focused regions, which in turn creates a more compact shape model due to the a priori weighted regions. 4.2. Compactness The compactness of a model measures the dimensionality reduction of the most important information perceived by the model and are shown as plots for the pelvic and femoral bone models for the RSM and FSMs in Fig. 5. The improved compactness of the Focused approach was found to be most apparent when using only 3 modes for both the femur and acetabulum models in the segmentation process. This compact FSM (hereon referred to as ‘‘Focused Compact’’ in subsequent tables) provided equivalent reconstruction accuracy as the RSM that utilised 15 and 7 modes for the pelvis and femur respectively. In the next section, we show that the accuracy improvements shown by the generalisability and compactness are translated into practical improvements in accuracy of the bone segmentations of the hip joint. 4.3. Quantitative analysis

Fig. 6. The acetabular validation MASD errors (in mm) for 25 healthy subjects in the joint area. The statistics were computed from the automatic segmentations of 25 weDESS and 10 truFISP unilateral MR examinations using Algorithm 4.

was found to have similar improved accuracy, but reverted to an equivalent accuracy of the RSM when using a very large proportion of shape representation within the RSM. This is consistent with the

Figs. 6 and 7 show performance of the FSMs as the surface distance errors between the manual and the bone segmentations for all images within the BCI areas of the hip joint. The acetabulum had a surface error of 0:75  0:20 mm in the joint area, which was 0.35 mm lower than the RSM and was found to be statistically significant (p < 2  106 ). The femur showed a less improvement in the surface error by 0.15 mm than the RSM (p < 2  106 ) having a surface distance of 0:55  0:18 mm. The validation results for the masked weDESS and truFISP segmentations are shown in Tables 2 and 3 respectively. Representative examples of the bone segmentation results using the Focused Joint model are shown in Fig. 8 for three truFISP images. A table of the segmentation pipeline parameters used in this work are given in Table 4. The average total computation time on a 3.2 GHz Intel Xeon desktop PC was found to

574

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

Table 2 Validation results for the weDESS unilateral MR dataset within the masked region (given in Fig. 2) of the joint area.

Robust Focused Joint Focused BCI Focused Compact

Bone

DSC mean

DSC median

DSC Max

Hausdorff distance mean (mm)

Hausdorff distance median (mm)

MASD (mm)

MASDMedian (mm)

Pelvis Femur Pelvis Femur Pelvis Femur Pelvis Femur

0.94 0.98 0.95 0.98 0.95 0.98 0.94 0.98

0.94 0.98 0.95 0.98 0.95 0.98 0.94 0.98

0.97 0.99 0.97 0.99 0.97 0.99 0.96 0.99

7.72 8.53 7.75 9.05 7.66 9.92 7.91 10.17

7.37 5.21 7.67 5.36 7.43 5.41 7.80 7.99

0.71 0.34 0.69 0.33 0.68 0.35 0.78 0.42

0.66 0.29 0.67 0.29 0.66 0.28 0.73 0.34

Table 3 Validation results for the truFISP unilateral MR dataset within the masked region (given in Fig. 2) of the joint area.

Robust Focused Joint Focused BCI Focused Compact

Bone

DSC mean

DSC median

DSC max

Hausdorff distance mean (mm)

Hausdorff distance median (mm)

MASD (mm)

MASDMedian (mm)

Pelvis Femur Pelvis Femur Pelvis Femur Pelvis Femur

0.92 0.96 0.92 0.97 0.91 0.96 0.91 0.96

0.92 0.96 0.91 0.97 0.91 0.97 0.91 0.96

0.96 0.98 0.96 0.98 0.95 0.98 0.95 0.98

7.24 9.92 7.29 9.68 7.86 10.18 7.80 10.74

7.33 5.97 6.84 5.40 7.97 7.81 7.94 8.74

0.95 0.55 0.99 0.47 1.06 0.51 1.08 0.57

0.95 0.45 1.05 0.43 1.07 0.44 1.14 0.47

be 6.0 ± 0.9 min with 1.3 ± 0.5 min for pre-processing, 3.1 ± 0.3 min for the bone segmentation and the remaining time taken up by the initialisation and free deformation stages of the pipeline. 4.4. Weights analysis Finally, an example of the importance of down-weighting less important bone areas is shown in Fig. 9, where a region-based RSM is utilised to simulate the Focused Joint model using the robust PCA, so that only the joint area was maintained. It was found that this type of model does not perform as well as the standard RSM approach. It had a weDESS (mean, median) DSC of ð0:93  0:03; 0:93Þ and ð0:97  0:01; 0:98Þ, and truFISP (mean, median) DSC of ð0:90  0:03; 0:89 and ð0:95  0:02; 0:96Þ for the pelvis and femur respectively. A similar effect was observed in the FSM approach (not shown) when the weight of the other parts of the bone was further down-weighted from a value of 0.1 to a value of 1  106 . 5. Discussion The FSMs were found to have a highly compact shape representation that typically required only a few modes to achieve the equivalent accuracy as the RSM for segmentation of the bony components of the hip joint in small FoV MR images. Alternatively, lower reconstruction errors and higher accuracy were obtained for bone segmentations when using the same precision of shape representation as the RSM. This improvement in the BCIs of the hip joint (around 0.2–0.4 mm on average), most notably for the acetabulum, could be beneficial for ensuring successful clinical adoption of bone segmentation algorithms for early OA assessment, since the hip cartilages can be of the order of 1 mm in thickness. 5.1. Generalisability Fig. 6 shows that the FSMs have improvements in the rim of the acetabulum, an important area in measuring femoroacetabular impingement, a condition implicated as a precursor in the development of early OA (Ganz et al., 2003). The smaller improvement in the

femur was likely due to both approaches achieving very good accuracy (mean DSC of 0.98) and always having a larger FoV than the acetabulum in most truFISP images. The better accuracy of the Focused Joint model over the Focused BCI model in the BCI region shows that the region is dependent on the adjacent bony areas for localisation and pose, such as the rim of the acetabulum, and the sacrum. There was no significant difference between the RSM or the compact FSM (having three modes) during segmentations of either bone. 5.2. Quantitative analysis The validation results between the RSM and Focused Joint model in Tables 2 and 3 show that the Focus Joint model compares favourably to the RSM for the weDESS pelvis and truFISP femur results despite the inclusion of some regions of the bone that are outside the focusing area of the FSMs due to the rectangular nature of the validation mask. These include areas such as the bone adjacent to the sacrum, opposing-side of the pelvis from the acetabulum, the femoral neck and trochanter region. In terms of MASD validation in the joint area, the Focused Joint model accuracy over the RSM was found to be significantly higher (p < 2  106 ) for both bones (see Figs. 6 and 7). Overall the approximate improvement is 0.2–0.4 mm in the joint area for both bones, which could be clinically important in imaging analyses of the hip cartilage given their thickness is of the order of 1 mm. The in-house implementation of the RSM (having minor improvements in the form of robust shape constraints) also had MASDs similar to the work of Schmid et al. (2011) on the hip, which achieved a MASD of approximately 1.1 mm for the bones. The Focused Joint model MASDs also compare well with respect to their work. 5.3. Compactness The FSMs were also found to preferentially order the eigenmodes of the shape representation to the up-weighted regions. The eigen-decomposition process generally orders the resultant eigenvectors or eigen-modes based on the magnitude of their representation within the model, i.e. the largest variations have the largest eigenvalues. In the RSM, the modes are thus represented as the

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

575

Fig. 9. An example of how the down-weighting of the non-joint region of the hip joint, as opposed to the binary nature of the weights in the RSM, can be beneficial to the segmentation process. (a) shows the RSM approach for a unilateral weDESS image shown in coronal view with the contour representing the segmentation achieved by this method. (b) shows the result of a region-based RSM approach to simulate the Focused Joint model with the robust PCA, i.e. areas that are weighted less than unity (see Fig. 1) are always marked as unknown. (c) shows the result of the Focused Joint model.

the primary eigen-mode of the FSM, in contrast to the RSM. This could be readily used in optimisers that incorporate a few shape model parameters in addition to other parameters during the segmentation process, which would improve the convergence and performance of the process and ensure subtle modes (responsible for ROI variations) are not lost either by being present in the top 1–3% of the shape representation or being misclassified as noise. An example of this can be seen Fig. 4(a) where reduced reconstruction error occurs in the RSM by the inclusion of a few extra modes in addition to 20 most significant modes of the model. 5.4. Weights analysis Finally, the performance of the region-based RSM was less favourable when compared to the RSM, especially for the pelvis, and was likely due to the same reason. Both the RSM and the region-based RSM have the same modes present in each model and the only difference is the amount of missing data during the projection stage. In contrast, the modes in the FSMs are ordered according to the weights including any subtle modes that might be important (see Fig. 4(a)). 5.5. Summary In summary, the advantages of the FSM are as follows: Fig. 8. Three examples of the Focused Joint model bone segmentation results was overlaid on the original truFISP MR images in coronal view. (a)–(c) have DSCs of (0.97, 0.91), (0.97, 0.94) and (0.98, 0.94) for the (femur, acetabulum) within the joint region respectively, where each tile shows a different slice within the volume. Table 4 The parameters for the FSM approach. Parameter

Value

Profile spacing h Profile length ‘ Smoothing iterations R Maximum iterations i Maximum displacement m Vertices in model (femur, pelvis) N

1/2 Image spacing 24 5 24 2 mm (3309, 7281)

largest variations in the training shapes, regardless of the regions that might be of more importance clinically (see Fig. 5). In the FSMs, the variation in the highly weighted areas become the dominant representations within the model and the resulting eigen-modes are ordered by importance to the problem at hand. This is shown in Fig. 1, where the highly weighted BCI region is represented by

1. Higher accuracy in the joint area – When using an equivalent representation as the RSM (in terms of the precision of the model), a higher accuracy in terms of reconstruction and validation MASD errors are achieved (see plots of Fig. 4 and error statistics of Figs. 6 and 7). 2. Focused primary modes – The modes match the selected anatomical structures ensuring modes responsible for subtle variations in those regions are not lost (see Figs. 1 and 4). This could have future applications where primary eigen-modes are used for detailed study of object shapes. 3. More compact representations – The FSMs inherently compress more effectively owing to the weighting of important information a priori (see plots of Fig. 5 and the validation of the focused compact model in Tables 2 and 3). 4. Arbitrary weights in model – The FSMs allow weights for regions and their surrounding areas, as opposed to binary weights with the RSM, i.e. either missing or known (see Fig. 9 for example). There are two readily identified weaknesses with the FSMs when compared to the RSM. Firstly, it cannot handle surfaces with

576

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

arbitrary missing vertices (i.e. partial surfaces) during the training of the shape model without introducing errors into the model. The FSMs also revert to the mean shape as the weights are decreased, which is not the same as the RSM for missing parts that effectively interpolates the unknown data from the model of known data. Thus, the RSM is better suited for applications where the recovered missing parts are important clinically or when no imaging data is available (see for example Fig. 9), such as attempting to recover the unknown anatomical areas when only the other parts of the anatomy are known. Further work needs to be done in applying the FSM approach in segmenting the cartilages of the hip joint and to study pathological hip imaging data. Pathological imaging data is challenging because of the various joint abnormalities of OA, such as osteophytes and bone lesions, and are highly patient dependent. The present method will likely experience segmentation difficulties for severe joint abnormalities, leading also to decreased segmentation accuracy for the remaining healthy bone components. Specialised methods for lesions, such as the recent work of Dodin et al. (2012), could prove valuable. The development of a robust outlier detection scheme to utilise the weightings available from the focused models could also prove to be important in this area.

6. Conclusion This paper presented a weighted shape learning approach for deformable models applied to hip joint segmentation in 3D MR images. A spatially weighted PCA was used to build a focused shape model (FSM) that allowed selective up-weighting of the hip joint area to create a region-sensitive and region-compact shape representation for the joint. The up-weighting also resulted in an ordering of the model primary modes that was more related to the local anatomy than using standard approaches. The focused model achieved the same accuracy (in terms of generalisiability and surface distance error to manual contours) as the RSM when utilising only three modes for shape representation for each of the bones (see Tables 2 and 3). An approximate improvement of 0.2–0.4 mm for surface distance errors of both femoral head and acetabular bone regions was obtained when using the same shape representation as the robust model (see Figs. 6 and 7). This improvement could prove important clinically in imaging analyses of the hip cartilage, which can be of the order of 1 mm in thickness.

Acknowledgments This research was supported under Australian Research Council’s Linkage Projects funding scheme LP100200422 and partially funded by the Cancer Council NSW (Project Grants 07-06 and 1105). Appendix A A.1. Template Surface Optimisation Obtaining a set of surfaces with point correspondence can be very challenging and is itself a vigorous area of research. Some popular methods involve spherical parameterisations of the surfaces, such as using Spherical Harmonics Description (SPHARM) (Brechbühler et al., 1995) and then the shape-image technique of Davies et al. (2008). This approach is only valid for genus zero objects and has been successfully used in a number of applications (Chandra et al., 2012; Shen et al., 2012; Neubert et al., 2012). Another popular approach (not dependent on surface genii) utilises surface registration methods to determine correspondences, such as the work of Hufnagel et al. (2008) incorporating the EM Iterative Closest Points (ICP) algorithm of Granger and Pennec (2002). Recently, work on using particle-based systems and minimum descriptor lengths have also been developed by Datar et al. (2011) and Davies et al. (2010) respectively. Xia et al. (2011) introduced a simple semi-automatic approach to construct a set of surfaces with correspondence via a template mesh that is not dependent on the genus of the surface. We improve upon the work of Xia et al. (2011) in creating a fully automatic template surface correspondence optimisation approach to create the training surfaces for the shape models in this paper. The procedure in Algorithm 5 utilises a single template mesh that is fitted to all training labels, resulting in a set of surfaces that all inherently have point-to-point correspondence. It is worth noting that although this method does not guarantee optimal correspondence such as the method of Davies et al. (2010), it was found to be sufficiently accurate for the work presented in this paper. The method can be repeated with the (rescaled) mean surface of the resulting model of the first run as the template surface to reduce possible biases introduced in the model by selecting an initial template surface. Algorithm 5 is similar to the method of Hufnagel et al. (2008), but instead of computing surface-to-surface

Algorithm 5. Template surface optimisation 1: A template surface is produced from a selected case using marching cubes (Lorensen and Cline, 1987), smoothed (Taubin et al., 1996) and decimated (Garland and Heckbert, 1997) to achieve relatively even sampling while having approximately 4–7 K points. 2: for each case do 3: Register template surface to the label marching cube surface via the ICP (Besl and McKay, 1992) after centroid and centroid size alignment. 4: The marching cube surface can be decimated to reduce computation time of the ICP. 5: Place this surface into the labelled image and freely deform towards the image gradient. 6: end for 7: for optimisation iterations do 8: Construct an initial shape model using the standard PCA and surfaces from previous loop. 9: for each case do 10: Fit the surface using an ASM approach in order to correct the pose. Use a low precision of the model except for the final iteration. 11: Freely deform the resulting surface toward the image gradient along the vertex normals to get the final surfaces. 12: end for 13: end for 14: return Set of s surfaces having point-to-point correspondence. Build shape model from these surfaces.

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

registrations and iteratively refining the shape model, the non-rigid registration is computed as a surface-to-image registration problem. The latter has the advantage of being faster and simpler to implement if an ASM approach is to be used in subsequent segmentations. The authors found that the entire set of 56 surfaces were optimised in less than 30 min on a 3.2 GHz Intel Xeon desktop PC for each bone. The optimisation process has only two parameters, namely the capture-range c, which defines the profile length as x ¼ 8  c and the maximum allowed displacement as m ¼ c, and the number of optimisation iterations. The authors have found empirically that one optimisation iteration often suffices for most applications. Tuning the capture-range is only a concern given very thin structures and a value of three usually suffices.4 Appendix B. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.media.2014. 02.002. References Altman, R.D., Bloch, D.A., Dougados, M., Hochberg, M., Lohmander, S., Pavelka, K., Spector, T., Vignon, E., 2004. Measurement of structural progression in osteoarthritis of the hip: the Barcelona consensus group. Osteoarthr. Cartilage 12 (7), 515–524. Audenaert, E.A., Baelde, N., Huysse, W., Vigneron, L., Pattyn, C., 2011. Development of a three-dimensional detection method of cam deformities in femoroacetabular impingement. Skeletal Radiol. 40 (7), 921–927. Baka, N., Kaptein, B., de Bruijne, M., van Walsum, T., Giphart, J., Niessen, W., Lelieveldt, B., 2011. 2D-3D shape reconstruction of the distal femur from stereo X-ray imaging using statistical shape models. Med. Image Anal. 15 (6), 840–850. Barratt, D.C., Chan, C.S., Edwards, P.J., Penney, G.P., Slomczykowski, M., Carter, T.J., Hawkes, D.J., 2008. Instantiation and registration of statistical shape models of the femur and pelvis using 3D ultrasound imaging. Med. Image Anal. 12 (3), 358–374. Besl, P.J., McKay, H.D., 1992. A method for registration of 3-D shapes. IEEE Trans. Pattern Anal. Mach. Intell. 14 (2), 239–256. Bourgeat, P., Fripp, J., Stanwell, P., Ramadan, S., Ourselin, S., 2007. MR image segmentation of the knee bone using phase information. Med. Image Anal. 11 (4), 325–335. Brechbühler, C., Gerig, G., Kübler, O., 1995. Parametrization of closed surfaces for 3D shape description. Comput. Vision Image Understan. 61 (2), 154–170. Burstein, D., Bashir, A., Gray, M.L., 2000. MRI techniques in early stages of cartilage disease. Invest. Radiol. 35 (10), 622–638. Chandra, S., Dowling, J., Shen, K.-K., Raniga, P., Pluim, J., Greer, P., Salvado, O., Fripp, J., 2012. Patient specific prostate segmentation in 3-D magnetic resonance images. IEEE Trans. Med. Imaging 31 (10), 1955–1964. Cootes, T.F., Taylor, C.J., Cooper, D.H., Graham, J., 1995. Active shape models-their training and application. Comput. Vision Image Understan. 61 (1), 38–59. Dalvi, R., Abugharbieh, R., Wilson, D., Wilson, D.R., 2007. Multi-contrast MR for enhanced bone imaging and segmentation. In: Annual International Conference of the IEEE Engineering in Medicine and Biology Society. 2007, pp. 5620–5623. Datar, M., Gur, Y., Paniagua, B., Styner, M., Whitaker, R., 2011. Geometric correspondence for ensembles of nonregular shapes. Med. Image Comput. Comput.–Assist. Intervention – MICCAI 14 (Pt 2), 368–375. Davies, R.H., Twining, C.J., Cootes, T.F., Taylor, C.J., 2010. Building 3-D statistical shape models by direct optimization. IEEE Trans. Med. Imaging 29 (4), 961–981. Davies, R.H., Twining, C.J., Taylor, C., 2008. Groupwise surface correspondence by optimization: representation and regularization. Med. Image Anal. 12 (6), 787– 796. Dice, L.R., 1945. Measures of the amount of ecologic association between species. Ecology 26 (3), 297–302. Ding, F., Leow, W.K., Howe, T.S., 2007. Automatic segmentation of femur bones in anterior-posterior pelvis X-ray images. In: Kropatsch, W.G., Kampel, M., Hanbury, A. (Eds.), Computer Analysis of Images and Patterns, Lecture Notes in Computer Science, vol. 4673. Springer, Berlin, Heidelberg, pp. 205–212. Dodin, P., Abram, F., Pelletier, J.-P., Martel-Pelletier, J., 2012. A fully automated system for quantification of knee bone marrow lesions using MRI and the osteoarthritis initiative cohort. J. Biomed. Graphics Comput. 3 (1), 51. Dong, X., Ballester, M.A.G., Zheng, G., 2007. Automatic extraction of femur contours from calibrated x-ray images using statistical information. J. Multimedia 2 (5). 4 The algorithm had an average DSC of 0.96 between the training surfaces and its corresponding CT training labels for both pelvis and femoral bones during the optimisation process.

577

Dowling, J.A., Lambert, J., Parker, J., Salvado, O., Fripp, J., Capp, A., Wratten, C., Denham, J.W., Greer, P.B., 2012. An atlas-based electron density mapping method for magnetic resonance imaging (MRI)–alone treatment planning and adaptive MRI–based prostate radiation therapy. Int. J. Radiat. Oncol. Biol. Phys. 83 (1), e5–e11. Duda, R.O., Hart, P.E., 1972. Use of the Hough transformation to detect lines and curves in pictures. Commun. ACM 15 (1), 11–15. Ecabert, O., Peters, J., Schramm, H., Lorenz, C., Von Berg, J., Walker, M., Vembar, M., Olszewski, M., Subramanyan, K., Lavi, G., Weese, J., 2008. Automatic modelbased segmentation of the heart in CT images. IEEE Trans. Med. Imaging 27 (9), 1189–1201. Fripp, J., Crozier, S., Warfield, S., Ourselin, S., 2010. Automatic segmentation and quantitative analysis of the articular cartilages from magnetic resonance images of the knee. IEEE Trans. Med. Imaging 29 (1), 55–64. Ganz, R., Parvizi, J., Beck, M., Leunig, M., Nötzli, H., Siebenrock, K.A., 2003. Femoroacetabular impingement: a cause for osteoarthritis of the hip. Clin. Orthop. Relat. Res. (417), 112–120. Garland, M., Heckbert, P.S., 1997. Surface simplification using quadric error metrics. In: Special Interest Group on Computer Graphics and Ineractive Techniques (SIGGRAPH). ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, pp. 209–216. Gower, J., 1975. Generalized procrustes analysis. Psychometrika 40 (1), 33–51. Heimann, T., Wolf, I., Meinzer, H.-P., 2002. Active shape models for a fully automated 3D segmentation of the liver – an evaluation on clinical data. In: Larsen, R., Nielsen, M., Sporring, J. (Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI, Lecture Notes in Computer Science, vol. 4191. Springer, Berlin, Heidelberg, pp. 41–48. Heimann, T., Meinzer, H., 2009. Statistical shape models for 3D medical image segmentation: a review. Med. Image Anal. 13 (4), 543–563. Heimann, T., Wolf, I., Meinzer, H.-P., 2006. Active shape models for a fully automated 3D segmentation of the liver – an evaluation on clinical data. In: Larsen, R., Nielsen, M., Sporring, J. (Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI, Lecture Notes in Computer Science, vol. 4191. Springer, Berlin, Heidelberg, pp. 41–48. Hodler, J., Trudell, D., Pathria, M.N., Resnick, D., 1992. Width of the articular cartilage of the hip: quantification by using fat-suppression spin-echo MR imaging in cadavers. Am. J. Roentgenol. 159 (2), 351–355. Hufnagel, H., Pennec, X., Ehrhardt, J., Ayache, N., Handels, H., 2008. Generation of a statistical shape model with probabilistic point correspondences and the expectation maximization–iterative closest point algorithm. Int. J. Comput. Assisted Radiol. Surg. 2 (5), 265–273. Ibanez, L., Schroeder, W., Ng, L., Cates, J., 2005. The ITK Software Guide: The Insight Segmentation and Registration Toolkit, second ed. Kitware, Inc.. Lamecker, H., Seebaß, M., Hege, H., Deuflhard, P., 2004. A 3D statistical shape model of the pelvic bone for segmentation. In: Proceedings of the SPIE, vol. 5370. pp. 1341–1351. Lorensen, W.E., Cline, H.E., 1987. Marching cubes: a high resolution 3D surface construction algorithm. Special Interest Group on Computer Graphics and Interactive Techniques (SIGGRAPH), vol. 21. ACM Press/Addison-Wesley Publishing Co., pp. 163–169. Martin, S., Troccaz, J., Daanen, V., 2010. Automated segmentation of the prostate in 3D MR images using a probabilistic atlas and a spatially constrained deformable model. Med. Phys. 37 (4), 1579–1590. Masjedi, M., Marquardt, C.S., Drummond, I.M.H., Harris, S.J., Cobb, J.P., 2012. Cam type femoro-acetabular impingement: quantifying the diagnosis using three dimensional head-neck ratios. Skeletal Radiol., 1–5. Mattes, D., Haynor, D., Vesselle, H., Lewellen, T., Eubank, W., 2003. PET-CT image registration in the chest using free-form deformations. IEEE Trans. Med. Imaging 22 (1), 120–128. Naish, J.H., Xanthopoulos, E., Hutchinson, C.E., Waterton, J.C., Taylor, C.J., 2006. MR measurement of articular cartilage thickness distribution in the hip. Osteoarthr. Cartilage 14 (10), 967–973. Neubert, A., Fripp, J., Engstrom, C., Schwarz, R., Lauer, L., Salvado, O., Crozier, S., 2012. Automated detection, 3D segmentation and analysis of high resolution spine MR images using statistical shape models. Phys. Med. Biol. 57 (24), 8357– 8376. Nguyen, N.T., Laurendeau, D., Branzan-Albu, A., 2007. A new segmentation method for MRI images of the shoulder joint. In: Proceedings of the Fourth Canadian Conference on Computer and Robot Vision. IEEE Computer Society, Washington, DC, USA, pp. 329–338. Nishii, T., Sugano, N., Sato, Y., Tanaka, H., Miki, H., Yoshikawa, H., 2004. Threedimensional distribution of acetabular cartilage thickness in patients with hip dysplasia: a fully automated computational analysis of MR imaging. Osteoarthr. Cartilage 12 (8), 650–657. Pearson, K., 1901. On lines and planes of closest fit to systems of points in space. Philos. Mag. 2 (6), 559–572. Pollard, T.C.B., Gwilym, S.E., Carr, A.J., 2008. The assessment of early osteoarthritis. J. Bone Joint Surg. British vol. 90-B (4), 411–421. Schmid, J., Kim, J., Magnenat-Thalmann, N., 2011. Robust statistical shape models for MRI bone segmentation in presence of small field of view. Med. Image Anal. 15 (1), 155–168. Seim, H., Kainmueller, D., Heller, M., Lamecker, H., Zachow, S., Hege, H.-C., 2008. Automatic segmentation of the pelvic bones from CT data based on a statistical shape model. In: Proceedings of the First Eurographics conference on Visual Computing for Biomedicine. Eurographics Association, Aire-la-Ville, Switzerland, pp. 93–100.

578

S.S. Chandra et al. / Medical Image Analysis 18 (2014) 567–578

Shen, K.-k., Fripp, J., Mériaudeau, F., Chételat, G., Salvado, O., Bourgeat, P., 2012. Detecting global and local hippocampal shape changes in Alzheimer’s disease using statistical shape models. NeuroImage 59 (3), 2155–2166. Shepherd, D., Seedhom, B., 1999. Thickness of human articular cartilage in joints of the lower limb. Ann. Rheum. Dis. 58 (1), 27–34. Skocˇaj, D., Leonardis, A., Bischof, H., 2007. Weighted and robust learning of subspace representations. Pattern Recognit. 40 (5), 1556–1569. Taubin, G., Zhang, T., Golub, G.H., 1996. Optimal surface smoothing as filter design. In: Proceedings of the 4th European Conference on Computer. ECCV. SpringerVerlag, pp. 283–292. Thomaz, C., Giraldi, G., da Costa, J., Gillies, D., 2010. A Simple and Efficient Supervised Method for Spatially Weighted PCA in Face Image Analysis. Technical report, Department of Computing, Imperial College London. Thomaz, C.E., Giraldi, G.A., 2010. A new ranking method for principal components analysis and its application to face image analysis. Image Vision Comput. 28 (6), 902–913. Tustison, N.J., Avants, B.B., Cook, P.A., Zheng, Y., Egan, A., Yushkevich, P.A., Gee, J.C., 2010. N4ITK: improved N3 bias correction. IEEE Trans. Med. Imaging 29 (6), 1310–1320. Vercauteren, T., Pennec, X., Perchant, A., Ayache, N., 2009. Diffeomorphic demons: efficient non-parametric image registration. NeuroImage 45 (1, Suppl. 1), S61– S72.

Xia, Y., Chandra, S.S., Salvado, O., Fripp, J., Schwarz, R., Lauer, L., Engstrom, C., Crozier, S., 2011. Automated MR hip bone segmentation. In: Digital Image Computing: Techniques and Applications (DICTA). pp. 25–30. Yokota, F., Okada, T., Takao, M., Sugano, N., Tada, Y., Sato, Y., 2009. Automated segmentation of the femur and pelvis from 3D CT data of diseased hip using hierarchical statistical shape model of joint structure. In: Yang, G., Hawkes, D., Rueckert, D., Noble, A., Taylor, C. (Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI, Lecture Notes in Computer Science, vol. 576. Springer, Berlin/Heidelberg, pp. 811–818. Yushkevich, P.A., Piven, J., Hazlett, H.C., Smith, R.G., Ho, S., Gee, J.C., Gerig, G., 2006. User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. NeuroImage 31 (3), 1116–1128. Zheng, G., Schumann, S., Ballester, M.A.G., 2010. An integrated approach for reconstructing a surface model of the proximal femur from sparse input data and a multi-resolution point distribution model: an in vitro study. Int. J. Comput. Assisted Radiol. Surg. 5 (1), 99–107. Zheng, Y., Barbu, A., Georgescu, B., Scheuering, M., Comaniciu, D., 2008. Fourchamber heart modeling and automatic segmentation for 3-D cardiac CT volumes using marginal space learning and steerable features. IEEE Trans. Med. Imaging 27 (11), 1668–1681. Zoroofi, R.A., Sato, Y., Nishii, T., Sugano, N., Yoshikawa, H., Tamura, S., 2004. Automated segmentation of necrotic femoral head from 3D MR data. Comput. Med. Imaging Graphics 28 (5), 267–278.