Automatic Four-Chamber Segmentation Using Level-Set Method and Split Energy Function

Article information

Abstract

I. Introduction

Computed tomography (CT) is an imaging procedure that uses X-ray technology to produce tomographic images of specific objects. To produce CT volume of data, the patient is placed inside a tube. This tube emits X-rays toward the center of the cylinder. The X-rays pass through the body, and the intensity is measured on the other side. Then, reconstruction is performed to actually obtain a three-dimensional (3D) image. CT distinguishes bones better than organic tissues. The muscles and cavities of specific organs are not well differentiated, both appearing in similar gray tones on the CT scan.

The heart is a muscular organ about the size of a closed fist that functions as the body's circulatory pump. It takes in deoxygenated blood through the veins and delivers it to the lungs for oxygenation before pumping it into the various arteries, which provide oxygen and nutrients to body tissues by transporting the blood throughout the body. The heart is located in the thoracic cavity medial to the lungs and posterior to the sternum.

Computed tomography angiography (CTA), magnetic resonance imaging (MRI), and single-photon-emission computed tomography (SPECT) are widely used, non-invasive cardiac imaging modalities. Compared with other modalities, CTA can provide more detailed anatomic information of the heart chambers, vessels, and coronary arteries due to its higher spatial resolution [1]. Therefore, CTA is widely used in image segmentation [2] because it presents more detailed anatomic information of organs. Cardiac CTA is commonly employed for cardiac diagnosis. The disorders of the blood vessels of the heart often cause cardiovascular diseases [3], and heart segmentation from CTA has been used to detect them.

Several methods for automatic heart segmentation have been proposed. Rousseau and Bourgault [4] presented a heart segmentation method that uses an iterative Chan-Vese algorithm [5]. They used L1 fidelity term for computational efficiency instead of L2 fidelity which is a classic term. However, this approach extracted only the whole heart. Van Assen et al. [6] proposed an automatic segmentation method using a 3D active shape model in cardiac MRI. In that method, relative gray-level differences instead of absolute gray values are used for the classification of the 3D regions of interest. However, this method requires further generalization of the fuzzy inference system [7], and it segments only the LV and myocardium. Mitchell et al. [8] developed a model-based segmentation method using a 3D active appearance model. This method does not require any interactive user input of cardiac MRI and echocardiographic temporal image sequences. This method does require a number of training data sets from manually traced segmentation results. Ecabert et al. [9] proposed automatic segmentation of four chambers using a statistical geometry model and training meshes from cardiac CTA images. This method requires well-defined training data sets, and it takes too much time and effort to generate a template mesh. Funka-Lea et al. [10] developed a method to automatically isolate the outer boundary of the entire heart using balloon expansion and the graph-cut technique. This approach extracts only the whole heart, not four chambers. Grady et al. [11] evaluated three techniques, graph cuts, isoperimetric minimization, and random walks for the segmentation of four chambers. However, these methods extract four chambers with the input of seed points for each chamber, and a lot of seed points are required for accurate segmentation results.

This paper presents a method that automatically extracts four chambers by extracting the whole heart, separating the left and right sides of the heart, and splitting the atrium and ventricle regions in cardiac CTA efficiently. First, we smooth the images by applying filters to remove noise. The heart regions in cardiac CTA are inhomogeneous due to the irregular characteristics of angiography, making segmentation more difficult. Therefore, the process of image smoothing is necessary for more accurate segmentation of the heart. Second, the volume of interest (VOI) of the heart is detected using k-means clustering [12]. In this step, we utilize the intensity information of the image to extract the whole heart coarsely, and it is used for seed volumes in the next step. Third, we detect seed volumes using a geometric analysis based on anatomical information and separate the left and right sides of the heart with power watershed. Because the left and right sides of the heart have no distinct separator, it is difficult to split the left and right regions from the whole heart. To separate the heart into the left and right regions, we exploit power watershed [13] which is powerful to extract a region from an object that has weak edges or no edge. Finally, we refine the left and right sides of the heart using Malladi's level-set method [5] for more accurate segmentation, and we extract the atrium and ventricle from the left and right sides of the heart using the split energy function based on the shape and intensity of the heart. This function allows the heart to be robustly divided into chambers even in the case of having a weak edge.

