Automated Segmentation of Real-time 3d Echocardiography Using an Incompressibility Constraint

Coronary Heart Disease (CHD) is characterized by reduced blood flow and oxygen supply to the heart muscle, resulting from the occlusion of one or more major coronary arteries. CHD remains the most prevalent cause of mortality in developed countries and represents one of the major burdens on the healthcare systems today. Approximately every 25 seconds, an American will suffer from a coronary event, and about every minute, someone will die from one. According to the 2009 American Heart Association (AHS) report (Lloyd-Jones et al., 2009}, an estimated of 16.8 million American adults have CHD (extrapolated to US population in 2009 from National Health and Nutrition Examination Survey (NHANES) 2005-2006) and about 20% of total deaths in the United States are caused by CHD.


Introduction
Coronary Heart Disease (CHD) is characterized by reduced blood flow and oxygen supply to the heart muscle, resulting from the occlusion of one or more major coronary arteries. CHD remains the most prevalent cause of mortality in developed countries and represents one of the major burdens on the healthcare systems today. Approximately every 25 seconds, an American will suffer from a coronary event, and about every minute, someone will die from one. According to the 2009 American Heart Association (AHS) report (Lloyd-Jones et al., 2009}, an estimated of 16.8 million American adults have CHD (extrapolated to US population in 2009 from National Health and Nutrition Examination Survey (NHANES) [2005][2006] and about 20% of total deaths in the United States are caused by CHD.
The advancements in noninvasive imaging techniques, such as real-time three-dimensional (RT3D) echocardiography, have enabled physicians to detect CHD in its earliest and most treatable stage. With cardiac imaging, physicians are able to evaluate essential global and local functional parameters, such as ejection fraction (EF), wall thickening, strain/strain rate, and etc. The rapid progress in cardiac imaging, however, has led to new challenges in handling of huge amount of image data involved in comprehensive functional patient studies. Manually analyzing these data sets becomes a formidable task for cardiologists, radiologists, and technicians in order to interpret the data and derive clinically useful information for diagnosis or decision support for surgical and pharmacological interventions. Also, manual analysis is subjective and therefore compromises the accuracy and reproducibility of quantitative measurements.
The abovementioned reasons have triggered a great demand for computerized techniques to automate the analysis of cardiac images. Various image-processing tasks need to be performed in order to recover diagnostically useful information, among which myocardial segmentation is one of the most important tasks. Myocardial segmentation aims to delineate the endocardial (ENDO) and epicardial (EPI) boundaries from cardiac images. Accurate segmentation of myocardial boundaries is essential for deriving cardiac global functional parameters such as ventricular mass/volume, ejection fraction, and wall thickening. It is also a prerequisite step for accurate myocardial deformation analysis.

Background
The heart is composed of four chambers: left atrium (LA), left ventricle, right atrium (RA), and right ventricle. Fig. 1(a) shows a long-axis cross-section of the heart. As shown in Fig. 1(a), the atria and ventricles are surrounded by muscle tissues called myocardium. Contraction and relaxation of the muscle fibers in the myocardium cause the pumping of the heart.  www.intechopen.com Automated Segmentation of Real-Time 3D Echocardiography Using an Incompressibility Constraint 87 Fig. 1(b) shows a short-axis cross-section of the heart. The grey region represents the myocardium. The inner surface of the myocardium is called endocardium, and outer surface is called epicardium. As the LV is a powerful blood pump for the systemic circulation, the delineation of the LV endocardium and epicardium is often the object of interest for cardiologists, particularly because it is an important step for analyzing LV anatomical structure, quantifying cardiac function, and estimating LV motion.
Automatic LV segmentation algorithms can be broadly categorized into region-based and boundary-based methods. Region-based methods exploit the homogeneity of spatially dense information, e.g. pixel-wise grey level values, to produce segmented images. Examples of region-based segmentation methods include thresholding, region-growing, Markov random field-based approaches, and etc. In contrast to the region-based methods, boundary-based methods rely on the pixel-wise difference to guide the process of segmentation. They try to locate points of abrupt changes in the grey tone images. Examples of boundary-based methods include edge-detection, Hough transform, boundary-based snakes, and etc. However, in numerous medical imaging modalities, the LV boundaries cannot be reliably and accurately detected using the algorithms which rely only on boundary or edge information. Reasons for this include significant signal loss, noise, complicated geometric shapes. These problems are ever present for acquired in echocardiography, where the boundary detection problem is further complicated by the presence of attenuation, speckle, shadows, confusing anatomical structures.
In an effort to overcome these difficulties, various sources of image information have been used in the segmentation process. Most commonly used ones include grey level distribution, gradient, texture, and phase information. The choice of which information to use depends on imaging modality, image quality, and specific applications.

Grey level distribution
Incorporation of imaging physics as prior knowledge has been proven to be useful for cardiac segmentation (Paragios, 2002;Pluempitiwiriyawe et al., 2005, Chen et al., 2008. The random intensity distribution in ultrasound images, also known as speckle, is caused by the interference of energy from randomly distributed scatters, too small to be resolved by the imaging system. When an acoustic pulse travels through tissue or any medium, backscattering from the scatters in the range cell contributes to the returned echo. This contribution to the echo from the scatters in the range cell has been treated as a random walk because the locations of the scatters are considered to be random. The backscattered echo, therefore, is complex valued with a real part X and imagery part Y . The Rayleigh distribution has been proven to be one of the most popular models in ultrasound image analysis. For example, a Rayleigh distribution was incorporated into a region-based level-set approach for the segmentation of ultrasound images (Sarti et al., 2005). Also Cohen et al. proposed to use Rayleigh distribution as an image formation model in a maximum likelihood motion estimation scheme for noisy ultrasound images (Cohen and Dinstein 2002). Steen et al. proposed to use Rayleigh distribution for anisotropic diffusion-based edge detection (Steen and Olstad 1994).
While Rayleigh distribution is extensively used in ultrasound image processing, it is not always not norm because it is valid only in the special case of a large number of randomly distributed scatters. In this condition, the speckle is called fully developed. It has been shown that the echo envelop has a Rayleigh distribution when signal noise ratio (SNR)is approximately 1.91. Deviations from such special scattering conditions are called pre-Rayleigh condition when SNR <1.91, and post-Rayleigh condition when SNR>1.91 (Tuthill et al., 1988). The K-distribution has been proposed to model the pre-Rayleigh condition by accounting for low effective scatter density (Shankar et al., 1993;Shankar, 1995), and Rice distribution has been proposed to model the post-Rayleigh condition by accounting for a coherence component due to the presence of a regular structure of scatters within the tissue (Insana et al., 1986;Tuthill et al., 1988). However, the Rician family fails to model the pre-Rayleigh condition, and K-distribution model does not take into account the post-Rayleigh condition. Generalized K-distribution (Jakeman, 1999), homodyned K-distribution (Dutt and Greenleaf, 1994), and Rician Inverse of Gaussian distribution (Eltoft, 2003) have been proposed as general models to encompass pre-Rayleigh, Rayleigh, and post-Rayleigh conditions. Unfortunately, the complex nature of these models limited their practical applications.
Recently, Nakagami distribution has been proposed as a simple generalized model to collectively represent pre-Rayleigh, Rayleigh, and post-Rayleigh distributions (Shankar, 2000). It is a two-parameter distribution with analytical simplicity, which makes it relatively easy to estimate it parameters. The probability density function (pdf) of the Nakagami distribution is given by where    is the Gamma function,  is the shape parameter, and  is the scaling parameter. For 1   the pdf reduces to Rayleigh condition. For 1   , the pdf can be described as post-Rayleigh, which is similar to Rician distribution, while for 1   , the pdf can be described as pre-Rayleigh.
As shown in (Shankar, 2000), the shape parameter and scaling parameters can be obtained from the moments of the envelop as follows      (Davignon et al., 2005).
Apart from the theoretical models mentioned above, empirical models have also been used in ultrasound image segmentation. For example, Tao et al. modeled the ultrasound speckle with a Gamma distribution and incorporated it into a tunneling descent optimization framework to overcome local minima (Tao and Tagare, 2007). Xiao et al. used a log-normal distribution for modeling speckle in breast images (Xiao et al., 2002). Qian et al. incorporated a log-normal distribution into a level-set framework to segment rat echocardiographic images with large dropout (Qian et al., 2006).

Gradient
Intensity gradient looks for intensity discontinuities between subregions corresponding to different tissue types. Intensity gradient can be computed from an image with standard differential operators (e.g. Sobel operator). A voxel is considered to be a boundary voxel if its intensity gradient is above a threshold. Gradient-based segmentation, such as such as gradient-based snake (Kass et al., 1987), Active Shape Model (ASM) (Cootes et al., 1995), and level-set approaches (Malladi et al., 1995), has been extensively used in cardiac applications (Lynch et al., 2006;Lynch et al., 2008;Assen et al., 2008;Chen et al., 2002;Corsi et al., 2002).
One limitation of intensity gradient is that it may suffer from spurious, missing, and discontinuous edges. This is especially evident in EPI segmentation because there are usually no clear and visible edges between the myocardium and background. In addition, intensity gradient is suboptimal for ultrasound images because a boundary response in ultrasound images is anisotropic and depends on the relative orientation of the transducer to the boundary. Therefore, intensity gradient information is often used in conjunction with other image features (Pluempitiwiriyawej et al., 2005;Chen et al, 2003).

Texture
Although no formal definition of texture exists, it can be loosely defined as one or more local patterns that are repeated in a periodic manner. Texture analysis has been attempted in numerous ways in medical image analysis, particularly for ultrasound images. The three principal approaches to describe texture are image pyramids, random fields, and statistical models.
The idea behind image pyramids is to generate a number of homogeneous parameters that represent the response of a bank of filters at different scales and possibly different orientations. Different filters have been proposed including Gabor filters (Xie et al., 2005;Zhan and Shen, 2006), Gaussian derivatives (Chen et al., 1998) and wavelet transforms (Mojsilovic et al., 1998;Yoshida et al., 2003). The idea behind random fields is that the value at each pixel/voxel is chosen by two-dimensional (2-D)/ three dimensional (3-D) stochastic process. Given a type of pdf of this stochastic process, one can estimate the value at a particular pixel/voxel given the values of other pixels/voxles in its neighborhood. The most commonly used random field is Markov Random Field (MRF) (Bouhlel et al., 2004). Statistical models analyze the spatial distribution of grey level values, by computing local features at each point in the image, and deriving a set of statistics from the distributions of the local features. Statistical approaches yield characterizations of textures as smooth, coarse, grainy, etc. Major statistical texture descriptors include co-occurrence matrix (Nicholas et al., 1986, Sahiner et al., 2004, Kuo et al., 2001, auto-correlation (Chen et al., 2000), edge frequency, run length, Law's texture energy, and fractal texture description. Several researchers have proposed to extract fractal texture features (Wu et al., 1992;Lee et al., 2003;Chen et al., 2005). In addition, several attempts have been made to combine multiple texture measures to improve discrimination abilities (Hao et al, 2001;Christodoulou et al., 2003;Stoitsis et al., 2006).

Phase
The local phase provides an alternative to intensity gradient to characterize structures in an image (Kovesi, 1996;Boukerroui et al., 2001;Mulet-Parada and Noble, 2000;Hacihaliloglu et al., 2008). Phase-based methods postulate that feature information is encoded at points where phase congruency is maximized, i.e. all the Fourier components are in phase. Generally, phase is estimated by quadrature filter bank. Thus, there is a link between phasebased methods and other wavelet methods. Phase-based methods are suggested to be more robust than intensity gradient for ultrasound images because intensity gradients in ultrasound images depend on the relative orientation of the transducer to the boundary. In addition, the presence of speckles and imaging artifacts might cause the variation of intensity gradient of equally significant features in the data. One limitation of phase-based methods, however, is that speckle also has its own phase signatures, and therefore an appropriate spatial scale has to be selected.

Incompressibility of myocardium
The heart is a remarkably efficient and during mechanical pump composed of complex biological materials. The main structural elements of the myocardium are inter-connected networks of muscle fibers and collagen fibers, as well as matrix that embeds them. The fibers are generally tangential to the ENDO and EPI surfaces, following a path of a righthanded helical geometry. The interstitial fluid carries only hydrostatic pressure, which in turn, is affected by the length and configuration changes of the fibers. These cause pressure gradients, which may result in the flow of the matrix. However, the fluid within the tissue is negligible for the duration of a cardiac cycle because the myocardium as low permeability. Consequently, the myocardium can be assumed to be nearly incompressible (Glass et al., 1990).
A few independent studies have been performed to quantitatively analyze the change of myocardial volume (MV) over an entire cardiac cycle. For example, Hamilton et al. performed experiments on frogs, turtles, and dogs, and discovered a relatively consistency of the MV during a cardiac cycle (Hamilton et al., 1932). Hoffman et al. used canine Dynamic Spatial Reconstructor data to study the changes in MV, obtaining a relatively constant volume that was consistent with Hamilton's findings (Hoffman and Ritman, 1987;Hoffman and Ritman, 1985). More recent work done by King et al., who used freehand three dimensional (3-D) echocardiography to measure MV and mass, showed a 1.1% difference of volume between end-diastole(ED) and end-systole (ES) (King et al., 2002). Bowaman et al. found a variation of around 5% from ED to ES by using highresolution magnetic resonance (MR) imaging (Bowman and Kovacs, 2003). O'Donnell also analyzed the myocardial volume using MR imaging, and found the difference of ED and ES to be around 2.5% (O'Donnell and Funka-Lea, 1998). The common conclusion of these studies is that the myocardial is nearly incompressible and the variation is less than 5% during a cardiac cycle. This incompressibility property of the myocardium is used an as important cardiac structure constraint which is taken into consideration in our approach, as will be detailed in Section 4.3.

General framework
In this section, we formulate our segmentation problem in a maximum a posterior (MAP) framework combining image-derived information and the incompressibility property of the myocardium. Let I be a 3-D cardiac image, in S be the ENDO surface, and out S be the EPI surface. A MAP framework that combines image information and incompressibility constraint can be expressed as follows

PS S
is the incompressibility constraint which keeps the MV enclosed by in S and out S nearly constant during a cardiac cycle. Since the myocardium is only nearly incompressible, it is more reasonable to define it within a probabilistic framework rather than imposing a deterministic constraint.

Data adherence
Region-based deformable models have been successfully applied in the segmentation of images with weak boundaries, due to their robustness (Chan and Vese, 2001). In this work, we evolve a three-phase region-based deformable model based on the statistical intensity distribution from RT3D echocardiographic images.
To evolve a region-based model, we first need to determine the intensity distribution of each region within a cardiac image. It is obvious that the entire cardiac image is partitioned by in S and out S into three regions, namely, LV blood pool, LV myocardium, and background, as shown in Fig. 2.  Fig. 2. The short-axis view of a heart. Two dotted circles are ENDO-and EPI contours, which partitioned the entire image into three regions: LV blood pool, LV myocardium, and background. The background is inhomogeneous because it contains more than one tissues, such as RV blood pool, RV myocardium, lung air, and liver.
The LV blood pool and myocardium are homogeneous, and therefore can modeled with a single pdf. In this work, we use Nakagami distribution (Shankar, 1995)  where l  is the shape parameter of the Nakagami distribution, and l  is the scaling parameter. Equation (2) describes the intensity distribution for LV blood pool when l=1, and intensity distribution for LV myocardium when l=2.
Unlike LV blood pool and myocardium, the background (see Fig. 2) because it contains more than one tissues (RV blood pool, RV myocardium, and liver), and therefore modeling it with a single distribution function would be insufficient because it contains a wide range of intensities. To handle this problem, we use a mixture model (McLachlan and Peel, 2000) to describe the intensity distribution of the background.  Fig. 3 (c). The first peak corresponds to RV blood pool and lung air, while the second peak corresponds to RV myocardium and liver.  (3)

Incompressibility constraint
It is mentioned in Section 3 that the MV is nearly constant during cardiac cycle. Therefore, we assume that the MV has a Gaussian distribution   V is the average volume, which can be calculated from the manual segmentation of the first frame. As mentioned in Section 3, the MV changes by less than 5% during a cardiac cycle.

Echocardiography -New Techniques 94
Assuming that 0 0.025V is the maximum deviation and using the three  rule, we have the following relationship: Therefore, the incompressibility constraint is encoded into a probability function which favors small variance of the MV while penalizing large deviations.
Combining Equations (1), (3), and (4), we arrive at the following optimization problem The maximization of Equation (5) where in n and out n are the normals of in S and out S respectively, and  is the propagation time step.
To implement the evolution of Equations (6) and (7), we embed the surface in S and out S in two higher dimensional functions where  is the gradient operator. Hence, the level set evolution equation for Equations (8) and (9)

Experimental setup
The ultrasound data were acquired using Philips echocardiographic system with a 4 MHz X4 xMatrix transducer. The transducer consists of 3000 miniaturized piezoelectric elements and offers steering in both azimuth and elevation of the beam, permitting real-time volumetric image acquisition and rendering (Philips, 2005). In this work, we acquired 11 sequences of canine images, and each sequence consists of 20-30 frames per cardiac cycle depending on the cardiac rate. Hence, we ran our program with a total of 286 sets of volumetric data.

Quantitative measures
To quantify the segmentation results, we used two distance error metrics and two area error metrics, namely, mean absolute distance (MAD), Hausdorff distance (HD), percentage of true positives (PTP), and percentage of false positives (PFP).
Let A and B be two surfaces from automatic and manual segmentation, respectively.

Experimental results
In Fig. 4, we show the long axis view of the automatically segmented ENDO-and EPI surfaces at frames 2, 4, 6, and 8 during ventricular systole. Fig. 4. Long axis view of segmented ENDO-and EPI contours at frames 2, 4, 6, and 8 during cardiac systole.
In Fig. 5, we compare the segmentation results of the EPI surface with and without incompressibility constraint. For a fair comparison, we used the same data adherence term for both cases. We observed that while the ENDO border was correctly detected even without the incompressibility constraint, the EPI contour leaked into other tissues, such as liver, that look like myocardium. This is because the LV myocardium/background contrast is lower than the LV blood pool/myocardium contrast, especially for the lateral and inferior sectors of myocardium. This low contrast obscures the exact location of the EPI boundary, making EPI segmentation more challenging than ENDO segmentation. When the incompressibility constraint was applied, however, the coupling of ENDO and EPI contours forced the evolution of the EPI contour to be consistent with that of the ENDO contour, thus preventing the leakage of the EPI contour. As explained in Section 1, another challenge in the segmentation of the EPI boundary is the presence of the LV/RV myocardium junctures. The intensity similarity between the LV and RV myocardium makes the EPI boundary ambiguous at the juncture. When the incompressibility constraint was not applied, the EPI contour leaked out to segment RV myocardium. When the incompressibility constraint was applied, however, the LV myocardium was successfully separated from the RV myocardium at the juncture. Fig. 6 compares the MV for an entire sequence from manual segmentation, automatic segmentation with incompressibility constraint, and automatic segmentation without incompressibility constraint. We took the mean value from three experts as the MV from manual segmentation. We observed that the variation of MV for manual segmentation is 92. . Also, we noticed that the MV obtained without incompressibility constraint is much larger than that from manual segmentation. This is because the EPI contour leaked into the background, leading to the over-estimation of the MV. Fig. 6. Comparison of MV for an example sequence from manual segmentation, automatic segmentation with incompressibility constraint, and automatic segmentation without incompressibility constraint.
In addition, we performed Bland-Altman analysis (Bland and Altman, 1986) to assess the agreement of the MV measurements from manual segmentation and automatic segmentation. The Bland-Altman analysis revealed small bias and good coherence between the MV measurements from manual segmentation and automatic segmentation with incompressibility constraint. As shown in Fig. 7, the bias is -0.2%, and 95% of the computer measurements of MV can be expected to differ from expert measurements by less than 6.4% below and 5.9% above the mean. In comparison, when the incompressibility constraint was not applied, the MV largely deviated from the MV from manual segmentation (bias = 18.8%, 95% confidence interval = [2.4%,35.2%]).  Table 1 shows that the ENDO boundaries were detected with sufficient accuracy even if the incompressibility constraint was not applied. This is because that the ENDO boundaries are relatively clear compared to the EPI boundaries. However, we observed from Table 2 that for EPI boundaries, when the incompressibility constraint was not applied, the automaticmanual MAD was 2.78mm larger than the manual-manual MAD, the automatic-manual HD was 3.82mm larger than the manual-manual HD, and the automatic-manual PTP&PFP were 18%-20% worse than the manual-manual PTP&PFP. This is because of the leakage problem of EPI boundaries without incompressibility constraint. When the incompressibility constraint was applied, however, the automatic-manual MAD decreased by 3.38mm, the automatic-manual HD decreased by 4.11mm, and the automatic-manual PTP&PFP decreased by 16%-19%. We can see from Tables 1 and 2 that the automatic algorithm produced results with comparable accuracy to a manual segmentation. Furthermore, we observed from Tables 1 and 2 that the variability of manual-manual segmentation was smaller for the ENDO surface than for the EPI surface. This is probably because the EPI boundaries are more ambiguous for observers to detect, which was also the reason why we have multiple observers instead of a single one.

Conclusion
In this chapter, we have presented a novel approach to segmenting the full myocardial volume from cardiac MR and ultrasound images. Motivated by the incompressibility property of myocardium during a cardiac cycle, we coupled the propagation of the ENDO and EPI surfaces according to image-derived information that maximized the piecewise homogeneities of a cardiac image, as well as the incompressibility constraint which constrains the variation of myocardial volume within 5%.
To validate our approach, we computed the MAD, HD, PTP, and PFP between the contours from automatic algorithm and manual segmentation. We observed that when the incompressibility constraint was not applied, the EPI boundaries leaked into other tissues. When incompressibility constraint was applied, however, the MAD, HD, PTP, and PFP were significantly improved.