## RESEARCH ARTICLE

# Quantification of Epicardial Fat by Cardiac CT Imaging

**Giuseppe Coppini**

^{1}**,**

**Riccardo Favilla**

^{1}**,**

**Paolo Marraccini**

^{1}**,**

**Davide Moroni**

^{*, 2}**,**

**Gabriele Pieri**

^{2}^{2}Institute of Information Science and Technologies (ISTI), Italian National Research Council (CNR), Pisa, Italy

### Article Information

#### Identifiers and Pagination:

**Year:**2010

**Volume:**4

**First Page:**126

**Last Page:**135

**Publisher Id:**TOMINFOJ-4-126

**DOI:**10.2174/1874431101004010126

#### Article History:

**Received Date:**5/12/2009

**Revision Received Date:**24/2/2010

**Acceptance Date:**1/3/2010

**Electronic publication date:**27/7/2010

**Collection year:**2010

**© Coppini**

*et al*.; Licensee*Bentham Open*.open-access license: This is an open access article licensed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted, non-commercial use, distribution and reproduction in any medium, provided the work is properly cited

## Abstract

The aim of this work is to introduce and design image processing methods for the quantitative analysis of epicardial fat by using cardiac CT imaging.

Indeed, epicardial fat has recently been shown to correlate with cardiovascular disease, cardiovascular risk factors and metabolic syndrome. However, many concerns still remain about the methods for measuring epicardial fat, its regional distribution on the myocardium and the accuracy and reproducibility of the measurements.

In this paper, a method is proposed for the analysis of single-frame 3D images obtained by the standard acquisition protocol used for coronary calcium scoring. In the design of the method, much attention has been payed to the minimization of user intervention and to reproducibility issues.

In particular, the proposed method features a two step segmentation algorithm suitable for the analysis of epicardial fat. In the first step of the algorithm, an analysis of epicardial fat intensity distribution is carried out in order to define suitable thresholds for a first rough segmentation. In the second step, a variational formulation of level set methods - including a specially-designed region homogeneity energy based on Gaussian mixture models- is used to recover spatial coherence and smoothness of fat depots.

Experimental results show that the introduced method may be efficiently used for the quantification of epicardial fat.

**Keywords:**Image segmentation, level set methods, epicardial fat, cardiac CT.

## 1. INTRODUCTION

Epicardial fat, as other visceral fat localizations, is correlated with cardiovascular disease, cardiovascular risk factors and metabolic syndrome [1, 2]. The important characteristic of epicardial fat is the anatomical localization and its functional relationship with coronary arteries and myocardium [3]. It is postulated that epicardial fat may work as an endocrine organ secreting hormones, cytokines and chemokines that may play a role in the atherogenesis [4, 5]. This hypothesis is partially supported by few epidemiological studies. Moreover many concerns remain about the method for measuring epicardial fat, its regional distribution on the myocardium and the accuracy and reproducibility of the measurements. At present, three techniques appear suitable for the quantification of epicardial fat, namely echocardiography, Magnetic Resonance Imaging (MRI) and Computed Tomography (CT). All of them have been used in medical studies (see e.g. [6-9]).

With respect to other imaging modalities, however, CT may provide a more accurate evaluation of fat tissue due to its higher spatial resolution compared to ultrasound and MRI.

In addition, CT is widely used for evaluating coronary calcium score, which is an independent predictor of cardiovascular prognosis [10], and to obtain a non-invasive study of coronary morphology. Current MultiDetector Computed Tomography (MDCT) systems allow the visualization of the entire heart volume during an entire heart cycle with high density resolution. MDCT imaging is thus opened to a wide spectrum of cardiac applications. As pointed out by other researchers [11], different pieces of information about coronary morphology, myocardium structure and function may be simultaneously collected leading to a more complete assessment of the functional status of the heart. Quantitative evaluation of epicardial fat may add a prognostic value to cardiac CT examinations with a potential improvement of its cost-effectiveness. Unfortunately, evaluating epicardial fat by direct user interaction is a long and tedious task prone to strong inter- and intra-observer variability. As a matter of fact, accurate and reproducible measurements on reasonably large populations are heavily hampered.

