Shape estimation of subcutaneous adipose tissue using an articulated statistical shape model

The quantification of adipose tissue can give insights into obesity and the associated diseases. The segmentation of this tissue is however difficult, particularly due to the complex geometry and the local appearance similarities between the subcutaneous and visceral types. Shape priors can be used to regularise the segmentation of geometries with small variation in shape. However, human bodies are articulated and substantially different across subjects. In this paper, a novel method is proposed for the segmentation of the subcutaneous fat layer in full-body magnetic resonance imaging scans. The proposed method is based on a statistical shape model of the whole body surface, which is learned from geometric scans. The body model is factorised into pose and shape deformations, which allows a compact parametrisation of large variations in human shape. The proposed method is applied in the segmentation of five magnetic resonance imaging data-sets. The experiments show that the proposed model can be used to effectively segment the subcutaneous fat geometry in subjects with different body mass indices. The incorporation of the statistical model in the algorithm regularises the segmentation, and establishes correspondences between the subcutaneous fat layer of the geometries across subjects. The registration of the fat layer with a common geometry could facilitate the statistical analysis of the shape distribution across the different geometries.


Introduction
The modeling and analysis of fat distribution in the human body is challenging, and has a wide range of applications. The distribution of adipose tissue in the human body is of increasing interest in many research studies with metabolic background (Thomas, Fitzpatrick, et al. 2013). Visceral adipose tissue (VAT) has shown to be metabolically highly active while subcutaneous adipose tissue (SCAT) was assumed to have little clinical meaning. However, as SCAT also contributes to obesity, its quantification as well as the determination of its distribution across the body is essential for phenotyping of subjects at increased risk for metabolic diseases. For these reasons, SCAT segmentation is important for supervising changes in body composition during lifestyle intervention and weight control. The modeling and simulation of the human body fat requires an accurate geometric representation of the fat layer. In particular, the geometry of the subcutaneous fat layer, which exhibits viscoelastic characteristics, should be included in the modeling and analysis of body dynamics.
The recent advances in medical imaging have allowed a reliable quantification of whole body adipose tissue. Amongst the different imaging modalities, computed tomography (CT) imaging and magnetic resonance imaging (MRI) have been widely used by clinicians for the quantification of adipose tissue. However, research in this area rarely differentiates between SCAT and VAT tissue, or does so only in the abdominal area.
The geometric extraction of the subcutaneous fat layer in the human body is a challenging task. Manual segmentation of the fat layer from biomedical image data can be extremely tedious due to the size of the image data and the complexity of the SCAT geometry. The similar intensities of the subcutaneous adipose tissue and visceral adipose tissue also makes their automatic classification very challenging. Factors such as image noise, fuzzy object edges and intensity inhomogeneity, inherent to imaging techniques such as MRI, also affect the accuracy of the segmentation. Although several techniques have been designed for the automatic quantification of the adipose tissue, they are often applied to single or regional slices of the image dataset.
In this paper, we present an algorithm that estimate the shape of the subcutaneous layer from whole-body MRI scans. The main novelty with respect to previous work is the usage of a statistical, articulated model of the body surface, created from thousands of 3D surface scans (Anguelov, Srinivasan, et al. 2005). Unlike previous methods which exploit shape priors or other appearance models (Chen and Radke 2009;Cremers, Osher and Soatto 2006;Kainmueller 2015;Rousson and Paragios 2002), we model the whole body surface with a compact parametrization that decouples body model deformations due to pose, and deformations due to body shape. This has three main advantages. Firstly, by having a strong model of the exterior surface, every data point fit by the model helps in the inference of the configuration of the full body geometry. This helps to deal with noise, large slice gaps, or fuzzy object edges in the scans. Secondly, by fitting a model with a common geometry to all MRI scans, we can establish correspondences between SCAT geometries across different subjects. This will facilitate the creation of a statistical model of the SCAT distribution across the population, which can be useful for the analysis of SCAT between the datasets. Subsequently, the usage of an articulated model allows us to deal with the inevitable differences in pose existing in different scans.

