A reproducible semi-automatic method to quantify the muscle-lipid distribution in clinical 3D CT images of the thigh

Many studies use threshold-based techniques to assess in vivo the muscle, bone and adipose tissue distribution of the legs using computed tomography (CT) imaging. More advanced techniques divide the legs into subcutaneous adipose tissue (SAT), anatomical muscle (muscle tissue and adipocytes within the muscle border) and intra- and perimuscular adipose tissue. In addition, a so-called muscle density directly derived from the CT-values is often measured. We introduce a new integrated approach to quantify the muscle-lipid system (MLS) using quantitative CT in patients with sarcopenia or osteoporosis. The analysis targets the thigh as many CT studies of the hip do not include entire legs The framework consists of an anatomic coordinate system, allowing delineation of reproducible volumes of interest, a robust semi-automatic 3D segmentation of the fascia and a comprehensive method to quantify of the muscle and lipid distribution within the fascia. CT density-dependent features are calibrated using subject-specific internal CT values of the SAT and external CT values of an in scan calibration phantom. Robustness of the framework with respect to operator interaction, image noise and calibration was evaluated. Specifically, the impact of inter- and intra-operator reanalysis precision and addition of Gaussian noise to simulate lower radiation exposure on muscle and AT volumes, muscle density and 3D texture features quantifying MLS within the fascia, were analyzed. Existing data of 25 subjects (age: 75.6 ± 8.7) with porous and low-contrast muscle structures were included in the analysis. Intra- and inter-operator reanalysis precision errors were below 1% and mostly comparable to 1% of cohort variation of the corresponding features. Doubling the noise changed most 3D texture features by up to 15% of the cohort variation but did not affect density and volume measurements. The application of the novel technique is easy with acceptable processing time. It can thus be employed for a comprehensive quantification of the muscle-lipid system enabling radiomics approaches to musculoskeletal disorders.


Introduction
Loss of muscle function, as a result of diseases such as neuropathies and myopathies on one hand, or sarcopenia, the age-related loss of muscle mass and function on the other, results in reduced function. Muscle diseases also play a role in the pathogenesis of osteoporosis, falling incidence and bone frailty and fractures. Muscle biopsy is the gold standard for muscle assessment but is invasive and only evaluates a small sample, not always representative, of the relevant muscle tissue. Two different approaches are currently available for in vivo 3D muscle and fat evaluation. The first consists of a qualitative or, at best, semi-quantitative grading of muscle structures based on a washed-out and moth-eaten muscle appearance in MRI [1][2][3] or CT [4,5] images. The second approach consists of a quantitative measurement of muscle volume or cross-sectional area (CSA) in both MRI or CT images and density measurements with CT [6][7][8][9][10] While MRI has advantages in the qualitative assessment of muscle structures due to its superior soft tissue contrast, quantitative MRI analyses are affected by complex artifacts such as proton spin inhomogeneities, which usually depend on the specific MR scanner and acquisition sequence. In contrast, CT is less affected by technical variations [11], and provides a higher spatial resolution, which is an important advantage for fine-grained measurements as 3D texture. In addition, muscle density cannot be measured by MRI.
For a radiomics [12] approach to the muscle-lipid system (MLS), a reproducible, robust and fast segmentation and quantification method is required, to exploit more diagnostic relevant information from CT images. This study specifically targeted the thigh instead of the whole upper or even combined upper and lower legs, because many CT scans originally performed to determine bone mineral density of the hip extend to the upper to mid shaft of the femur.

Tissue composition of the thigh
For a better understanding of the specific segmentation and quantification techniques described in the following sections, a brief overview of the various tissues and compartments of the thigh will be given. The outer layer of the leg is the subcutaneous adipose tissue (SAT) with the skin as an outer surface and the deep fascia (F) as an inner surface (Fig 1). The volume of interest (VOI) inside the fascia (VOI IF : IF = Intrafascia) can be separated into the femoral bone, anatomical muscles (M) consisting of muscle tissue and adipocytes within the muscle border and perimuscular adipose tissue (PAT) separating the anatomical muscles. In the present work, no attempt was made to separate individual muscles. Therefore, M will always refer to the combination of all muscles. Adipose tissue (SAT and PAT) consists of adipocytes (fat cells) containing lipids.
Muscle tissue can be further differentiated into myocytes (muscle fibers), which may contain intramyocellular (IML) [13], and extramyocellular lipids (EML). In a CT image, IML can be measured indirectly only via the CT muscle density, which will decrease with increasing IML [14]. EML are adipocytes embedded within M among the muscle fibers. EML is typically present as adipocyte clusters, which, depending on their size and the spatial resolution of the CT image, can be directly segmented because their CT value differs from that of the surrounding muscle tissue. According to reference [15] IMAT (intermuscular AT) will include both PAT and EML in the present study (Fig 1), although some authors [16] refer to PAT only.