Although the quantitative analysis of epicardial fat may thus have important clinical implications, the medical imaging literature on this topic is still limited. Actually, the extraction of clinical parameters for the evaluation of epicardial fat poses several non-trivial problems related to the areas of image segmentation, image reconstruction and shape modeling for the geometrical and densitometric characterization of the 3D models [12] representing fat depots.

Although the normal CT attenuation range for fat tissues is known to be in the interval -190<HU< - 30 [13], such information alone is not sufficient to obtain accurate segmentation of fat depots. Indeed, such interval is somewhat dependent from the CT scanner and varies also across individuals. For this reason, the attenuation range for fat should be assessed on the specific study under examination.

In addition, care must be taken in using the attenuation range for threshold-based segmentation, since spurious regions should be previously eliminated.

For this reason, in [14] a method is proposed for the segmentation of abdominal fat tissue from CT images. The method classifies each pixel into three classes, corresponding to *fat*, *non-fat* and *background*, using a feature-vector paradigm. In particular, a wide range of Laws' features and Gabor texture features are extracted. After assessing the discriminative power of each feature, a restricted subset was chosen and employed for segmentation through a hierarchical multi-class fuzzy affinity method.

Very similarly, in [15], the work in [14] is extended to the quantification of pericardial fat. The authors added a stage for the ad hoc labeling of organs and tissues in the cardiac CT (namely lungs, heart, descending aorta and subcutaneous regions). A major limitation of the approach is that it is mainly 2D-based and does not provide any insight in the regional distribution of the pericardial fat.

In [16, 17], a computer-aided method for the quantification of pericardial fat using non-contrast CT images is presented. The method uses a preprocessing step to remove all other structures apart from the heart, using a region growing segmentation strategy. Then, an experienced user is required to scroll through the slices between upper and lower heart limit and, if the pericardium is visible, he is required to place from 5 to 7 control points. Catmull-Rom cubic spline functions are then automatically generated to obtain a smooth closed pericardial contour. Identification of the fat inside this contour is finally achieved by thresholding.

The recent contribution [18] uses the same preprocessing step presented in [16], but adds a high level step for the identification of the pericardium, thus avoiding manual placement of spline control points. Indeed, in each slice, the segmentation algorithm sweeps the anterior region of the heart from 0 to 180 degrees with rays having the heart centroid as origin. The higher intensity points along each ray are stored and used to reconstruct the epicardial contour, employing some robust estimation method to deal with the quite strong amount of noise and outliers. Although this approach may be interesting (and, indeed, may be integrated with the methods presented in this paper), the reported results are not impressive (with only 4 out of 10 images correctly segmented in a fully automatic way) and, so, some improvement is still necessary.

For what regards the estimation of visceral fat by other imaging modalities, in [19], an approach for abdominal fat estimation from MRI data is presented. The approach is based on two point Dixon imaging technique and consists in three stages for a) 3D phase unwrapping, b) image intensity inhomogeneity correction and c) final morphon-based registration and segmentation of fat tissue.

With respect to the related works presented above, the aim of this paper is to focus on image processing methods for the quantitative analysis of epicardial fat from volumetric non-contrast CT data acquired for calcium scoring, by extending the previous work [20]. Particular care is given throughout the paper to the identification of possible sources of inter- and intra-observer variability. Besides the description of a procedure for identifying the Volume of Interest (VoI) in the CT scan, the main contribution of the paper is a two step segmentation method suitable for the analysis of epicardial fat. In the first step of this method, an analysis of epicardial fat intensity distribution is carried out in order to define suitable thresholds for a first rough segmentation. In the second step, an advanced segmentation method, based on level set, is used to recover spatial coherence and smoothness of fat depots. Among several frameworks for level sets, the geodesic active contour framework [21] has been selected, exploiting a variational formulation without contour re-initialization. Since intensity and region homogeneity are crucial for the classification of voxel belonging or not to fat tissue, a specially-designed region homogeneity term based on Gaussian mixture models has been incorporated in the framework. Experimental results show that the introduced segmentation method may be efficiently used for the quantification of epicardial fat.