Related work
There has been a recent increased use of medical imaging as a diagnostic tool. The different imaging modalities such as MRI have been widely used for anatomical shape and functional analysis. Several techniques have been designed for the extraction of various geometries from biomedical datasets. It is however a great challenge for current state of the art techniques such as the statistical methods in (Michailovich, Rathi and Tannenbaum. 2007;Wang, Hu and Pheng 2012), alignment based techniques (Christensen, Joshi and Miller 1997;Maintz and Viergever 1998;Seixas, Damasceno, et al. 2007), graph based techniques (Boykov and Funka-Lea 2006) and deformable models (Malladi, Sethian and Vemuri 1995;Caselles, Kimmel and Sapiro 1997;Paragios and Deriche 2002;Cremers, Rousson and Deriche 2007;Yeo, Xie, et al. 2011) to segment the 3D geometry of the subcutaneous fat layer from MRI dataset.
Several techniques (Jin, Imielinska, et al. 2003;Roullier, Cavaro-Menard, et al. 2006;Leinhard, Johansson, et al. 2008;Positano, Christiansen, et al. 2009;Sussman, Yao and Summers 2010;Zhou, Murillo and Peng 2011;Mosbech, Pilgaard, et al. 2011;Thomas, Fitzpatrick, et al. 2013) have been designed for the automatic quantification of the adipose tissue of single or regional slices of the image dataset. Such segmentation techniques are typically based on the coherence of fat tissue intensity and its spatial continuity, using algorithms like fuzzy clustering, thresholding and various heuristics. However, there are only a few articles in the literature which applied segmentation algorithms on whole body adipose tissue. In (Brennan, Whelan, et al. 2005), simple techniques such as image thresholding and heuristics about voxel label connectivity are used to segment the adipose tissue from whole body MRI. The different adipose tissues, i.e. subcutaneous and visceral, are not categorized in this technique. The fuzzy clustering is used in (Wurslin, Machann, et al. 2010) for adipose tissue segmentation. The MRI dataset is partitioned into different anatomic sections, and an active contour model is used to delineate the subcutaneous and visceral adipose tissue at the abdominal section. An interesting SCAT segmentation approach which deals with the contact between parts in the MR datasets is proposed in (Fouquier, Anquez, et al. 2011). Different regions of the body are roughly segmented by analysing the area inside the body for each axial slice. However, the curve based segmentation can be sensitive to extreme variations of shape, in which some assumptions (i.e. smaller waist than hips) required by the technique can be violated.
In (Wald, Schwarz, et al. 2011), a technique was designed to quantify the amount of different adipose tissues from whole body MRI data. The fat tissue is initially extracted by automatic thresholding of the intensity. A statistical shape model of the abdomen is then used to quantify the subcutaneous and visceral adipose tissue in the abdominal region. In (Kullberg, Johansson, et al. 2009), a complicated technique based on the reconstructed fat and water MR signal components is used in the labelling of the different parts of the body. This approach requires a mixture of segmentation techniques tailored for the different regions of the body to segment the fat voxels. Different intensity thresholds, and morphological techniques such as dilation and erosion are used to segment some parts of the dataset. The voxels denoting the pelvis section are segmented by a fuzzy clustering technique and a statistical model of the pelvis shape. Intensity thresholding is used in (Muller, Raudies, et al. 2011) to identify voxels corresponding to fat tissue. The SCAT voxels are then delineated by iteratively extracting the voxels which are within an intensity range and adjacent to the labelled voxels. This technique however has difficulties dealing with thin structures, and regions with fuzzy edges and intensity inhomogeneity.
In general, our approach is different from the existing techniques as it uses an explicit model of the full body geometry. Instead of applying a variety of techniques to identify specific sections of the body, the proposed approach automatically aligns the full body surface to a generic template, in a coherent manner across body regions. The prior knowledge about body shape helps in overcoming the difficulties due to noise or lack of contrast, since errors in the segmentation would cause unlikely configurations according to the body model.
The model based techniques have been widely used for biomedical segmentation. Conventional model based segmentation techniques of medical images focused on modeling the variation of specific anatomical geometries delineated from the medical image dataset. In (Cootes, Taylor, et al. 1995), the training shapes are represented using control points, and principal component analysis (PCA) is used to model the shape variation. Various groups (Rousson and Paragios 2002;Cremers, Osher and Soatto 2006;Kim, Cetin and S. 2007) have incorporated shape prior into the level set technique. In general, several of the shape models are based on statistical assumptions that the training shapes are Gaussian distributed. Such segmentation techniques are limited to the extraction of simple geometries, and are not appropriate for the segmentation of the complex geometry of the fat layer. In (Cremers, Osher and Soatto 2006;Kim, Cetin and S. 2007;Chen and Radke 2009), the kernel density estimation (KDE) is applied to model the shape distribution, which allows the model to handle a larger variation of shapes. However, the complexity of the fat geometry makes it difficult for the non-parametric technique to model the shape variation.
The human geometry is articulated, and a small variation in pose can generate a substantial variation in shape. The high articulability of the body makes it difficult to model the exterior geometry, as well as the geometry of their subcutaneous fat layer with such shape priors. Also, the statistical modeling of the anatomical geometries require a large training dataset consisting of manually segmented images, which are not easily available for complex geometries such as that of the subcutaneous fat layer. In this paper, a statistical body shape model is learned from geometric scans, and the statistical model is factorized into pose and shape, which allows the modeling of a large variation of geometries. The recent technique in (Kainmueller 2015) is particularly relevant, in which deformable statistical models are applied to the segmentation of MRI scans. However, the segmentation algorithm was applied to single joints. The articulated model proposed in this paper is applied to the segmentation of the whole body surface from MRI data, which is more challenging.