Materials and methods
The integrated approach can be divided into four steps: (A) definition of a global VOI analysis of the thigh, (B) segmentation of the fascia to delineate SAT and IMAT, (C) segmentation of muscle and lipid VOIs within VOI IF and (D) definition of features characterizing the musclelipid distribution within VOI IF . CT scans of the proximal femur, obtained earlier [17] to measure bone mineral density in subjects with high fracture risk, were used. Subjects were scanned on top of a calibration phantom to calculate BMD from the measured CT values. The femoral bone was segmented as part of a previous study [18].

A. Analysis VOI
The thigh surface was first segmented using a threshold-based volume growing. The threshold, which was empirically selected as 70%-threshold between highest and lowest CT values, was determined from the CT value spectrum of a sphere automatically positioned on the thigh surface.
The muscle-lipid analysis was carried out at the level of the upper femoral shaft (VOI US ; Fig  2). Segmentation of the VOI US started with the automatic determination of two points: the first was the center of mass of the femur in plane A, which was perpendicular to the femoral shaft axis and intersected the mid-height of the lesser trochanter. The second point was the voxel of the femoral head at a maximal distance from A. The projected distance between the two points onto the scanner z-axis was used as an anatomic size-and pose-specific distance d. VOI US consisted of n = 0.5 d / s (s: CT slice thickness) slices. The most proximal slice was defined as the distal end of the acetabulum.