The paper is organized as follows. In Section 2, the steps leading to an user-assisted segmentation of epicardial fat and to the extraction of proposed parameters for its quantification are described, while in Section 3, after describing the considered acquisition protocol, the results of the processing are reported step by step. Section 4 ends the paper with remarks and directions for future work.

## 2. METHODS

The procedure that has been devised for the measurement of epicardial fat consists in four steps, namely Identification of the VoI (Section 2.1), Segmentation of Epicardial Fat (Section 2.2), Segmentation Refinement (Section 2.3) and, finally, Extraction of Clinical Parameters (Section 2.4).

### 2.1. Identification of the Volume of Interest (VoI)

The goal of this section is to identify the Volume of Interest (VoI) on which to perform the quantitative analysis of epicardial fat. In particular, the first stage is to identify the anatomical structures in the CT scan corresponding to lungs, liver, diaphragm, sternum and ribs, in order to remove them from the VoI, since they do not contribute to epicardial fat. In addition, since the basal part of the heart (i.e. the volume containing the atria and the great vessels) is known to be more complex, its inclusion for the analysis of epicardial fat would require either strong manual intervention (leading to scarcely reproducible results) or extremely refined and knowledge-based segmentation methods. Thus, for the scopes of the present work, we focus on the analysis of the epicardial fat localized below the atrioventricular sulcus, i.e. in the regions corresponding to left and right ventricle. Indeed, the atrioventricular sulcus is an easily discernible pseudo-planar anatomical feature that can be certainly identified by experts, with very low inter-observer variability.

It is thus assumed that i) the acquired CT scan has been preprocessed by an expert, ii) the atrioventricular sulcus has been identified and iii) the volumetric dataset has been resampled along planes parallel to the sulcus. Formally, denoting with *I* = *I* (*x*, *y*, *z*) the volumetric dataset (where (*x*, *y*, *z*) ranges through the domain, namely a rectangular box ), it is assumed the the planes with equation *z* = const are parallel to the intraventricular sulcus.

Under this assumption, some assisted image editing is still necessary for the identification of the VoI. Indeed the region corresponding to the heart should be separated from the nearby tissues. Since in computed tomography, the attenuation coefficient of lung tissue is well separated from the one corresponding to muscles and fat, lungs do not pose any problem, since they can be disregarded by simple thresholding. On the other hand, the intensity values of liver and diaphragm are similar to values found in the VoI. Thus, they cannot be separated merely on the basis of the intensity value. Another similar problem consists in the separation of fat depots near the chest wall from the epicardial fat. Also in this case, the separation cannot be based just on the local image appearance, but some complex vision system is required.

For the scopes of the present work, a simple form of manual intervention is introduced for these separation problems. Namely, an expert observer is required to scroll the slices between the atrioventricular sulcus and the apex and to place some control points on the pericardium. Such control points are automatically used to draw a natural cubic spline, which represents the closed pericardium contour, on each slice. A rasterization procedure is finally applied to convert the stack of splines into the volume V ⊂ B they bound. The volume V represents the identified VoI.

### 2.2. Segmentation of Epicardial Fat

After having defined the VoI, a first rough estimation of epicardial fat is obtained by applying a threshold-based approach.

As reported in Section 1, absolute Hounsfield Unit (HU) values of voxels correspond directly to tissue property and, thus, epicardial fat may be identified by thresholding the VoI with the interval corresponding to adipose tissue (i.e. -190<HU<- 30 ). However, such nominal interval should be corrected on the basis both of the actual CT scanner and of the acquisition protocol (especially in case of contrast-enhanced images). In addition, it has been reported that the attenuation range varies also across population [15]. Thus, a tool has been designed for guiding the selection of a suitable interval [*T*_{0} , *T*_{1}] to achieve the segmentation of epicardial fat, expressed as the following binary image *J* :

(1) |

In particular, if the segmentation obtained by the default interval ([*T*_{0} , *T*_{1}] =[-190, - 30], using values in HU) is not satisfying, the cardiologist may proceed to the fine tuning of the the interval, for example by “probing” a patch of the image corresponding to epicardial fat. Indeed, the analysis of the histogram of the patch may be used to define a custom interval [*T*_{0} , *T*_{1}]. More precisely, the following procedure is used:

- Select a voxel
*P*= (*x*_{o},*y*_{o},*z*_{o}) in the interior of the region corresponding to epicardial fat. - Consider a cube
*C*_{o}centered in*P*with edge length ℓ. - Compute the mean
*µ*and the standard deviation*σ*of the volumetric dataset*I*restricted to the cube*C*_{o}. - Define the thresholds in the interval [
*T*_{o},*T*_{1}] by:

where*ξ*is a positive real number. Notice that the interval is chosen with center*µ*and width proportional to*σ*.

### 2.3. Segmentation Refinement

The thresholding operated by Equation 1 provides a segmentation of epicardial fat which may suffer from spurious oscillations in the contour. A level set method is used to achieve both contour regularization and better adherence to image edges. Among several frameworks for level sets, we work slice-by-slice in the geodesic active contour framework [21], exploiting a variational formulation without contour re-initialization. Since the ultimate goal is to classify voxels inside the VoI with the labels fat and non-fat on the basis of their intensity, regional information is crucial and, thus, a region homogeneity term has been incorporated in the framework. More formally, consider for each slice z = z* _{k}* = cost , the region of interest R = V ᑎ {z = z

*} and let*

_{k}*I*(

_{k}*x, y*) =

*I*(

_{k}*x, y, z*) be the induced two-dimensional image. Let {R

_{k}^{+},R

^{-}} be a partition of R into two disjoint regions corresponding to non-fat and fat respectively. Notice that the slice index

*k*has been omitted to avoid cumbersome notation. Let be

*C*the contour separating one region from the other. We look for an optimal

*C*by minimizing a suitable energy functional

*E*.

In the level set formulation, the contour *C* is represented as the zero-level set of a scalar function

(2) |

In addition, the region R^{+} and (resp. R^{-}) is precisely the set on which Ф assumes positive (resp. negative) values.

Artificial time is then introduced to evolve the contour towards a minimum value of the energy *E*. Taking as initial guess the partition given by the thresholding in Equation 1, we set:

(3) |

Thus, the initial value Ф_{0} at time *t* = 0 is set equal to the *Signed Distance Function* (SDF) from the contour. Of course several other choices for the function Ф_{0} may be used to encode the same initial contour, however -as it is well known- initializing the level set function with a SDF is suitable for achieving stability in numerical approximations.

For example, a SDF is almost everywhere differentiable (actually as we will see in Section 2.3.3), reducing numerical inaccuracy that may show up during contour evolution.

A composite energy term was associated with Ф_{t}, expressing the goodness of contour fitting to the image data and inner conditions of homogeneity and stability:

(4) |

Each of the summands in Equation 4 achieves a special effect and, for sake of clarity, they will be discussed separately in the sections below.

#### 2.3.1. The Geodesic Active Contour Term

As usual, let the *edge indicator function**g* be defined as:

(5) |

where *G _{σ}* is a Gaussian kernel with variance

*σ*

^{2}.

The first summand *E*_{GAC} in Equation 4 is defined as:

(6) |

where the integral is extended to the whole region R, *α* is a positive constant and *δ* is the Dirac delta function. The term *E*_{GAC} corresponds to the length of the contour *C* weighted by the edge indicator function. Requiring that the length of the contour is minimal is a well-known condition for achieving smoothness, since spurious contour oscillations are wiped out. Minimizing *g -weighted length* achieves instead *simultaneously* contour smoothness and adherence to image data. In particular, notice that a piece of contour crossing an homogeneous region is discouraged, since the region is characterized by high values of *g* (namely *g* ≈ 1). By contrast a piece of contour on an ideal edge costs nothing, since *g* = 0 there.

#### 2.3.2. The Region Homogeneity Term

The *E*_{RH} term is specially designed for favoring the selection of homogeneous regions. Namely, we suppose that the intensities of voxels belonging to R^{+} and R^{-} are drawn from two distinct probability density functions, being *p*^{+} and *p*^{-} respectively.

After a careful analysis of the histograms of the region of interest, it has been judged convenient to model both distributions by using Gaussian Mixture Models (GMM), with *N ^{+}* and

*N*components respectively:

^{-}(7) |

