ieee transactions on medical imaging 1 123

IEEE TRANSACTIONS ON MEDICAL IMAGING

123

1


2

Body Composition Assessment in Axial CT Images using FEM-based Automatic Segmentation of Skeletal Muscle Karteek Popuri, Dana Cobzas, Nina Esfandiari, Vickie Baracos, Martin Jägersand

Index Terms—Finite Element Method (FEM), Thoracic CT, Abdominal CT, Muscle segmentation

I. I NTRODUCTION Body composition, i.e., the proportion of fat and muscle tissues in the human body is related to the risk factors associated with a host of medical conditions such as growth failure in children, obesity, cachexia syndromes (in chronic disease of lung, liver, heart or kidney), malnutrition, lipodystropy, metabolic syndrome and frailty. In particular, body composition has important implications for cancer patients. It has been found that the presence of a relatively high body fat content makes the patients more vulnerable to the onset or recurrence of several types of cancers [1]. On the other hand, sarcopenia, a wasting syndrome which involves the loss of muscle tissues has been correlated with poor response to chemotherapy treatment and reduction in overall survival of the patients [2]. The muscle and fat tissues are target locations for the water- and fat-soluble drugs respectively used for cancer treatment. Consequently, the proportions of these tissues are believed to determine the chemotherapy toxicity and efficacy. Similarly, in patients suffering from diabetes, body composition has been linked with metabolic alteration such as insulin resistance [3]. Therefore, it is of considerable interest to understand the complex relationships between body composition and the various health ailments. Revolutionary advances in body composition research were brought about by the introduction of computed tomography K. Popuri is with the Department of Computing Science, University of Alberta, Canada. e-mail:[email protected]. D. Cobzas and M. Jägersand are with the Department of Computing Science, University of Alberta, Canada. Nina Esfandiari, Vickie Baracos are with the Department of Oncology, University of Alberta, Canada.

Original

Manual segmentation

Thresholded

Thoracic CT Abdominal CT

Abstract—The proportions of muscle and fat tissues in the human body, referred to as body composition is a vital measurement for cancer patients. Body composition has been recently linked to patient survival and the onset/recurrence of several types of cancers in numerous cancer research studies. This paper introduces a fully automatic framework for the segmentation of muscle and fat tissues from CT images to estimate body composition. We developed a novel finite element method (FEM) deformable model that incorporates a priori shape information via a statistical deformation model (SDM) within the templatebased segmentation framework. The proposed method was validated on 1000 abdominal and 530 thoracic CT images and we obtained very good segmentation results with Jaccard scores in excess of 90% for both the muscle and fat regions.

Fig. 1. Illustration of the challenges in abdominal and thoracic CT segmentation. Muscle tissue (red), fat tissue (blue) and thoracic cavity (green). It can be seen that segmentation solely based on thresholding the muscle tissue [−29 150] and fat tissue [−190 − 30] HU ranges results in a lot of errors due to the significant overlap of intensities between the muscle tissues and neighboring organs.

(CT) imaging technique, which has a very high precision and specificity for different tissues in the human body [4], [5]. In body composition studies using CT images, it is sufficient to acquire 2D cross-sectional images taken at specific skeletal landmarks instead of whole body 3D scans. This is because the proportion of muscle and fat tissues at these specific skeletal landmarks correlates well with the whole body muscle to fat ratio [6], [7]. Further, this localized image acquisition also prevents the patients from unnecessary radiation exposure. In current practice, abdominal and thoracic 2D CT images taken at the 3rd lumbar vertebra (L3) and the 4th thoracic vertebra (T4) respectively have been widely used for body composition analysis [8], [9]. These studies rely on the manual segmentation of muscle and fat tissue regions from CT images using pre-defined windows of Hounsfeld units (HU, units of radiation attenuation) for each tissue. However, the manual segmentation of large databases of CT images used in these studies is not practical and hence automatic segmentation methods are needed. The segmentation of the fat region from CT images using automatic methods [10], [11] is relatively straightforward due to the unique HU range of the fat tissue [−190 − 30]. But, the automatic segmentation of the muscle region is quite challenging as there exists a significant overlap between the HU ranges of the muscle tissue [−29 150] and surrounding organs (see Figure 1). Very few works exist on muscle segmentation from CT images [12], [13]. Moreover, these works only consider the segmentation of a specific muscle group as opposed to the


segmentation of the total muscle region captured in the image. Our method is the first work to address the segmentation of the whole muscle area contained in a CT image and thus fully automate the estimation of body composition given the abdominal and thoracic axial CT images corresponding to the L3 and T4 vertebrae locations respectively. In this work, we propose a template-based segmentation method that uses a non-uniform finite element method (FEM) mesh to segment the muscle region in abdominal L3 and thoracic T4 images. In abdominal CT images, the muscle region has a well defined shape (see Figure 5a) except for the variation between different patients and the variability due to operator errors in the manual identification of the L3 vertebrae. But, this is not valid for thoracic CT images where the shape of the muscle does not conform to a specific class of shapes (see Figure 5b). However, in thoracic CT images the thoracic cavity containing the lungs, ribs and other organs (see Figure 1) exhibits a consistent shape across the patient population, modulo the inter-patient variation and the variability in the manual identification of the T4 vertebrae. Hence, in the case of abdominal CT images we utilize the a priori shape knowledge about the muscle region to disambiguate the muscle tissue from the nearby organs that have overlapping intensities. Whereas, in thoracic CT images we use a shape prior model of the thoracic cavity for discriminating the muscle tissue from the organs inside the thoracic cavity.