B. Segmentation of the fascia
Accurate fascia segmentation is essential for the separation of SAT from muscle and IMAT. However, fascia segmentation is difficult since it is a very thin structure with low contrast relative to the surrounding AT. Pure muscle tissue has higher CT values than AT. However, muscle tissue with high lipid infiltration has CT values close to those of AT. Muscle had also to be distinguished from blood vessels, edema, dermis and genital organs through the following steps: 1. The first step was a gross identification of muscle, based on its CT appearance: a grade G in the range [1,3] was first set by the operator where G = 3 denotes moth-eaten and washedout muscle structures (Fig 3a and 3b). In addition, a contrast value (C), relative to AT and water was assigned to each voxel (v) according to a linear scale where C = 0 for a CT value equal to the average CT value of the adipose tissue (CT AT ) and C = 1 for a CT value equal to that of water. CT AT was determined in a VOI AT resulting from simple volume growing using a CT value range , which started within the sphere on the thigh surface. CT H2O was determined from the in-scan calibration phantom.
In order to classify a voxel v as potential muscle, four empirical conditions had to be fulfilled: (1) a minimum distance of 10 voxels from the body surface; (2) connection with the femur; (3) a contrast C above the minimum contrast C min characterizing the washed-out appearance; and (4) an aggregation α above the minimum aggregation α min characterizing the moth-eaten appearance, where: Θ denotes the Heaviside function and v' is a voxel in the 26-neighborhood of the voxel v. The aggregation value is the ratio of neighboring voxels with the contrast above C min . From the resulting voxel clusters only the two largest connected components within VOI US were kept. According to the position of their center of mass relative to the femur, the two components were classified as right (VOI RC ) or left (VOI LC ). VOI RC and VOI LC were independently processed by simple morphological closing using a z-axis elongated elliptical structure Step 1: Connected components of right leg (VOI RC ) for G = 1 (a) and G = 3 (b) (G: Grade of muscle appearance, see text).
Step 2: Probable muscle voxels within processed left leg (VOI LC ) (c) and compactification resulting in VOI IFA (d and e) (IFA: Approximation of IF).
Step 3: refinement of VOI IFA . Candidate voxels close to FA that may belong to VOI IF (f and g), seed points to select clusters that belong to VOI IF (h), clusters that were not selected stay porous (i) and are removed by morphological smoothing yielding the final result (j). User-dependence is indicated by blue color.
https://doi.org/10.1371/journal.pone.0175174.g003 element to account for the cylindrical shape of the fascia and the predominant extension of the muscles along the leg axis, and contour-filling.
2. The processed VOI RC and VOI LC exclude the dermis, connective tissue and genital organs, whose CT values are typically similar to those of muscle. Within each of these VOIs, a simple contrast value threshold was used to identify probable muscle voxels (Fig 3c). These voxels were compacted (Fig 3d and 3e) resulting in a structure named VOI IFA , an approximation of the intrafascia VOI. Its surface was named fascia approximation (FA).
3. FA was already a very good approximation of the fascia when the amount of IMAT was low, i.e. in younger healthy subjects, in which the fascia is usually in direct contact with the muscle. However, with increasing adipose tissue infiltration, the amount of IMAT increases and the fascia is often bordered by adipose tissue on either side, SAT on the outer and IMAT on the inner side. In many areas, the CT contrast of the fascia is very low, preventing the direct detection of a closed surface. Instead, in our study FA was used as a start point to select additional voxels, which were likely to belong to VOI IF . Specifically, voxels, which were connected to FA, were not connected with the dermis and were located at a maximal distance of 15 voxels from FA, were identified as candidates (Fig 3f and 3g). Connection to the dermis was defined as voxels with C > 1 connected with the thigh surface. Based on volume growing, connection to FA was ensured by a local noise adaptive threshold. Namely, a voxel v was accepted by the volume growing if f was true: where g depended on the mean CT value (mean), and the standard deviation (SD) in the 26-neighborhood of v: The 2SD factor accounted for noise in the CT data. Next, starting from a seed point set by the operator in the space between FA and a local group of the new candidate voxels, 'rays' were sent out isotropically (26 directions). If 65% or more of these rays were reflected by one of the candidate voxels or by FA, the status of the voxel was changed to 'shielded' [18] (Fig 4). All shielded and candidate voxels were merged with VOI IFA (Fig 3i), which was then filtered for porous or non-compact structures by a morphological smoothing. The resulting VOI was called intrafascia (VOI IF , Fig 3j), its surface defined the 3D fascia. As shown in Fig 3h, often multiple seed points were required to process the complete dataset. The operator could also select not to add certain conglomerates of candidate voxels to VOI IF . 4. Optionally, a fascia-specific modification of the explicit deformable model by Mastmeyer [19] could be utilized to attract the fascia to the nearest maximum under constraint of a given regularity of the 3D surface. This step was useful, if neither strong edema nor very washed-out muscle structures were present in the image, which was typically the case for younger patients.