where are the mean and standard deviation of the *i*-th Gaussian distribution and
are the weights of the mixture.

Assuming that voxels are independent and identically distributed, the conditional probability of the image *I _{k}* given a partition {R

^{+}, R

^{-}} is given by:

(8) |

Applying the logarithm to Equation 8, we have:

(9) |

This in turn suggests to define the following energetic term for region homogeneity (see also [22]):

(10) |

where *β* is a positive constant and *H* denotes the Heaviside function, defined as:

(11) |

#### 2.3.3. The SDF Constrain Term

For the introduction of the last, technical summand *E*_{SDF}, recall that we chose Ф_{0} to be a SDF, for several numerical advantages. In usual level set formulation, Ф_{t} ceases to be a SDF as the time evolves. This can lead to less confident estimations of the derivatives in usual finite difference schemes and, eventually, the functions Ф_{t} may develop shocks. To avoid this, it is customary to periodically re-initialize Ф_{t} to a SDF. Several strategies exist for achieving this goal, for example by applying some iso-contour algorithm and, then, computing the SDF again or by solving the so-called re-initialization equation, i.e. another partial differential equation. However, none of this method is appealing, because it is difficult to decide when and how often to perform re-initialization. The strategy we adopted consists in softly constraining Ф_{t} to be a SDF during the whole evolution [23].

To discuss this issue, we first recall that a SDF Ψ satisfies:

(12) |

in every point having one and only one closest point and, thus, generically. It is not difficult to show that the converse also holds, i.e. any function Ψ satisfying is a SDF. This differential characterization of signed distance functions is suitable for the present purpose. Actually defining

(13) |

we have a measure of how close Ф is to a signed distance function. A positive constant γ is introduced to manage the tradeoff with the previous energy terms. During contour evolution, the *E*_{SDF} the term prevents Ф_{t} to move away from the space of signed distance functions. In this way, numerical inaccuracy are avoided, retaining the good features of a SDF, without any need of other numerical remedies.

#### 2.3.4. Energy Minimization

The minimization of the functional *E* is performed by gradient descent, using the initialization described in Equation 3:

(14) |

where is the first variation of the functional *E* By linearity:

(15) |

It is not difficult to compute these first variations, which are given by:

(16) |

(17) |

(18) |

In the numerical discretization, all these quantities are approximated by a first order standard finite difference scheme.

### 2.4. Extraction of Clinical Parameters

After having applied the previous method to each slice, the direct computation of global parameters for the evaluation of epicardial fat is straightforward. In particular, the area *A _{k}* of the region corresponding to epicardial fat may be computed in the slice

*z*=

*z*by the formula:

_{k}(19) |

where Ф is the optimal function found minimizing the functional in Equation 4 and *H* is the Heaviside function defined in Equation 11. Notice that there is no need to define a crisp segmentation, i.e. for example by extracting a binary mask from Ф . Using a mollified Heaviside function, partial volume effects near the interface may be also controlled. Let ∆ be the interslice spacing in the dataset. The total *epicardial fat volume**V*_{fat} may be estimated by the sum:

(20) |

where the summation runs over all the slices between the atrioventricular sulcus and the heart apex.

For clinical applications, it is convenient to normalize the total fat volume by the Body Surface Area (BSA), thus defining a *fat index*:

(21) |

Besides total volume, we started to analyze the regional distribution of epicardial fat depots. To this end, the volume of fat surrounding each ventricle is estimated. For this analysis, the user is asked to locate, in the original images, the plane separating the left from the right ventricle. Such a plane can be easily and reproducibly identified according to a well-defined anatomical landmark, i.e. the interventricular sulcus.

## 3. EXPERIMENTAL SECTION

### 3.1. Experimental Data

The procedure was tested on CT images obtained during standard calcium scoring CT studies. Cardiac CT examinations were performed using a GE Light- Speed 64 VCT scanner. Patients were pretreated with proponolol (2 mg i.v.) and isorbide denitrate (0.6 mg i.v.) as clinically indicated. After the scout of the thorax, the heart volume was scanned to obtain the images for calcium scoring. Slices with voxel size 0.625 × 0.625 × 2.5 mm were so obtained without contrast medium injection. Ten patients have been used for testing the procedures described in this paper. Examples of acquired images are given in Fig. (**1**).