The remainder of this paper is organized as follows. The next section describes the proposed method of automatic segmentation of four chambers in CTA. This procedure consists of four processing steps. Section III presents the results of the proposed method obtained using a clinical dataset. In Section IV, we summarize and discuss the results.

II. Methods

1. Pre-processing

First, we smooth the input CTA. In general, there is much noise in cardiac CTA, and it may be unclear. Therefore, image smoothing is essential to segment the heart region. Although there are many methods of image smoothing, we use Gaussian filtering [14] because the heart has inhomogeneous regions that nearly becomes homogeneous by applying Gaussian filter. Then, we use anisotropic diffusion filtering [1516], which minimizes total variation (TV), to preserve the edge while smoothing the original image and preserve finer detailed structures in images. The equation of anisotropic diffusion filter is expressed as

minTV=∫Ωux2+uy2dxdy,(1)

where u is an image, u_x and u_y are the derivatives of u w.r.t. x and y, respectively. This method reduces image noise without removing important parts, such as edges or features comparison with isotropic filters, e.g., Gaussian filter. To minimize Equation (1), we use the Euler-Lagrange equations.

4. Separation of Left and Right Heart Regions

In this step, we split the heart into the left and the right regions from the seed volumes of the previous step. It is hard to divide the heart into the left and the right regions automatically because there is ambiguity in the boundaries or difference between them. Therefore, we exploit power watershed [13].

A power watershed [13] is a segmentation algorithm used to solve an energy minimization problem based on graph theory. In graph-based image segmentation, an edge weight is assigned to encode nodal similarity such that edges with high weights are considered to be strongly connected, and edges with low weights are regarded as weakly connected or nearly disconnected nodes. One common equation for generating weights from image intensities is expressed as

(3)ω_ij = exp(−β(∇I)²),

where ∇I is the normalized gradient of the image I, and β is a weighting coefficient. Here, β was experimentally determined as 1.

For image segmentation in two classes, given foreground FO and background BA seeds, the energy model of the power watershed is given by

x=arg minx∑eij∈EωPijxi-xjq+∑viωPFixiq+∑viωPBixi−1q,(4)

s.t.x(FO) = 1, x(BA) = 0,

si=1 if xi≥12,0 if xi≤12,

where x represents the target configuration, and p and q are values that provide a reference for the different algorithms. This model is represented as a power watershed in p → ∞ and q = 2 [17]. Therefore, we set the foreground seed volume as the candidate region of the LV (left vectricle) that was detected in the previous step, and we set the background seed volume as the candidate region of the RV (right vectricle) that was also detected in the previous step. Then, we can separate the left (i.e., LV and left atrium [LA]) and right (i.e., RV and right atrium [RA]) parts of the heart using the power watershed framework. To apply the power watershed algorithm to extract the left and right heart regions, we detect the seed volume of each side of the heart. We set the seed volumes to detect the left and right regions by comparing x coordinates of the whole heart mask volume.

5. Fine Segmentation of Left and Right Heart Regions Using Level-Set Method

In this step, we segment the left and right sides of the heart accurately using the previous rough segmentation results with Malladi's level-set method [5]. We use the 3D surface and the cross-sectional curve, which are called the level-set function ϕ and the zero level-set, respectively, to extract the boundary of heart. All points that lie in the zero level-set x(t) have a zero level-set function value. Therefore, the following equation holds:

(5)ϕ(x(t),t) = 0.

Taking the time derivative of both sides by applying the chain rule,

(6)ϕ_t + ∇ϕ(x(t),t) · x'(t) = 0.

Function F represents the speed moving in the direction of the normal vector away from the curve; therefore, x'(t) · n = F, where n=∇ϕ∇ϕ. When ϕ(x,t = 0) is given, the evolution equation for ϕ can be derived as

(7)ϕ_t + F|∇ϕ| = 0.

Therefore, the evolving boundary Γ(t), according to the time and speed function F, is calculated as