3

(ILR) space embedding of the muscle shape was employed for the segmentation of extensor and flexor muscles from MRI images via a convex energy minimization framework. Apart from statistical shape models, other approaches such as multi-atlas fusion and fuzzy C-means clustering have also been explored for the segmentation of psoas major [21] and paraspinal [22] muscle groups respectively from CT images. In this work, we take a template-based segmentation approach where a binary template defining an initial shape is deformed via non-rigid or deformable registration to match the region of interest in the input image. The desired segmentation boundary is then implicitly defined by the contour of the initial shape and the deformation field estimated between the template and the input image. Existing works on templatebased segmentation either use a non-parametric representation of the deformation field [23] or parametrize the deformation field using B-spline basis functions [24], [25] on a uniform discretization of the image domain. Shape models can be directly built from a training set of deformation fields and they are commonly referred to as statistical deformation models (SDM) [26]. Previously, SDMs have been explored for introducing shape prior constraints into the template-based segmentation process [27], [28]. In fact, in our earlier work we had used a B-spline based SDM with great success for the segmentation of muscle tissue from abdominal CT images [29]. B. Our contributions

A. Related work We briefly review the relevant literature on the segmentation of muscle and/or fat tissues from medical images. Further, as we adopt a template-based approach for image segmentation in this paper, a short survey of the current works on templatebased segmentation is also presented. As mentioned before, owing to the distinct HU range of the fat tissues, CT images are a natural choice for the estimation of fat content in the human body [10], [11]. However, there have been efforts to delineate fat tissue regions from magnetic resonance images (MRI) as well, using user-defined thresholds [14], [15] or automatic image-adaptive thresholding techniques [16], [17]. We note that these works additionally sub-segment the fat region into subcutaneous fat tissue (SFT) and visceral fat tissue (VFT) respectively. The general focus of the existing works on muscle segmentation has been to extract a specific group of muscles from whole body CT or MRI scans using statistical shape models. Kamiya et al. devised a rule-based expert system for the segmentation of the psoas major [12] and recuts abdominis [13] muscles from CT images, where the muscle shape was approximated by a simple quadratic function. A more elaborate representation of the muscle shape can be considered through the use of a point distribution model (PDM) constructed from a set of training shapes with manually annotated landmark correspondences. Such PDMs have been employed for providing shape information during the Markov random field (MRF) based segmentation of the calf muscle [18] and the snake based segmentation of the quadratus lumborum muscle [19] from MRI images. Alternatively, in [20] an isometric log-ratio

•

•

•

The main contribution in this paper is the introduction of a completely automatic segmentation framework for muscle and fat tissues in abdominal and thoracic CT images. Such a framework has not been reported before to the best of our knowledge. Further, our proposed segmentation framework handles both abdominal and thoracic CT images in a unified manner via the shape modeling approach. This paper is an extended version of our previous works [29], [30]. We developed a computationally efficient finite element method (FEM) based registration framework to solve the template-based segmentation problem. Earlier works on template-based segmentation parametrized the deformation field on a uniform mesh using either the nonparametric [23] or B-spline [24], [25] models. This is inefficient because the deformation field is computed with the same accuracy everywhere, even though detailed deformations are only needed along the contour of the initial shape in the template. Our proposed FEM-based approach remedies this issue by employing a non-uniform mesh well adapted to the contour of the initial shape in the template and through the use of Lagrange basis functions instead of the B-spline basis functions to parametrize the deformation field. Another contribution of our work is that we incorporate a SDM using the Lagrange basis parametrization into the FEM-based segmentation framework to better capture the shape deformations of the muscle and thoracic cavity regions respectively. As the SDM based on Lagrange basis parametrization is constructed from the deformations of a


4

non-uniform mesh, it encodes the shape variations more compactly than the B-spline based SDMs [27], [28], [29] improving the overall efficiency of our proposed FEMbased framework. • Recent publications using CT based body composition analysis include several hundred images, often done in duplicate, requiring many weeks of work by manual segmentation [31], [32], [33]. In this context, an important contribution of this work is the validation of the proposed segmentation framework on a large database of abdominal and thoracic CT images. This demonstrates that the proposed automated segmentation method can be employed to greatly reduce the time taken to do body composition analyses on large databases. The rest of the paper is organized as follows. In section II we describe the proposed template-based segmentation methodology. In section III we discuss the automatic segmentation framework for abdominal and thoracic CT images. In section IV we provide implementation details. In section V we present an evaluation our automatic CT segmentation system. Finally, in section VI we present the conclusion. II. S EGMENTATION VIA FEM- BASED DEFORMABLE REGISTRATION WITH A G AUSSIAN SDM The theoretical framework for our FEM-based templatebased segmentation method is presented in this section and an overview of the proposed system is shown Figure 2. Given a binary template (see Figure 2c), the method computes optimal segmentation of a binary image (see Figure 2b). This segmentation (see Figure 2d) is achieved by optimally deforming the template such that it matches the input image. Image deformations are defined using a FEM-based deformable registration framework that is described in Section II-A. This framework is adapted to template-based segmentation as described in Section II-B. For improving the convergence and to add robustness of the underlying minimization problem we constrain the FEM mesh deformations using a Gaussian SDM described in Section II-D. A. FEM-based deformable registration framework Given input I : Ω → R and template IT : ΩT → R images, where Ω, ΩT ⊂ Rd , the task of deformable or non-rigid registration is to find a dense deformation field U : ΩT → Rd such that the input image warped using the deformation field, I(x + U(x)) is similar to the template image IT . In a FEMbased framework, the deformation field U is approximated as a linear combination of a set of basis functions {φn }N n=1 : U(x; Θ) =