Method
The presented algorithm aims to obtain a geometric segmentation of the SCAT layer (see Fig. 1) which is consistently parametrized across subjects. This can be achieved by fitting an articulated body model to the body surface of the subject, as well as to the interior surface of the SCAT. The body model is fit to the exterior (interior) surface by estimating the model parameters that maximize the (negated) image gradient evaluated along the normal of the body model surface.

Volumetric image data
The MRI data was captured with a 1.5 T scanner (Magnetom Sonata, Siemens Healthcare) after overnight fasting. The subjects were in prone position with extended arms. Axial T1-weighted fastspin-echo images were recorded from feet to hands applying the following measurement parameters: 10-mm-thick sections, 10-mm gap between sections, 256 × 192 matrix, 450-530mm field of view depending on the extension of the subject, and five sections per acquisition.

Statistical body model
The main novelty in this work is the usage of a statistical body model of the exterior body surface for regularizing the image alignment. The model used is based on (Anguelov, Srinivasan, et al. 2005), in which deformations are parametrized by pose (related to the configuration of the subject's articulations) and shape (related to the intrinsic geometry of a given subject). This model has been previously used in the alignment of textured surface scans (Bogo, Romero, et al. 2014), but to the best of our knowledge not in biomedical image data such as MRI. While the main challenge of dealing with surface scans lies on identifying correspondences between the model and the scan surface, the alignment of geometry to MRI data is more difficult as the scan surface (i.e. geometry represented by the edge voxels) itself has to be estimated from the intensity images.
The body model transforms a triangulated mesh representing an average person in a rest pose into a mesh representing a specific subject in an arbitrary posture. These meshes are represented by their vertex positions T and M in Eq. 1. However, the model describes the translation-free 3 × 3 deformation matrices applied to two edges per triangle. The linear operator A computes edge vectors as differences between adjacent vertices. The deformations to be applied to each template triangle f are decomposed into skeleton deformations B(θ) parametrized by the body part relative orientations θ, shape deformations D(β) parametrized by shape parameters β, and pose dependent deformation Q(θ), which models non-skeleton deformations like muscle bulging The components B, D and Q learned from thousands of surface scans as explained in (Anguelov, Srinivasan, et al. 2005), and are not further optimized in the segmentation process presented here.
The edges corresponding to the resulting triangles (AM ) f represent the transformed, translationfree differences between neighboring vertices. Vertices M are obtained by stitching the triangles in a least-squares sense (minimizing the difference between the left and right hand side of Eq. 1). Note that this stitching process implicitly solves for the translation of each vertex by imposing neighboring edges to be connected.

Exterior fat geometry
Let I : Ω → R, Ω ⊂ R 3 be an image function mapping the image domain Ω into scalar image intensities, which represents the MRI scan. In order to evaluate I in arbitrary points, we perform cubic spline interpolation, with zero-padding in the exterior. We define a closed 2D surface embedded in 3D space S : Φ → R 3 , Φ ⊂ R 2 , as the boundary to be fitted to the body, with normals N : Φ → R 3 . The surface S should have three characteristics. Firstly, the normals N should be aligned to strong negative image gradients −∇I (positive in the case of the interior surface , where intensities should increase along the normal). Secondly, S should resemble a human body, and therefore be close to the space of body shapes that can be generated with the body model. Thirdly, S should be smooth. The three characteristics are expressed as objective functions in the remainder of the section and are incorporated in the segmentation technique in a weighted manner where the continuous surface S has been replaced by its discretization into vertex positions v, done according to the triangulation of the model template T . The vertices v which minimize Eq. 2 are regarded as the resulting exterior surface registration. Optimization is performed with a newton conjugate-gradient method (Nash 1984) in a staged manner. In subsequent stages, the regularization weights α c and α s are reduced progressively. To obtain the gradients of E cpl and June 9, 2016 Computer Methods in Biomechanics and Biomedical Engineering cmbbe1 E smooth , we make extensive use of the auto-differentiation library Chumpy 1 with performancecritical parts implemented explicitly in C, in which the derivatives of complex functions can be computed using the chain rule. The objective E align is differentiated following (Kimmel 2003). The model parameters are initialized with an average body shape and an average prone pose, as can be seen in Fig. 1.

Image alignment
The alignment of the surface normals with image gradients is a problem that has been tackled previously, e.g. (Holtzman-Gazit, Kimmel, et al. 2006). Here we use a similar formulation, but we include a multiscale approach that allows our energy to have larger basins of attraction around the minima, making it suitable for gradient-based optimization. The energy takes the following form in the continuous domain and in the discrete version in which G σ is a Gaussian function with standard deviation σ, w σ is a weight applied to each scale of the Gaussian, x is a point inside the area differential da, A is the total area of the surface S, and v and n are, respectively, point samples on the surface S and normal samples on N , with corresponding weight w i related to the area associated to that point. In order to obtain point correspondences across subjects, we use the triangulation of the model template T as the coherent point samples, which are assumed to be roughly uniformly distributed across the surface. This leads to constant wi A = 1 {v} (where {v} denotes the cardinality of the set of vertices), which has no effect in terms of energy minimization. The weights w σ vary across stages in the optimization, starting with uniform weights and then increasing the weights associated to low σ. We estimate the derivative of E align inspired by the derivation in (Kimmel 2003)

Coupling external surface to the body model
Factors such as noise, sparse sampling and lack of contrast of the voxels can affect the alignment of the geometry. Therefore, regularization components are incorporated to constrain the deformation of the geometry. The first regularization term encourages v to resemble a human surface. This is achieved by minimizing the distance between the triangles containing vertices v and the model triangles resultant of the optimized shape β and pose θ parameters

Geometric smoothness
Subsequently, a smoothness regularization is incorporated in the segmentation model. While Eq. 5 requires the segmentation to be close to the model surface, it does not require it to be smooth. A simple way to encourage smoothness of the geometry is to include a small penalty for distances Figure 1. Segmentation of the SCAT layer from MRI scan: (from left to right) exterior geometry initialization, segmented exterior geometry, segmented interior geometry, orthogonal slices of the segmented geometries between neighboring vertices, particularly when their associated image gradient is small in which N e(i) is the set of neighboring vertices of vertex i according to the template triangulation.

Interior fat geometry
Once the vertices v that minimize the exterior geometry objective E ext have been computed (see second column in Fig. 1), the subcutaneous inner surface, represented by its vertices v can be aligned to the image geometry. The objective E int to be optimized is very similar to E ext , with two exceptions. As the normals of the interior surface should be aligned with positive image gradients (dark to bright voxels), the image intensities should be reversed in equation 4, leading to the alignment function of E align (I) = E align (−I). The objective E int for the alignment of the interior geometry can therefore be defined as The coupling to the model E cpl is replaced by a term penalizing differences between the interior triangles (Av ) f and their exterior counterparts (Av) f , defined as This term expresses that exterior and interior geometry should have similar shapes. Similar to the exterior geometry, a weighted sum of the energies is optimized The optimization is also performed in a staged manner, with a newton conjugate-gradient method.
In particular, the interior segmentation is initialized with the exterior geometry which is slightly displaced along the normal direction, i.e. v = v − n in the experiments, in which is a small constant. Since the head, hands and feet scans tend to be noisy and have little amount of fat, they are initialized as the vertices of the exterior geometry, i.e. = 0, and their effect in E align is suppressed.

Results and Discussion
In this section, we show the experimental results of the segmentation algorithm applied to five MRI scans from subjects with varying body mass index (BMI). The proposed method is evaluated qualitatively and quantitatively by comparing the extracted geometries with manually segmented slices of the MRI datasets. To our knowledge, there is no relevant benchmark in the literature with ground truth subcutaneous fat annotated in whole-body MRI scans in which segmented SCAT geometries could be evaluated. One of the objectives of designing the segmentation technique is to apply the segmentation to a large corpus of scans to facilitate the creation of such benchmark. With the current experimental results, we would like to show the potential of the proposed method to the computer vision, graphics and medical communities.  Figure 1 depicts the segmentation of the subcutaneous fat layer from the MRI dataset of BMI 23 using the proposed model. The parameterized model is initialized with a mean shape and prone June 9, 2016 Computer Methods in Biomechanics and Biomedical Engineering cmbbe1 pose from the training dataset. The shape and pose of the initial model is substantially different from those of the MRI geometry. As shown in Figure 1, the proposed model can easily deal with the arbitrary initialization, as the shape β and pose θ of the model are varied such that the vertices v are aligned to the edge voxels representing the exterior SCAT layer in the MRI dataset. This is a crucial factor for dealing with large sets of heterogeneous MRI scans. Figure 3. Segmentation of the subcutaneous fat layer from MRI scans (from left to right, increasing BMIs 19,22,23,25,35): the top row depicting coronal slices, and the bottom row depicting the manually labelled slices (yellow) and the corresponding axial slices segmented (orange) using the proposed model. Figure 2 depicts the corresponding axial slices of the segmented geometries of the subcutaneous fat layer of MRI dataset BMI 23. The delineation of the shapes shown in this figure are very satisfactory, particularly the arm segmentation in the first axial slice given the low contrast of the image. This is due to the usage of a full body model, which exploits observations in well imaged areas to influence the segmentation in others. However, some slight misalignments can be observed as well. First, there is an inherent tradeoff between data fidelity and regularization; high values of α c in Eq. 2 and Eq. 9 effectively give more importance to be coherent with the body model than to be aligned with large gradients. This can result in good segmentations in challenging areas (like the arms) but not fully accurate in other ones. Also, as has been discussed in previous work (Kainmueller 2015), we believe part of this problem is the low resolution of the triangulated mesh, as can be observed in the last axial slice. This could be changed at the expense of higher computational costs. Figure 3 shows the performance of the proposed model across different subjects with varying BMIs. As depicted in the coronal slices in the top row of Figure 3, the delineated geometries are smooth despite the large slice gaps of the MRI dataset. This is because of the regularization of the shape model. The bottom row of Figure 3 depicts the corresponding axial slices of the segmented geometries using the statistical shape model and the manually labelled slices of the subcutaneous fat layer of the mid-section of the different dataset. Table 4 shows the accuracy of the segmentation of the exterior and interior geometries of the SCAT of the different MRI datasets. The Dice similarity coefficient (DSC) (Dice 1945;Polak, Zhang and Pi 2009) which measures the overlap between the manually labelled voxels and the segmented voxels is computed for the segmentation. The sensitivity, specificity and overall accuracy are also measured to evaluate the performance of the segmentation technique. The overall accuracy denotes the measure for the correctly segmented voxels based on the number of voxels in the image. In particular, the accuracy of the segmentation is defined as accuracy = T N +T P T N +F P +T P +F N , in which T N/T P denote the true negative/positive voxels, and F N/F P denote the false negative/positive voxels. The sensitivity and specificity are computed as the percentages of foreground and background voxels which are correctly segmented as foreground and background voxels, and are defined as sensitivity = T P T P +F N and specif icity = T N T N +F P . As shown in Table 4, the proposed model achieved a high accuracy in the segmentation of the exterior and interior SCAT geometries.   Figure 4 depicts the subcutaneous fat thickness of the different datasets represented as a color map on the segmented exterior surface. There is a substantial variance in the fat thickness, both in the quantity and distribution of fat in the different geometries. The proposed method, applied to a large set of MRI scans, could be used for obtaining a statistical model of the SCAT distribution across population.

Conclusion
A method for segmenting and registering the subcutaneous adipose tissue from MRI scans has been proposed. The proposed algorithm incorporates a statistical full-body model, learned from thousands of 3D surface scans, to regularize the segmentation objective and establish correspondences between the subcutaneous fat layer across subjects. The factorization of the statistical model into pose and shape enables the segmentation of subjects with varied body shapes and different poses. The qualitative results on MRI scans of subjects with varied BMIs are presented, which show the potential of the proposed algorithm for analysing the statistics of subcutaneous fat in a large corpus of MRI scans.