(8)Γ(t) = {(x,y)|ϕ(x,y,t) = 0}.

Malladi's level-set method is adopted in this paper for level-set propagation. The speed function is formed by combining the curvature speed component, which is inversely proportional to the curvature of the zero level-set curve, with the propagation speed component, F_prop, which is the intensity value of the speed image as follows:

(9)ϕ_t + F_prop(1−εK) |∇ϕ| = 0,

where ε controls the relative importance of the curvature component for level-set propagation. The zero level-set curve moves proportionally to the propagation speed, which is very slow close to high gradient areas and fast in low gradient areas. The propagation speed term makes the zero level-set curve propagate until it reaches the edges of anatomical structures. The curvature speed term regularizes the zero level-set curve by smoothing the high curvature region. Because of the normal vector, n→=∇ϕ, the curvature K is calculated as the divergence of the unit normal vector:

K=∇·∇ϕ∇ϕ=-ϕxxϕy2-2ϕxϕyϕxy+ϕyyϕx2ϕx2+ϕy23/2.(10)

The initial heart boundary acquired in the previous step is placed as the initial contour of the level-set method, which evolves according to Equation (9). The refined heart segmentation is achieved after level-set propagation. When one contour lies completely inside the other contour, we select the outermost contour, which leads to the same effect of hole filling.

6. Detection of Atrium and Ventricle from Left and Right Heart Region

In this step, we detect the atrium and ventricle from the left and the right heart regions using the split energy function. Since it is difficult to identify valves in a cardiac CTA, we find the separation plane approximately located in a valve. To find the separation plane, we first detect the rough location of the plane that separates the right side of the heart. It is easy to search for a separation point of the right side of heart because there are intermediate slices where RV and RA are separated. We examine the slices to determine whether the right side of the heart is separated or not using 2D connected component analysis. After the separation point of the right side of the heart is detected, we detect the separation point of the left side of the heart by calculating an orthogonal projection point onto the line, orienting the first principal axis, which lies on the first coordinate and has the greatest variance by the projection of the data, and passing through the center of the left side of the heart.

Although the separation points are detected, they are only approximated candidates since only some of them are the exact points of the separation plane of the left and right sides of the heart. Therefore, we accurately detect the separation planes of the left and right sides of the heart by devising and minimizing a split energy function. We propose a split energy function as follows:

(11)E_sep (p,θ,φ) = ω_areaE_area(p,θ,φ) + ω_intE_int(p,θ,φ)

where E_sep is the separation energy function, p is a voxel position, and θ and φ are the angles with the second and third principal axis, respectively. Here, ω_area and ω_int are weights (ω_area + ω_int = 1), and E_area is the energy term of the area, defined as

Eareap,θ,φ=∫ΩMHMxdx,(12)

where HMx=1,  Mx≥10,  otherwise,

where M(_x) is the mask volume of the left or right side of the heart, and Ω_M is a plane that passes through point p with its normal vector oriented toward the angles, θ and φ. Here, E_int is the energy term of the intensity, defined as

Eintp,θ,φ=∫ΩMHMxIxdx∫ΩMHMxdx(13)

where I(_x) is the intensity value of the CT image, and E_int is regarded as the mean intensity of the intercepted plane with the mask volume. We obtain the separation plane by minimizing the separation energy function. Once p, θ, and φ are detected, we obtain the normal vector and the position of the plane. This plane is used to split the left or right sides of the heart into the atrium and ventricle.

III. Results

We tested the proposed heart segmentation method using an Intel Core2 Quad desktop system with a 3.4-GHz processor and 16 GB of main memory. We implement our method using Microsoft Visual Studio 2008 without the help of other open-source libraries. We prepared 20 clinical scan datasets, each acquired from a different patient. All the datasets were obtained using a cardiac CT scanner by Siemens. The number of images per scan ranged from 183 to 285. Each image had a matrix size of 512 × 512. The voxel size was 0.36 × 0.36 × 0.70 mm. We extracted the iso-surface using the marching cube algorithm [18] and checked the computational time for each step.