C. Muscle and lipid VOIs
In order to facilitate a comprehensive analysis of the relation between muscle and lipids, various VOIs were determined. With the exception of VOI SAT , all of them were located within VOI IF . Specifically an anatomical muscle VOI M and a muscle tissue VOI MT were defined.
1. VOI SAT was obtained by subtracting VOI IF and all voxels connected with the dermis from VOI US . The spectrum of CT values of VOI SAT (Fig 5) was used to define an AT threshold (T AT ), which differed from CT AT used for the fascia segmentation. Typically, SAT is very homogenous but may contain edema and blood vessels that have higher CT values than  Semi-automatic method to quantify muscle-lipid distribution in clinical CT images adipose tissue. Therefore high and low CT values of the spectrum with frequencies lower than 30% of the most frequent CT value were cut. In order to determine T AT from a symmetric spectrum to prevent an influence of noise on the mean value, the lower tail of the spectrum was also cut resulting in a spectrum SAT Ã . T AT was defined as mean + 2SD from the mean and the standard deviation of the SAT Ã CT value histogram (Fig 5 left).
2. VOI M was determined by subtracting VOI PAT from VOI IF under the assumption that PAT but not EML was connected with SAT. VOI PAT was defined by connection with the fascia under the condition that the CT value of the voxel under consideration was smaller than T AT .
3. For the definition of VOI MT , a Gaussian mixture model (GMM) was employed in combination with a Levenberg-Marquard optimization algorithm [20] to fit two Gaussian curves to the CT value spectrum of the combined SAT and IF VOIs (VOI SAT [ VOI IF = VOI MLS ) (Fig 5,  right). One fit curve represented AT and one MT. For AT, height, peak and width of the SAT Ã distribution were used for initialization of the fit procedure. For MT the ratio of the number of voxels of VOI AT versus those of VOI MLS was used to initialize the height. A phantom based CT value was used to initialize its peak. VOI MT was then defined by 3D volume growing inside VOI IF starting from seed points defined as voxels with CT values higher than the peak b of the fitted MT curve. Voxels were included if their CT values were higher than that of the intersection a between AT and MT distributions (Fig 5, right).
4. Finally, in order to exploit the fact that the CT values of MT also reflect IML, similar to [21] an abstract description for the muscle-lipid distribution was formulated. For each voxel of VOI IF , a muscle concentration between 0% and 100% was determined using T AT and T HDM with T AT defining a muscle concentration of 0% and T HDM , a muscle concentration of 100% (HDM: high density muscle). It was determined as mean CT value minus the SD determined in a sphere placed in the muscle tissue of young climbers. Results from 30 subjects were averaged and resulted in a CT value of 35 HU. In order to account for calibration differences among scanners, T HDM was finally defined as 35 HU + CT H20 . Thus, apart from the constant of 35HU, it only depended on the water value of the in-scan calibration phantom. The two CT value thresholds T HDM and T AT were used to define 6 bins with decreasing "muscle concentrations". Bin 6 (B6) defined HDM for ! 100% muscle concentration, B5 muscle concentration for 75-99%, B4 for 50-74%, B3 for 25-49%, B2 for 0-24% and intermuscular adipose tissue (IMAT) for 0%.

D. Feature extraction
Numerous image features can be calculated in the segmented VOIs to describe the musclelipid distribution. Only few features quantifying volume, density and structure will be discussed in the present report to investigate the reproducibility and stability of the segmentation and the feature extraction.
Densities were determined either from the segmentation-based or from GMM-based procedures. Segmentation-based densities were defined as average CT values in a VOI after subtraction of CT H2O . GMM-based densities were defined as peaks of IF or MLS distribution curves again after subtraction of CT H2O i.e. D IF GMM (MT) or D MLS GMM (MT). rV defined a volume of a VOI relative to VOI MLS or VOI IF . A more advanced 3D descriptor was the average grain size G avg as determined by granulometry [22]. Porous and atrophic MT and HDM should have a smaller G avg . Fractal dimension (FD) [23] was used as a measure of 3D texture roughness. The 3D texture of dystrophic muscles should have a rougher appearance. Further, sphericity C of M and MT was quantified. C should be smaller for irregular surfaces and elongated muscle shape in case of pathologies.

E. Validation
Validation of the segmentation stability and features included determination of the intra-and inter-operator precision errors and the impact of noise and of CT H2O derived from the in-scan calibration phantom.
25 datasets from elderly subjects (age: 75.6 ± 8.7) of the EFFECT study [30] with a washedout and moth-eaten muscle appearance and high levels of IMAT and edema were processed. In these subjects, FA differed largely in most cases from the true fascia and an operator-interaction was needed in such cases. For the precision analysis, all 25 datasets were analyzed three times on different days by the same operator for intra-operator-precision and once by three different operators after an initial training for inter-operator-precision. Precision errors were calculated as CV RMS values [24]. In order to assess the diagnostic value of a feature, the CV RMS was compared to the biological variation of the 25 subjects measured as percentage coefficient of variation divided by a factor of 100 (CV 100 ).
Gaussian noise was added to the CT datasets acquired at 120 kV with 170 mAs to simulate a lower exposure of 135 mAs. Finally, CT H2O obtained from the in-scan calibration phantom was changed by ± 5 HU.

Fig 6 shows segmentation results
for three different patients, in which muscles were heavily infiltrated by adipose tissue and fascia not always in direct contact with the muscle. Thus, segmentation required operator interactions. However, even in these challenging cases the hierarchical segmentation process limited the required operator interaction to a few clicks by setting some seed points as shown in Fig 4f. Depending on the segmentation complexity the complete processing time per CT dataset was around 12 min (i5 processor 5GHz, 4GB RAM). Table 1 shows results for intra-and interoperator analysis precision errors. There were no repeat scans of patients. For comparison, Table 1 also shows the CV 100 of the same feature from the 25 subjects. All CV RMS results were below 2% and even below 0.6% in most cases reflecting the low impact of the operator interactions. Most of the CV RMS results were comparable to the CV 100 value.  Table 2 shows the percentage changes after simulating an exposure level of 135 mAs, which would double the noise when compared to the 170 mAs exposure, which was actually used for the CT acquisition. Density and volume measurements of VOI SAT and VOI IF were not affected by higher noise. For most other features, the noise related change was about one magnitude higher than the reanalysis precision error at 170 mAs. Highest noise related changes were measured for the volume of VOI IMAT relative to the volume of VOI IF and for the average grain size of VOI MT For the features or segmentations underneath the fascia that were calibrated by water calibration or more specifically CT H2O , minor variations (± 5 HU) in the CT H2O value, which in our case was obtained from the in-scan calibration phantom, caused a much larger effect than a increase in noise. Apart from the density of HDM, a small change (± 5 HU) in the CT H2O  Semi-automatic method to quantify muscle-lipid distribution in clinical CT images value caused density differences of more than 10% (Table 3) although relative to CV 100 results were comparable.

Discussion
To our knowledge, this is the first study presenting a hierarchical 3D approach to segment and quantify the different muscular and lipid components of the thigh. One critical step is segmentation of the fascia, which separates SAT from IMAT and muscle from edema, blood vessels, genital organs etc. Calibration based on the relatively homogeneous subcutaneous adipose tissue and the CT value of water, which, in the present work, was determined using an in-scan calibration phantom, is another important characteristic of our method. The third critical step is the attempt to characterize muscle in different manners: as anatomical muscle, as muscle tissue, and in an abstract way-comparable to grey matter concentration in neuroimaging [25]as muscle concentrations. The most difficult step was the segmentation of the fascia covering the muscles, which, in elderly subjects, is no longer in tight contact with the muscle surface. The main purpose of most of the preparatory steps illustrated in Fig 3a-3i was the identification of a search space for voxels likely to belong to VOI IF . A second purpose was the differentiation of external genital organs, vessels, dermis and edema from the moth-eaten and washed out muscle as shown in Fig 3. CT value differences, for example between edema and washed out muscle structures, are small and therefore techniques based on a simple global threshold to separate muscle and adipose tissue will have severe limitations. Our method based on fascia segmentation and integration of anatomical knowledge does not have such inconvenients (Fig 6).
A limitation to our study is the lack of comparison to a gold standard such as manual segmentation. However, manual segmentation is highly fastidious and therefore almost not feasible in 3D and not used by radiologists who prefer to perform semi-quantitative grading based on gross muscle CT appearance [1][2][3]. We also have not applied our method to patients with diseases that largely destroy muscle tissue. In this case, the separation of adipose and muscle tissue as shown in Fig 5 may fail. Further, we have not analyzed any longitudinal data in order to quantify the effect of age, disease or treatment related morphological changes on segmentation. Finally, we have not analyzed the impact of spatial resolution on segmentation. CT images used here had a slice thickness of 1 mm and an in plane voxel size of (0.8 mm) 2 and were reconstructed with a medium kernel. This provides an adequate balance between spatial resolution, noise and radiation exposure. Different settings of CT acquisition and reconstruction parameters would be required to evaluate the impact of spatial resolution.
The very low intra-and inter-operator precision errors of our method, below 1% (Table 1) is due to the high degree of automation of the analysis. Manual interactions as those shown in (Fig 3h) and pre-grading of muscle structures (Fig 3a and 3b) had minimal impact on precision. However, we did not carry out a precision study including patient repositioning, which Semi-automatic method to quantify muscle-lipid distribution in clinical CT images may increase precision error, since muscle shape may change across multiple scans. Repeating CT acquisitions in humans is difficult due to ethical issues. Precision error has to be put in perspective with the changes to be measured, either longitudinally to assess changes in a given subject or cross-sectionally to compare different subjects. Table 1 shows that the reanalysis precision errors were comparable with 1% of the variation of the corresponding measurement in 25 elderly female subjects, again a very good result.
It remains speculative whether other segmentation methods [26][27][28][29] would improve the results obtained here. When applying a hierarchical approach, selection of the different steps and their order remains subjective. To our knowledge, starting with a segmentation of the fascia, an anatomical boundary at the superficial aspect of the muscles is an innovative method. However, a comparison of the different segmentation methods remains to be done.
Existing techniques for muscle measurements on CT images are restricted to measure muscle volume or area and density. Several authors used fixed CT value thresholds for muscle segmentation [30][31][32]. Use of fixed thresholds to carry out muscle segmentation does not address differences in muscle density caused by variable degrees of muscle lipid infiltration, which vary widely among muscles.
Use of fixed CT values or Hounsfield units also does not address variation with scanner type and manufacturer. Although all clinical CT scanners are routinely calibrated to water, deviations of ±10 HU are frequent. In addition, a second value is required for calibration of muscle and adipose tissue. In the field of osteodensitometry [33], specific in-scan phantoms are used for the calibration of CT to measure bone mineral density.
In the method developed here, calibrated CT values were used for two different purposes. In the segmentation process CT H2O , the CT value of the water insert of the calibration phantom and the mean CT value of adipose tissue were used to define a contrast scale, which was primarily used to define probable muscle voxels and to define upper and lower limits for the local adaptive threshold based ray reflection model. However, the final fascia resulted from the combination of several advanced image-processing procedures rather than from simple thresholding.
For segmentation of VOIs inside the fascia, two CT values derived from image information were used, CT H2O and a second CT value defined precisely using the SAT CT value histogram (Fig 6). SAT is one of the most homogenous adipose tissues in the human body [34], and is well-suited for this purpose but, in order to exclude higher density voxels due to vessels and edema contained into the SAT, the histogram was trimmed, which resulted in a unique subject-specific threshold T AT . Finally, the threshold for 100% muscle depended only on CT H2O and a constant determined from muscle tissue of a group of young athletes.
Thus, our segmentation method uses a calibration that depends only on subject-specific internal CT values and CT H2O and is therefore independent from the scanner model and manufacturer. An accurate determination of CT H2O is critical: as shown in Table 3, a change of ± 5 HU in CT H2O results in a change in density values of up to 30%, emphasizing that deviations from the regular water calibration of clinical CT scanners must be controlled.
Another critical characteristic of segmentation and feature extraction is the sensitivity to noise. To test our method, we simulated a 100% noise increase by retrospectively adding Gaussian noise to the CT images. Table 2 shows that density and volume measurements of SAT were not affected by a 100% noise increase, which is very important with regard to the calibration. With the exception of rV IF of V IMAT , effects on muscle density and volume measurements were below 5%. Not surprisingly, an increase in noise had larger effects on structural measurements like 3D texture features. Automatic exposure control techniques, nowadays often used in clinical CT [35], can approximately guarantee similar noise levels between patients. A new data-driven method to identify and reduce the impact of technical variation-such as noise-on features is currently developed in our group for e.g. retrospective analyses of studies without automatic exposure control techniques.

Conclusion
We have developed an innovative semi-automatic approach for segmentation and quantification of the muscle-lipid distribution in CT images of the thigh. An important characteristic is the calibration combining subject-specific internal subcutaneous adipose tissue with an externally derived water calibration. Another characteristic is the hierarchical segmentation of the fascia consisting of a pre-grading based on established diagnostic criteria and a local adaptive threshold-based 3D ray reflection model which requires only minimal user interaction and offers an automatic refinement by an explicit deformable model. Extracted features included density-and volume-based muscle and lipid measurements as well as advanced 3D features for a detailed quantification of the muscle-lipid distribution underneath the fascia. The precision of the method was excellent and processing speed acceptable, enabling a comprehensive radiomics approach to musculoskeletal lipid distribution disorders from standard-of-care CT images. Results reported here concern the thigh but the same concept may be applicable to other body parts although some anatomy specific changes will be required in the segmentation step.