Fig. (4). First segmentation of the epicardial fat by thresholding: ( |

Fig. (5). Intensity distribution of the non-fat region from Fig. ( |

Fig. (6). Level set refinement of the epicardial fat segmentation: ( |

Fig. (7). Views of the global segmented volume: ( |

Fig. (8). 3D reconstruction of the segmented volume. |

### 3.2. Results

After the expert has identified the plane *H* on which the atrioventricular sulcus approximately lies (see Section 2.1), the original images acquired by the CT scanner were resampled along planes parallel to *H*. The volume was discretized by a grid of size (*n _{x}* ,

*n*,

_{y}*n*) = (512, 512, 62) . Since the spatial resolution of the datasets is high, no significant artefact was introduced by resampling. Examples of the resampled volumetric dataset are shown in Fig. (

_{z}**2**).

The volume was then edited as described in Section 2.1, obtaining as a result a VoI containing all the voxels that may contribute to epicardial fat volume. Some views of the VoI are shown in Fig. (**3**), together with the natural cubic spline curve delineated by an experienced observer.

The plane separating the left and right ventricle, identified through the interventricular sulcus, is also traced with some standard user interaction.

After identification of the VoI, the threshold-based segmentation was performed. For the choice of the interval [*T*_{0} ,*T*_{1} ], a voxel in the interior of the epicardial fat region was selected. A cube *C*_{0} having edge length *ℓ* = 6,25 mm was considered, obtaining as values *µ* = - 115 HU and *σ* ≈ 35 HU. The parameter ξ was manually set to 2. As expected, it resulted that *µ* is robust with respect to the choice of the voxel and of the length *ℓ*. Fig. (**4**) shows the segmentation obtained from the slices in Fig. (**3**) by thresholding according to the determined interval [*T*_{0} ,*T*_{1}]=[-185, -45] HU.

The level set refinement method was then applied slice-by-slice. Using the initial guess provided by thresholding, GMMs with *N*^{+} = *N*^{-} = 2 modes were fitted to the data, using the Expectation-Maximization algorithm [24]. An example is shown in Fig. (**5**).

Visual evaluation of the segmentation results lead to the choice of the parameters *α,β,γ* involved in the summands in Equation 4. In particular, the following values were used throughout the study:

*α*=2 *β*=10^{-3} *γ*=0.2

A time-step ∆t = 3 was deemed suitable and 30 iterations in the gradient descent brought to the results reported in Fig. (**6**). After applying such refinement to all the slices, the final segmentation of the epicardial fat could be inspected by using the conventional views reported in Fig. (**7**), where tissues of interest have been shaded. The geometry of epicardial fat was also reconstructed applying the marching cube isosurfacing algorithm. An example is shown in Fig. (**8**). For this particular patient, a *V*_{fat} =78 cm^{3} was computed, which leads to a Fat Index =41 cm ^{3} / m ^{2}.

## 4. CONCLUSIONS AND FURTHER WORK

In this paper we have introduced image processing tools for the assisted segmentation and analysis of epicardial fat. Actually epicardial fat has been shown recently to correlate with cardiovascular diseases and metabolic syndrome and, therefore, the interest of the medical community urges to develop reproducible and quantitative analysis methodologies.

In this work, after describing methods for the identification of the volume of interest and for the segmentation, parameters for evaluation of epicardial fat volume and regional distribution have been discussed. The presented preliminary results are quite encouraging.

In the framework of this research, the introduced tools are being refined and efforts are being spent both to improve the automation of the processing steps and to enrich the quantitative description by performing a regional analysis of fat thickness.

In addition, quantitative evaluation of the obtained results is currently in progress. In particular, a phantom is being considered with the purpose both of calibrating the CT scanner and of gathering some ground-truth information.

Finally, the main aim will be to correlate the extracted parameters with other clinically significant data, in order to better understand the association between epicardial fat and cardiovascular risk factors.

## ACKNOWLEDGEMENTS

This work has been partially supported by European Project HEARTFAID “A knowledge based platform of services for supporting medical-clinical management of the heart failure within the elderly population” (IST-2005-027107).