Figure 1A and 1B show the automatically segmented four chambers superimposed on axial CT images obtained from the first, second, and third datasets. Figure 1C shows the 3D surface rendering results of four-chamber segmentation. LV, LA, RV, and RA are accurately segmented in both the 2D and 3D views.

Figure 1

Result of four-chamber segmentation: left vectricle (red area), left atrium (green area), right vectricle (blue area), and right atrium (yellow area). (A) Superimposed image in the upper axial computed tomography (CT) slice. (B) Superimposed image in the lower axial CT slice. (C) The 3D surface rendering image.

To validate the results of the proposed heart segmentation method, we requested a cardiologist with 10 years of clinical experience to segment the left and right heart regions manually for each dataset. The cardiologist manually drew the boundary of the left and right heart regions in axial image slices. In this way, we obtained two manually segmented regions to serve as the ground-truth for the accuracy assessment of the proposed method.

Among the sequential processing steps in the proposed method, we evaluated the accuracy of the proposed method to segment each object including the left and right heart regions based on the following four evaluation metrics:

Efp=numVauto−numVauto∩VmanualnumVmanual×100%(14)

Efn=numVmanual−numVauto∩VmanualnumVmanual×100%(15)

Evol=numVautonumVmanual−1×100%(16)

Esim=1-2numVauto∩VmanualnumVauto+numVmanual×100%(17)

Table 1 summarizes the segmentation errors of our method for 20 datasets. Each value is the rate of mean ± standard deviation of all the datasets. In all of the datasets, E_sim ranged from 3.14% to 4.04%. The average value of E_sim was 3.43% ± 0.41% for all the datasets, indicating that the average difference between manual and automatic segmentation was less than 3.3%, approximately.

Table 1

Results of accuracy assessment of proposed method (unit: %)

We compared the segmentation accuracy of the proposed method with that of CVM [20] by measuring E_sim in (9). The average value of E_sim for CVM was 12.74% ± 7.40% for all the datasets, respectively. Specifically, some regions in the RV and RA were not included using the CVM. However, these were accurately segmented using the proposed method. Our proposed method exhibited the most accurate segmentation results with the lowest value of E_sim.

Table 2 summarizes the computational time for the evaluation of the performance of the proposed method; we measured the total processing time. The processing time, averaged over multiple tests for all the datasets, was 62.76 ± 1.09 seconds. It took 3.97 ± 0.21 seconds for pre-processing, 3.24 ± 0.28 seconds for the segmentation of the whole heart, 18.82 ± 1.68 seconds for the separation of the left and right sides of the heart, and 23.45 ± 1.42 seconds for the fine segmentation of the left and right sides of the heart and 2.82 ± 1.87 seconds for the detection of the atrium and ventricle from the left and right heart regions.

Table 2

Results of processing time of proposed method

IV. Discussion

In this paper, we propose an automatic method to segment four chambers by extracting the whole heart, separating the left and right heart regions, and splitting the atrium and ventricle regions from each side of the heart in cardiac CTA efficiently. First, we smooth the images by applying filtering to remove noise, namely Gaussian and anisotropic diffusion filtering. These types of filtering reduce image noise without removing important parts, such as edges or features. Second, the VOI of the heart is detected using k-means clustering. In this step, we utilize the intensity information of the image to extract the whole heart coarsely, and it is used for seed volumes in the next step. Third, we detect seed volumes using a geometric analysis based on anatomical information and separate the left and right heart regions with the power watershed algorithm. Because the left and right heart regions have no distinct separator, it is difficult to split the left and right regions from the whole heart. To separate the heart into the left and right regions, we exploit the power watershed algorithm. Finally, we refine the left and right sides of the heart using Malladi's level-set method for more accurate segmentation, and we extract the atrium and ventricle from the left and right heart regions using split energy. This function allows the heart to be robustly divided into chambers even in the case of having a weak edge.

We are planning to improve this method to operate robustly even when the CTA has poor contrast. Other image information, such as curvature, or the texture of the heart shape would be helpful to detect the missing part of the heart. We will devise a more robust energy model of the level-set method, which uses not only region-based information but also edge-based information and the split energy function.

Article information Continued