N X

Un φn (x)

∀x ∈ ΩT ,

(1)

n=1 Nd where Θ = [Un ]N . The basis functions {φn }N n=1 ∈ R n=1 are defined on a uniform or non-uniform tessellation of the template image domain ΩT given by the mesh M = N ({Pn }N n=1 , ∆) , where {Pn }n=1 denotes the nodes of the mesh and ∆ is the set of elements (triangles or rectangles - see Figures 3c, 3d). The deformable registration task is transformed

into finding the unknown nodal deformation field parameters Θ through the finite-dimensional multivariate minimization of an energy: Θ∗ = argmin ED (Θ; I, IT ) + αER (Θ),

(2)

Θ∈RN d

where ED is the data term which measures the similarity between the warped input and the template images, ER is the regularization term that enforces the smoothing constraints on the estimated deformation field and α is the regularization constant. Choosing the data term as the sum of squared differences (SSD) similarity measure and the regularization d R P T term as the diffusion regularizer 12 ΩT ∇Ui ∇Ui dx [34], i=1

we get the finite-dimensional formulations of these terms using the FEM approximation in Eq. 1 as: Z N X SSD (I(x + Un φn ) − IT (x))2 dx, (3) ED (Θ) = ΩT

diff ER (Θ) =

d X

n=1

N X

Uni Umi

i=1 m,n=1

Z

(∇φn · ∇φm ) dx,

(4)

ΩT

where Un = [Uni ]di=1 . In order to minimize Eq. 2, we use a semi-implicit fixedpoint iteration scheme that involves solving the following sparse system of N d linear equations at every step k (where IdN d is a N d × N d identity matrix): diff SSD (Θ(k) )+Θ(k) . )Θ(k+1) = −τ ∇Θ ED (IdN d +τ α∇Θ ER

(5) B. Template-based segmentation via FEM-based registration The above described FEM-based deformable registration methodology can be directly adapted for the template-based segmentation of a region of interest (ROI) from a given input image I. This requires choosing a binary image IT : ΩT → {0, 1} defining an initial shape of the ROI as the template. Then, the initial ROI in the template is deformed using FEM-based deformable registration to match the ROI to be segmented from the input image. In this paper, we address the muscle segmentation task by considering a templatebased segmentation method defined on binary input images I : Ω → {0, 1} (see Figure 2). In our method, the SSD data term in Eq. 3 is modified as: Z N X SSD MT · [I(x + Un φn ) − IT (x)]2 dx, (6) ED (Θ) = ΩT

n=1

where MT : ΩT → R is a mask defined on the template domain as: ( Hǫ (ΦT (x)) for Abdominal CT MT (x) = . (7) 1 for Thoracic CT

In the above, ΦT : ΩT → R is the signed distance transform of the binary template image IT and Hǫ (z) = 2 1 −1 z 2 (1 + π tan ǫ ) is the regularized Heaviside function. The mask MT is defined such that the data term is computed only within the mean muscle region shape in the abdominal CT


5

Template-based segmentation via FEM-based registration with a statistical deformation model

Binarization

Optimal parameters Θ∗

Gaussian deformation prior ESPCA (Θ, ΣΘ )

(a)

(b)

(c)

(d)

Fig. 2. Illustration of ROI segmentation in the abdominal and thoracic CT image cases. (a) Input image (b) Binarized input image I (c) Mean ROI (IT ) with a non-uniform discretization. Note the relatively higher density of nodes near the boundaries compared to elsewhere (d) Final estimated ROI I ∗ (green contour).

(a) Input image

(b) Template

(c) Uniform mesh (d) Non-uniform mesh

efficiency, as they can be naturally defined on non-uniform meshes that are adapted to the contours of the initial ROI shape in the template. We illustrate this fact using a simple example in Figure 3, where a “star” shaped ROI is segmented from a ∼ 200 × 200 input image. It can be seen that the ROI was successfully segmented using a Lagrange basis on a non-uniform mesh generated from just 36 nodes (see Figure 3d, Figure 3g), whereas the use of a B-spline basis on a uniform mesh with the same number of nodes (see Figure 3c, Figure 3f) results in a failed segmentation. A more quantitative comparison of the computational efficiency resulting from the use of Lagrange and B-spline basis functions is presented in Section V-A using real CT images.

(e) Initial tem- (f) Final contour (g) Final contour plate contour (B-spline) (Lagrange) Fig. 3. Comparison of FEM-based segmentation performance using a uniform (Bspline basis) VS non-uniform (Lagrange basis) mesh with the same number (N = 36) of nodes.

images, whereas in the thoracic CT images it is computed on the whole template domain. The template IT is then deformed through the energy minimization in Eq. 2 and the final segmentation label I ∗ : Ω → {0, 1} is given by inverse warping the template using the optimal nodal deformation parameters Θ∗ = [U∗n ]N n=1 as: I ∗ (x) = IT (x + U−1 (x; Θ∗ ))

∀x ∈ Ω.

(8)

C. FEM nodal basis functions We considered two types of basis functions for the FEMbased parameterization of the deformation field in Eq. 1, namely the piecewise-linear Lagrange basis functions (see Appendix A) [35] and the popular cubic B-spline basis functions [25]. However, we propose the use of Lagrange basis functions instead of the B-spline basis functions. This is because the Lagrange basis functions result in higher computational

D. Gaussian statistical deformation model Let us assume that we obtained a set of M nodal deforma(m) M tion field parameters {Θ(m) = [Un ]N n=1 }m=1 by registering (m) M M training images {I : Ω → R}m=1 to a template image IT , using the above discussed FEM-based deformable registration method. For introducing a priori shape knowledge into the future registration tasks, we now construct a statistical deformation model (SDM) from the nodal deformation field parameters {Θ(m) }M m=1 . Following [36], the space of deformation parameters is modeled using a multivariate Gaussian density N (Θ, ΣΘ ), with a sample mean Θ and a N d × N d sample covariance matrix ΣΘ . Further, the dominant modes of shape variation are computed using principal component analysis (PCA) and they are used to devise an additional shapebased regularization term as follows: 1 ||Θ − Θ||2 , γσ02 B = diag(η1 . . . ηK )[B1 . . . BK ]T , ESPCA (Θ) = ||B (Θ − Θ)||2 +

ΣΘ = (1 − γ)ΣΘ + γσ02 IdN d , ηk2 = ((1 − γ)σk2 + γσ02 )−1 − (γσ02 )−1 ,

(9)


where, {Bk }K k=1 is the PCA basis corresponding to the augmented covariance matrix ΣΘ . Here, σk2 are the eigen values of the matrix ΣΘ and γ, σ0 are constants. The above regularizer can be seen as imposing a shape prior on the deformation fields by penalizing deviations from the Gaussian SDM as opposed to strictly restricting the deformation fields to the span of the PCA basis. Incorporating the shape-based regularizer into Eq. 2 we obtain the following statistically constrained FEM-based deformable registration formulation: SSD diff Θ∗ = argmin ED (Θ) + αER (Θ) + βESPCA (Θ), (10) Θ∈RN d

where α, β are regularization constants. For minimizing the above energy, a semi-implicit scheme similar to equation Eq. 5 is devised by incorporating the additional shape-based gradient term ∇Θ ESPCA (Θ(k) ) as follows: diff SSD (Θ(k) ) (IdN d + τ α∇Θ ER )Θ(k+1) = −τ ∇Θ ED

− τ β∇Θ ESPCA (Θ(k) ) + Θ(k) . (11) The analytical expressions for the above gradient terms, SSD diff ∇Θ ED , ∇Θ ER , ∇Θ ESPCA are given in Appendix B. III. U NIFIED FRAMEWORK FOR AUTOMATIC SEGMENTATION OF MUSCLE AND FAT TISSUES We propose a unified framework for the segmentation of abdominal and thoracic images. The main idea of the proposed segmentation framework is to first determine a region of interest (ROI) in the given CT image, using the templatebased segmentation method presented in the Section II. The segmented ROI is then used to mask the CT image and the final muscle and fat tissue regions are determined using the pre-defined muscle and fat tissue HU ranges. The various steps of the proposed unified framework for abdominal and thoracic CT image segmentation are summarized below and illustrated in Figure 2: (0) ROI definition and binarization of input image: • Abdominal CT - The ROI corresponds to the muscle region. For segmenting the ROI, a binary version of the input CT image (see Figure 2b) is considered, which is obtained by setting the points in the image that lie within the HU range of the muscle tissue to 1. • Thoracic CT - Here, the ROI corresponds to the thoracic cavity and the binarized CT image (see Figure 2b) used for segmentation is obtained by setting the points in the input image that lie outside the HU ranges of the muscle and fat tissues to 1. (1) Computation of mean ROI shape: Given a training set of binary ROI shapes {I (m) : Ω → {0, 1}}M m=1 , an unbiased mean shape estimation approach [37] is followed to compute the mean ROI shape IT (see Figure 2c). The non-rigid registrations needed during the mean ROI shape computation are implemented using the FEM-based deformable registration method (see Section II-A). (2) Building a SDM of ROI shape: The possible variations of the ROI shape are represented by the set of nodal deformation field parameters {Θ(m) }M m=1 that are estimated through the FEM-based deformable registration of

6

the training shapes {I (m) }M m=1 with the mean shape IT . A Gaussian SDM is built to compactly encode these shape variations and the corresponding PCA-based regularizer ESPCA given in Eq. 9 is formulated for enforcing shape prior constraints in the next step. (3) ROI segmentation via SDM constrained FEM-based registration: In order to initialize the segmentation, the mean shape IT is affinely-aligned with the binarized image I. The mean shape IT is then deformed towards the binarized image I according to the minimization of the energy in Eq. 10 incorporating the shape-based regularizer ESPCA . The desired final ROI shape I ∗ (see Figure 2d) is determined by warping back the mean ROI shape IT using the inverse of the optimal deformation field (see Eq. 8). (4) Muscle and fat region segmentation: • Abdominal CT - The points in the input image lying within estimated ROI (muscle region) are selected to obtain the final muscle region and input image points not belonging to the final muscle region are selected to obtain the final fat region. • Thoracic CT - Both the final fat and muscle regions are determined by selecting the points in the input image that do not belong to the estimated ROI (thoracic cavity) using the known fat and muscle HU ranges. We note that, only the steps 0, 3 and 4 need to be performed for every given new input image, while the steps 1 and 2 are performed offline only once. IV. I MPLEMENTATION DETAILS The proposed FEM-based SDM constrained unified CT segmentation method was implemented using both the Bspline and the Lagrange basis functions in a multi-resolution framework, and was coded in MATLAB using the MEX facility. A. Non-uniform mesh generation The B-spline basis was defined on a uniform discretization of the template image domain. On the other hand, to take advantage of the Lagrange basis functions, a non-uniform mesh was generated on the template using an image-adaptive meshing strategy proposed by Yang et al. [38]. The central aspect of this meshing strategy is to place a relatively higher density of nodes in image regions with salient features as opposed to homogeneous regions with a uniform intensity profile. The salient features in the image are extracted based on a double derivative map of the image. The double derivative map is then converted into a binary image via halftoning. The “white” locations in the halftoned image correspond to the desired mesh nodes, which are input to Delaunay triangulation and refinement to generate the non-uniform mesh. The nonuniform meshes generated using this strategy on the binary templates containing the mean ROIs in the abdominal and thoracic cases are shown in Figure 2c. It can be clearly seen that a dense discretization is achieved only around the boundaries of the mean ROI, while the uniformly white/black regions in the template are discretized sparsely. Such an adaptive non-uniform mesh leads to a computationally efficient


parametrization of the deformation field using Lagrange basis functions. B. Multi-resolution strategy We followed the standard multi-resolution approach to perform template-based segmentation in a coarse-to-fine manner using 4 resolution levels. The input and template images at each successive resolution levels are obtained by successively downsampling the original images by a factor of two. For the B-spline basis, the uniform mesh at the next resolution was obtained by halving the mesh spacing at the current resolution. While for the Lagrange basis, a new non-uniform mesh was generated on the downsampled template image at each new resolution level. In both B-spline and Lagrange basis cases, a separate Gaussian SDM shape prior (see section II-D) was learned at each new resolution level. Starting at the coarsest resolution level with a zero initial deformation field, the estimated deformation field at the current level is multiplied by a factor of two and then upsampled to the next resolution level. It is then used to initialize the estimation of the nodal deformation parameters at the next resolution level and so on until the finest resolution level. V. E XPERIMENTS We evaluated the proposed FEM-based SDM constrained unified CT segmentation framework on a large database of CT images obtained from the Cross Cancer Institute (CCI), University of Alberta, Canada. The CCI dataset consisted of 1069 abdominal and 590 thoracic axial 512 × 512 CT images taken at the L3 and T4 level respectively from patients with head and neck cancers. The abdominal images were obtained from 670 patients and the thoracic images were obtained from 334 patients. Manual segmentations of the muscle and fat regions were available for all the images in the CCI dataset. The manual segmentation was performed by a single expert operator using Slice-O-Matic V4.3 software (Tomovision, Montreal, Canada). In the abdominal case, 69 images were used for training the SDM (see section II-D) and 1000 images were used for testing the proposed automatic segmentation framework. In the thoracic case, 60 and 530 images were used for training and testing respectively. In this section, we first compare the computational performance obtained by the B-spline and the Lagrange basis on a small subset of the CCI database. Then, we present the overall muscle and fat region segmentation results obtained on the complete CCI database. Both of the experiments below were run on a Intel i7 3.60 GHz machine with 64GB RAM. A. Comparison of FEM nodal basis functions The goal of this experiment is to determine the effect of the basis function choice on the computational performance of the proposed FEM-based SDM constrained segmentation method. We randomly selected 50 images from the test subset of the CCI database in each of the abdominal and thoracic cases. These images were segmented using both the Lagrange and the B-spline basis based implementations of the proposed

7

Evaluation Abdominal Thoracic Measure B-spline Lagrange B-spline Lagrange DOF 43520 8160 43520 3700 Time (sec) 14.79 ± 1.23 0.60 ± 0.09 16.06 ± 1.40 0.67 ± 0.44 Muscle Jaccard 90.16 ± 7.42 90.44 ± 7.49 91.34 ± 6.23 91.34 ± 6.19 (%) Fat Jaccard 91.20 ± 4.50 91.20 ± 4.50 89.74 ± 12.85 89.82 ± 12.57 (%) TABLE I C OMPARISON OF RESULTS OBTAINED USING B- SPLINE AND L AGRANGE BASIS . (A LL VALUES ARE REPORTED AS MEAN ± SD).

method. In Table I, we compare the B-spline and Lagrange basis performance using three different measures, the number of degrees of freedom (DOF) of the mesh, the computational time and the Jaccard score. The number of degrees of freedom (DOF) of a uniform or a non-uniform mesh is computed as twice the sum of the number of nodes used in the mesh at each multi-resolution level. The computational time corresponds to the time taken for solving the equation system in Eq. 11, which is the core optimization step of the proposed method. The Jaccard score quantifies segmentation accuracy by measuring the amount of overlap between the automatic and the manual labels. We performed pairwise t-test comparisons to determine the statistical significance of the differences in the performance of the Lagrange and B-spline basis respectively. It can be seen that, the use of the Lagrange basis results in a considerably lower number of DOFs (5 − 10 times lower) compared to the B-spline basis. This leads to the Lagrange basis being 10 − 20 times faster than the B-spline basis (significantly lower computational time, p < 0.01 ), while obtaining similar segmentation accuracies (no statistically significant differences found, p >> 0.01) for both muscle and fat regions on an average. Hence, as mentioned earlier in Section II-C, the implementation of the proposed method using the Lagrange basis is more computationally efficient. B. Validation on the CCI database The entire test subset of the CCI database consisting of 1000 lumbar and 530 thoracic images was segmented using the Lagrange basis based implementation of the proposed FEM-based SDM constrained segmentation method. Figure 4 shows the distribution (histogram) of segmentation accuracies obtained for the muscle and fat regions on the CCI test subset using our proposed method. For comparison, we also report the distribution of accuracies obtained with the direct thresholding-based segmentation of the CT images using the known muscle tissue [−29 150] and fat tissue [−190 −30] HU ranges. It can be seen that the proposed method on an average obtains high (> 90%) muscle and fat region segmentation accuracies in both the abdominal and thoracic CT image cases. Also, it can be noted that the proposed method performs poorly (< 80% Jaccard score) on only a small percentage (5 − 10%) of images in all the four different cases shown. As expected, the proposed method is able to achieve a vast improvement in the muscle region segmentation results compared to the

70 60 50 40 30 20 10 0

8

Proposed: 90.61 ± 7.42 Thresholding: 42.34 ± 8.87

0

10

20

30

40

50

60

70

No. of images (%)

No. of images (%)


80

90

80 70 60 50 40 30 20 10 0

100


0


0

10

20

30

40

50

60

70

No. of images (%)

No. of images (%)

80 70 60 50 40 30 20 10 0

80

10

20

30

40

50

60

70

80

90

100

90

100

(b) Thoracic muscle Jaccard scores (%)

(a) Abdominal muscle Jaccard scores (%)

90

100

(c) Abdominal fat Jaccard scores (%)

80 70 60 50 40 30 20 10 0


0

10

20

30

40

50

60

70

80

(d) Thoracic fat Jaccard scores (%)

Fig. 4. Thresholding-based vs FEM-based segmentation accuracies. Results from abdominal (N = 1000) and thoracic (N = 530) cases are shown. In (a)-(d), the distribution (histogram) of Jaccard scores is shown along with the mean ± SD for each of the thresholding-based and FEM-based methods.

baseline thresholding-based approach. Further, in the thoracic CT case, our proposed method obtains better segmentation accuracies even for the fat region. In Table II, we compare the tissue cross-sectional area estimates (cm2 ) obtained from the proposed automatic segmentation method and the manual segmentation method respectively. The tissue areas were computed by summing up tissue pixels and multiplying by the pixel surface area. The coefficient of variation (COV) between the proposed and the manual methods is between 2−6% which is comparable to the inter-operator COVs reported for manual segmentation [39], [40].

fail. In summary, the proposed segmentation method has demonstrated a very good segmentation performance on the CCI data set that consisted of a diverse collection of CT images. The images with major segmentation errors typically belonged to individuals who would be described clinically as cachexic (i.e., emaciated). Notably, these individuals lacked abdominal adipose tissue, subcutaneous adipose tissue and had very poor contrast between tissues as the wasting process leads to reduction of the attenuation values of muscle and increase of attenuation values of adjacent fat.

In Figure 5, we show sample visual results representative of the general segmentation performance of the proposed method in the abdominal and thoracic cases. In the abdominal case, the proposed method achieves very high (> 95%) muscle and fat region accuracies when the images do not have any abnormalities (Figure 5a, rows 1,2). The segmentation performance slightly degrades in the presence of pathologies such as kidney tumors (Figure 5a, row 3), that typically have overlapping intensities with the muscle tissue. A complete failure of the segmentation occurs only when either the muscle tissue has significantly atrophied resulting in a poor contrast between muscle and adjacent tissue (Figure 5a, row 4) or the individual is so emaciated that only a few flecks of fat maybe present (Figure 5a, row 5). In such situations, both manual and automated segmentation are extremely challenging. A similar trend is observed even in the thoracic case, where an excellent performance (> 95% Jaccard score) is observed in the absence of abnormalities (Figure 5b, rows 1,2), while the presence of lung tumors (Figure 5b, rows 3,4) results in relatively lower Jaccard scores for both the muscle and fat regions. In the thoracic case also, the presence of only a few flecks of sporadically found fat tissue (Figure 5b, row 5) causes the proposed segmentation method to completely

VI. D ISCUSSION We presented a fast and accurate method for the segmentation of muscle and fat tissues from CT images. We demonstrated the robustness of our method by validating it on a large database of CT images taken from cancer patients and obtained very good segmentation results. As the proposed method is completely automatic, it facilitates the undertaking of large scale cancer research studies that require the measurement of body composition.

L AGRANGE

A PPENDIX A NODAL BASIS FUNCTIONS

Consider a non-uniform mesh M = ({Pn }N n=1 , ∆) , where {Pn }N denotes the nodes of the mesh and ∆ is the set of n=1 elements. The Lagrange basis is defined as:   is linear within each adjacent elements,  φn (x) = 1 at each node Pn ,   0 at every other node Pm 6= Pn . (12)


Original CT

9

Muscle region

Fat region

96.34

96.79

96.46

95.10

90.10

Original CT

Muscle region

Fat region

95.56

99.06

98.09

97.45

90.44

93.60

75.59

71.10

88.01

0.00

91.99

31.08

79.01

74.96

22.12

(b) Thoracic

(a) Abdominal Fig. 5. CT muscle and fat region segmentation sample visual results. Manual label (green), automatic label (red) and overlap (yellow). The Jaccard scores (%) are mentioned above each of the respective images.

Abdominal Thoracic Proposed Proposed Manual COV Proposed Proposed Manual COV Jacc. (%) Area (cm2 ) Area (cm2 ) (%) Jacc. (%) Area (cm2 ) Area (cm2 ) (%) Muscle 90.61 ± 7.42 144.52 ± 38.01 143.68 ± 38.54 2.48 ± 3.36 92.01 ± 5.27 190.54 ± 46.82 182.66 ± 44.41 3.68 ± 4.24 Fat 90.72 ± 6.26 402.21 ± 227.62 386.87 ± 224.94 3.91 ± 3.49 90.96 ± 10.82 145.56 ± 80.59 138.91 ± 73.89 5.34 ± 11.86 Tissue

TABLE II M USCLE AND FAT TISSUE AREA ESTIMATES USING THE PROPOSED FEM- BASED AND MANUAL APPROACHES ON A BDOMINAL (N = 1000) AND T HORACIC (N = 530) IMAGES . T HE COEFFICIENT OF VARIATION (COV) BETWEEN THE PROPOSED AND MANUAL METHODS IS ALSO PRESENTED . A LL VALUES ARE REPORTED AS MEAN ± SD.

A PPENDIX B A NALYTICAL EXPRESSIONS FOR THE GRADIENT

TERMS

• The gradient of the data term in Eq. 6 is given as: SSD SSD (Θ)]N ∇Θ ED (Θ) = [∇Un ED n=1 , where : Z SSD [I(ξ(x; Θ)) − IT (x)] (Θ) = 2 · MT · ∇Un ED ΩT

∇I|x=ξ(x;Θ) φn (x) dx,

and ξ(x; Θ) = x +

N X

n=1

Un φn (x).

(13)

• The gradient of the diffusion regularizer in Eq. 4 is given as: diff diff (Θ)]N ∇Θ ER (Θ) = [∇Un ER n=1 , where : Z N X diff (Θ) = 2 Um ∇Un ER ∇φm · ∇φn dx. m=1

(14)

ΩT

• The gradient of the PCA-based shape regularization term in Eq. 9 is given as: ∇Θ ESPCA (Θ) = 2((γσ02 )−1 IdN d + BT B)(Θ − Θ). (15)


R EFERENCES [1] A. Kubo and D. Corley, “Body mass index and adenocarcinomas of the esophagus or gastric cardia: a systematic review and meta-analysis,” Cancer Epidemiology Biomarkers & Prevention, pp. 872–878, 2006. [2] W. Dewys, C. Begg, P. Lavin, P. Band, J. Bennett, J. Bertino et al., “Prognostic effect of weight loss prior tochemotherapy in cancer patients,” The American journal of medicine, pp. 491–497, 1980. [3] S. Otake, H. Takeda, Y. Suzuki, T. Fukui, S. Watanabe et al., “Association of visceral fat accumulation and plasma adiponectin with colorectal adenoma: evidence for participation of insulin resistance,” Clinical Cancer Research, pp. 3642–3646, 2005. [4] R. Baumgartner, K. Koehler, D. Gallagher, L. Romero et al., “Epidemiology of sarcopenia among the elderly in new mexico,” American journal of epidemiology, pp. 755–763, 1998. [5] R. Ross, D. Dagnone, P. Jones, H. Smith et al., “Reduction in obesity and related comorbid conditions after diet-induced weight loss or exerciseinduced weight loss in mena randomized, controlled trial,” Annals of Internal Medicine, pp. 92–103, 2000. [6] W. Shen, M. Punyanitya, Z. Wang et al., “Total body skeletal muscle and adipose tissue volumes: estimation from a single abdominal crosssectional image,” Journal of Applied Physiology, pp. 2333–2338, 2004. [7] M. Mourtzakis, C. Prado, J. Lieffers et al., “A practical and precise approach to quantification of body composition in cancer patients using computed tomography images acquired during routine care,” Applied Physiology, Nutrition, and Metabolism, pp. 997–1006, 2008. [8] C. Prado, V. Baracos, L. McCargar, M. Mourtzakis et al., “Body composition as an independent determinant of 5-fluorouracil–based chemotherapy toxicity,” Clin. Canc. Res., pp. 3264–3268, 2007. [9] C. Prado, J. Lieffers, L. McCargar, T. Reiman, M. Sawyer et al., “Prevalence and clinical implications of sarcopenic obesity in patients with solid tumours of the respiratory and gastrointestinal tracts: a population-based study,” Lancet Oncol., vol. 9, no. 7, pp. 629–635, 2008. [10] T. Yoshizumi, T. Nakamura, M. Yamane, A. Islam, M. Menju et al., “Abdominal fat: standardized technique for measurement at ct1,” Radiology, pp. 283–286, 1999. [11] B. Zhao, J. Colville, J. Kalaigian et al., “Automated quantification of body fat distribution on volumetric computed tomography,” Jour. of comp. assis. tomography, vol. 30, no. 5, p. 777, 2006. [12] N. Kamiya, X. Zhou, H. Chen, T. Hara, H. Hoshi, R. Yokoyama, M. Kanematsu, and H. Fujita, “Automated recognition of the psoas major muscles on x-ray ct images,” in IEEE EMBS, 2009, pp. 3557–3560. [13] N. Kamiya, X. Zhou, H. Chen, C. Muramatsu, T. Hara, R. Yokoyama, M. Kanematsu, H. Hoshi, and H. Fujita, “Automated segmentation of recuts abdominis muscle using shape model in x-ray ct images,” in IEEE EMBS, 2011, pp. 7993–7996. [14] S. Gronemeyer, R. Steen et al., “Fast adipose tissue (fat) assessment by mri,” Magnetic resonance imaging, pp. 815–818, 2000. [15] J. Mattei, Y. Fur, N. Cuge, S. Guis, P. Cozzone, and D. Bendahan, “Segmentation of fascias, fat and muscle from magnetic resonance images in humans: the dispimag software,” Magnetic Resonance Materials in Physics, Biology and Medicine, pp. 275–279, 2006. [16] J. Kullberg, H. Ahlström, L. Johansson, and H. Frimmel, “Automated and reproducible segmentation of visceral and subcutaneous adipose tissue from abdominal mri,” International Journal of Obesity, pp. 1806– 1817, 2007. [17] H.-P. Müller, F. Raudies, A. Unrath, H. Neumann, A. Ludolph, and J. Kassubek, “Quantification of human body fat tissue percentage by mri,” NMR in Biomedicine, vol. 24, no. 1, pp. 17–24, 2011. [18] C. Wang, O. Teboul, F. Michel, S. Essafi, and N. Paragios, “3d knowledge-based segmentation using pose-invariant higher-order graphs,” in MICCAI. Springer, 2010, pp. 189–196. [19] C. Engstrom, J. Fripp, V. Jurcak, D. Walker, O. Salvado, and S. Crozier, “Segmentation of the quadratus lumborum muscle using statistical shape modeling,” Journal of Magnetic Resonance Imaging, vol. 33, no. 6, pp. 1422–1429, 2011. [20] S. Andrews, G. Hamarneh, A. Yazdanpanah, B. HajGhanbari, and W. D. Reid, “Probabilistic multi-shape segmentation of knee extensor and flexor muscles,” in MICCAI. Springer, 2011, pp. 651–658. [21] S. Meesters, F. Yokota, T. Okada, M. Takaya, N. Tomiyama, J. Yao, M. Liguraru, R. Summers, and Y. Sato, “Multi atlas-based muscle segmentation in abdominal ct images with varying field of view,” 2012. [22] Y. Wei, X. Tao, B. Xu, and A. Castelein, “Paraspinal muscle segmentation in ct images using gsm-based fuzzy c-means clustering,” Journal of Computer and Communications, p. 70, 2014. [23] K. A. Saddi, C. Chefd’hotel, M. Rousson, and F. Cheriet, “Region-based segmentation via non-rigid template matching,” in MMBIA, 2007.

10

[24] X. Huang, D. Metaxas, and T. Chen, “Metamorphs: Deformable shape and texture models,” in Computer Vision and Pattern Recognition, 2004. [25] C. Li and Y. Sun, “Active image: A shape and topology preserving segmentation method using b-spline free form deformations,” in ICIP. IEEE, 2010, pp. 2221–2224. [26] D. Rueckert, A. Frangi, and J. Schnabel, “Automatic construction of 3D statistical deformation models using non-rigid registration,” in MICCAI, 2001, pp. 77–84. [27] X. Huang, Z. Li, and D. Metaxas, “Learning coupled prior shape and appearance models for segmentation,” MICCAI, pp. 60–69, 2004. [28] M. Taron, N. Paragios, and M. Jolly, “Registration with uncertainties and statistical modeling of shapes with variable metric kernels,” TPAMI, vol. 31, no. 1, pp. 99–113, 2009. [29] H. Chung, D. Cobzas, J. Lieffers, L. Birdsell, and V. Baracos, “Automated segmentation of muscle and adipose tissue on ct images for human body composition analysis,” in SPIE Medical Imaging, 2009. [30] K. Popuri, D. Cobzas, M. JaÌ´Lgersand, N. Esfandiari, and V. Baracos, “Fem-based automatic segmentation of muscle and fat tissues from thoracic ct images,” in ISBI, 2013. [31] L. Martin, L. Birdsell, N. MacDonald, T. Reiman, M. T. Clandinin, L. J. McCargar, R. Murphy, S. Ghosh, M. B. Sawyer, and V. E. Baracos, “Cancer cachexia in the age of obesity: skeletal muscle depletion is a powerful prognostic factor, independent of body mass index,” Journal of Clinical Oncology, pp. JCO–2012, 2013. [32] N. Fujiwara, H. Nakagawa, Y. Kudo, R. Tateishi, M. Taguri, T. Watadani, R. Nakagomi, M. Kondo, T. Nakatsuka, T. Minami et al., “Sarcopenia, intramuscular fat deposition, and visceral adiposity independently predict the outcomes of hepatocellular carcinoma,” Journal of hepatology, 2015. [33] Y. Hamaguchi, T. Kaido, S. Okumura, T. Ito, Y. Fujimoto, K. Ogawa, A. Mori, A. Hammad, E. Hatano, and S. Uemoto, “Preoperative intramuscular adipose tissue content is a novel prognostic predictor after hepatectomy for hepatocellular carcinoma,” Journal of hepato-biliarypancreatic sciences, vol. 22, no. 6, pp. 475–485, 2015. [34] J. Modersitzki, Numerical methods for image registration. Oxford Univ. Press, 2004. [35] O. Zienkiewicz and R. Taylor, The finite element method for solid and structural mechanics. Elsevier Butterworth-Heinemann, 2006. [36] T. Albrecht, M. Luthi, and T. Vetter, “A statistical deformation prior for non-rigid image and shape registration,” in CVPR, 2008, pp. 1–8. [37] A. Guimond, J. Meunier, and J. Thirion, “Average brain models: A convergence study,” CVIU, vol. 77, no. 2, pp. 192–210, 2000. [38] Y. Yang, N. Wernick, and G. Brankov, “A fast approach for accurate context-adaptive mesh generation,” IEEE transactions on image processing, vol. 12, no. 8, pp. 866–881, 2003. [39] E. Demerath, K. Ritter, W. Couch, N. Rogers, G. Moreno et al., “Validity of a new automated software program for visceral adipose tissue estimation,” Int. jour. of obesity, vol. 31, no. 2, pp. 285–291, 2006. [40] A. J. MacDonald, C. A. Greig, and V. Baracos, “The advantages and limitations of cross-sectional body composition analysis,” Current opinion in supportive and palliative care, vol. 5, no. 4, pp. 342–349, 2011.