An automatic technique for left ventricle segmentation from msct cardiac volumes

In this research, an automatic technique to segment the left ventricle from the heart information in multislice computed tomography images is proposed. A preprocessing stage is considered as a necessary preliminary task for diminishing the artifacts impact in the image analysis. With this idea, a similarity enhancement that combines a smoothed version of the original volume with a processed volume using mathematical morphology is used. This preprocessing approach is compared with respect to other strategies. After, a volume-of-interest is defined in order to isolate the cavity using two cropping planes detected with least squares support vector machines. Finally, the segmentations are obtained using both a region growing algorithm and a level sets algorithm. The robustness of each enhancement strategy is validated by performing the segmentation of images. This evaluation considered the Dice score, and both volume and surface errors. A clinical dataset from 12 patients is used in the inter- and intra subject evaluation. During intra-subject validation the proposed scheme achieves the best results, while a modified version of this scheme achieved the best performance during inter-subject validation.


Introduction
The left ventricle (LV) is considered the main heart cavity, therefore, the accurate description of its shape is important for assessing the cardiovascular function. Multi-dimensional medical imaging modalities have been developed to provide relevant information of the internal organs of humans. In 3D cardiac imaging, multislice computed tomography (MSCT) is an imaging modality useful for obtaining information about the morphology of the heart and vessels.
MSCT provides 3D images of the heart with sub-millimeter isotropic or anisotropic spatial resolution in temporal phases of each cardiac cycle [1]. MSCT images are strongly sensitive to changes in heart rate. The cardiac motion artifacts and slice overlaps are generated due to irregularities in ECG recording. The acquisition in inconsistent phases of the cardiac cycle leads to artifacts between stacks of axial  [2]. In addition, the noise (Poisson noise) in CT is directly related to the number of detected Xray photons. These detected photons could be modeled using a Poisson distribution [3].
Image segmentation is an important tool to analyze features related to anatomic tissue, and spatial distribution of functional regions. This tool enables simplification of 3D medical data without discarding important image features.
A big number of methodologies have been proposed for segmenting cardiac images and for extracting parameters describing the ventricular shape and function, thus increasing the frontiers of clinical diagnoses and research on cardiovascular diseases. In this sense, Uzunbas et al. [4], combine deformable models and a graph cut method for performing the semi-automatic left ventricle segmentation in 15 cardiac MRI sequences. The LV endocardium segmentation attained a mean and standard deviation of a Dice score of 0.82±0.06. An active contour model, used to detect the LV in cardiac CT images, was proposed in [5]. The model considers an external energy functional composes by an adaptive diffusion flow component and a localizing region intensity fitting. The accuracy of segmentation was evaluated using the Hausdorff distance and the volume overlap. A minimum distance of 7.21 was attained; meanwhile the maximum overlap was of 88.67%. A strategy that combines Hermite transform, active shape model and level sets to segment LV boundaries from tomographic cardiac studies is presented in [6]. A CT dual source scanner was used to acquire 28 studies from healthy subjects. These studies were also annotated. The Hausdorff distance, Dice score and ray feature error were used as quantitative metrics to evaluate the strategy. The results shown that the strategy accurately discriminated left ventricle.
This paper is an extension of the works reported in [7] and [8]. Essentially, the main contributions to the image enhancement scheme proposed are listed below: 1. Implementation of preprocessing strategies based on partials modifications of the proposed scheme. 2. Contrasting the proposed scheme with respect to less complex preprocessing techniques. 3. Application of least squares support vector machines (LSSVM) for: a) Volume-of-interest automatic definition. b) Automatic initialization of segmentation methods. 4. Evaluation of the preprocessing schemes considering inter and intra-subject variability.

Method
Ten preprocessing strategies are implemented and validated using clinical data. The validation is performed after segmentation. The preprocessed volumes are subsequently processed with two threedimensional segmentation methods. The first method is a region growing technique while the second corresponds to the level sets algorithm. The preprocessing strategy that provides the lower segmentation error should be selected as the best. Figure 1 shows an overview of the complete approach, including the preprocessing stage and segmentation algorithm.   Figure 1), can be consulted in [8].
2.1.2. Preprocessing stages. Two preprocessing stages are performed. First, the MSCT images are processed using nine enhancement strategies in order to minimize artifacts and Poisson noise. During the second stage, the volume-of-interest (VOI) is automatically established to exclude certain cardiac structures in the enhanced MSCT images sequences.

Enhancement strategies.
The proposed strategies denoted as S0 to S9 are described below.
2. Strategy S1: The MSCT sequences are processed using the scheme shown in Figure 1.
3. Strategy S2: The averaging filter, in Figure 1, is replaced by the unprocessed images. 4. Strategy S3: In this case, the averaging filter (in Figure 1) is replaced by a median filter. 5. Strategy S4: The averaging filter (in Figure 1) is replaced by an adaptive averaging filter. The smoothed volume is calculated by performing the following tasks: a) A gradient of the original volume is calculated. b) A slice located at the equator of the gradient volume is selected. The average gray level of voxels belonging to the inner region of LV is calculated. This average represents the threshold value that controls the averaging filter. c) The following rule is then applied: the gray level of the current voxel f(i, j, k) is replaced by the average gray level calculated in a 3D neighborhood of the current voxel if and only if f(i, j, k) is greater than the predefined threshold. 6. Strategy S5: The averaging filter (in Figure 1) is replaced by an anisotropic diffusion filter. 7. Strategy S6, S7 and S8: These strategies are obtained by replacing the similarity enhancement block of the diagram in Figure 1, by an anisotropic diffusion filter (S6), a Gaussian filter (S7) and a multi-scale Gaussian filter (S8), respectively. 8. Strategy S9: In the diagram of Figure 1, the averaging filter is replaced by the original unprocessed MSCT sequences while the top-hat filter is replaced by a gradient magnitude filter.

Volume-of-interest definition.
LV segmentation based on region growing techniques usually fails due to the fact that gray levels of different cardiac structures have almost similar values. To address this problem a VOI is defined. In this sense, two cropping planes are built and then used to exclude certain anatomical structures. The construction of each cropping plane requires at least two points located in the MSCT volume. The left atrium-LV joint (p1) and apex (p2) are the points for a cropping plane, whereas the aortic valve-LV joint (p3) and p2 are the respective points to the second plane. The planes constructed with points described above are henceforth referred as mitral plane and aortic plane, respectively. The automatic detection of such points is performed by applying an automatic learning approach. This approach considers a wavelet-based scale reduction followed by the application of a least square support vector machine (LSSVM) and it is described in [9].
The VOI is constructed as follows: a) The coordinates of p1, p2 and p3 are mapped to the enhanced volumes in full resolution size. b) The pairs of points p1-p2 and p3-p2 are used to construct the mitral and aortic planes. c) Plane orientation is determined by the normal vector obtained using these points.
The location of the planes is determined by the p1 and p3, respectively. d) Two linear discriminant functions divide the MSCT volume information using the mitral and aortic planes as hyperplanes decision surfaces, respectively.   6. Region growing algorithm. The complete segmentation algorithm is described in [7]. This regionbased method uses the hybrid-linkage region growing algorithm in order to group voxels into 3D regions. The 3D hybrid-linkage technique starts also with a seed that lies inside the region of interest and spreads to the p-connected voxels that have similar properties. This technique assigns a property vector to each voxel where the property vector depends on the neighborhood of the voxel.
2.1.7. Sparse field level sets algorithm. The level sets method allows the iterative deformation of certain geometric structures called snakes. The sparse field algorithm based on uniformity of regions is chosen to perform the LV segmentation to exploit their ability for both contraction and expansion. The ItkSnap implementation of the sparse field LS algorithm is used to perform the segmentation [10].
A voxel in an optimal location is needed for initializing the segmentation algorithms. This voxel represents the seed in the RG algorithm while in the LS algorithm represents the center of the initial shape. An additional LSSVM is trained and used for detecting the necessary voxel to initialize the algorithms. A 3-D iso-sphere is used as an initial shape in the LS algorithm. The sphere radius is fixed at 20 mm such that the iso-surface is fully contained within the LV for all enhanced volumes.
The threshold required for generation of the features image matches with the average gray levels contained in the interception of the image preprocessed and the aforementioned iso-sphere. By default, the ItkSnap considers a sigmoid function as a transfer function for performing the thresholding process.
The stopping criterion used by default in the ItkSnap is the number of iterations. It is known that there is a relationship between the spatial dimension of an image and the number of iterations required for attaining the segmentation of a particular structure. In consequence, we have chosen 512 as the maximum number of iterations, which is the spatial size of the MSCT slices [11].

Data source
In total, 44 multi-slice computerized tomography volumes are used to evaluate the proposed image enhancement scheme for 3-D segmentation of the left ventricle. The data used in this research are sequences 4-D (3-D + time) of cardiac images that have been acquired using two scanners. Each dataset consisted of 20 volumes representing anatomical information, for a complete cardiac cycle. The number of patients used in the experiments is twelve. There is a database for each patient, which will be denoted henceforth as DB1 to DB12. Each volume could contain between 148 and 326 slices. The slice thickness varies from 0.400 mm to 0.625 mm. In all volumes, the slices have a resolution of 512×512 pixels, but the pixel size varies from 0.273 mm to 0.488 mm. The voxel value is represented by 12 bits.

Validation process
The idea is to develop a methodology useful for evaluating the inter and intra-subject variability of the complete approach. In this sense, the images segmented using the current approach are compared with respect to segmentation manually traced by a cardiologist in the same images. The error metric proposed by Suzuki and Dice [7], are reformulated in 3D domain and then considered in the methodology as segmentation errors. The MSCT sequences are enhanced using the nine strategies. From these results, multiple left ventricle shapes are estimated using the segmentation algorithms. For each dataset, the segmentation errors are calculated. Finally, the best enhancement strategy is chosen from the segmentation that minimized the Suzuki and Dice metrics.
In the inter-subject validation eleven datasets are considered (DB1 to DB11). In order to obtain the ground truth, the cardiologist chooses the image volumes at end-diastolic and at end-systolic in the MSCT sequences. Then, the specialist guides a manual process to segment the 22 volumes chosen. These segmentations are the end-diastolic and end-systolic ground truths used as references to indirectly quantify the quality of the image enhancement scheme.
In the intra-subject validation a cardiologist segments, manually, the 20 volumes of the remaining database, DB12. These shapes are used as references during this validation.  Figure 2 shows the coronal view of a raw image and its smoothed versions after applying each of the filtering strategies. Figure 2(a) shows the unfiltered image where the stairstep artifact is visible. The Figure 2 also shows how each of the filtering strategies reduces, in some degree, the effect of the mentioned artifact.

Volume-of-interest results
The process used to define the VOI generates images in which the LV is isolated from other structures. Figure 3 shows the orthogonal views of the enhanced with strategy S1.  Table 1 presents the averages and standard deviations (μ ± σ) of errors for the metrics, Dice score (Ds), volume error (Ve) and surface error (Se), that quantify the performance of the filtering strategies after segmentation of databases used for inter-subject validation. The values in this table are for both segmentation algorithms. The best Dice score (Ds) was obtained by S9, followed closely by S1. S9 also obtained the lower volume and surface errors for both segmentation algorithms. The Dice score is highly correlated with the volume and surface errors. A closer look to this parameter suggests that the segmentations obtained after applying filtering strategies (from S1 to S9) produce better results than those obtained by the segmentation of volumes without any filtering (strategy S0). S0 is clearly suboptimal. Figure 4 shows the segmentations generated using a level sets (LS) algorithm, after applying S9 during inter-subject validation. Segmentations generated by LS were generally of better quality than those generated by the RG algorithm.  Table 2 presents a summary of the averages and standard deviations (μ±σ) of errors calculated for the filtering strategies after performing the segmentation of database DB12. Results show that the best performance is obtained using strategy S1, followed closely by S9. These two strategies generated both the best of Ds and the lower Ve and Se, for both algorithms. Considering the Ds, S8 was the strategy with lower performance followed by S4. The rest of strategies exhibit a behavior close to each other. Considering that S1 is the best strategy, Figure 5, shows the LVsegmentations generated by LS after applying S1. Figure 5. LV segmentations using S1 and LS during intra-subject validation.

Results discussion
The values of Dice coefficient, reported for S1 and S9, allow to affirm that they produce the threedimensional morphology of the LV with an excellent precision, when compared with both the manual segmentations and with the results of other researchers that report this metric.

Conclusions
The results obtained for strategy S0 compared with other filtering strategies suggest that application of an appropriate filtering strategy is necessary. The group of strategies that includes S6, S7 and S8 which are based on the application of a single filter technique, failed when compared with strategies S1 and S9 which are based on similarity enhancement. Note that, during intra-subject validation the strategy that attained the best results was S1; while S9 achieved the best performance during inter-subject validation. S9 was the most efficient strategy, as it has the lowest average computational cost. This was due to the substitution of the averaging filter by a gradient magnitude filter which is less complex than the filter used in S1 which is based on mathematical morphology.
In general, the 3D shapes representing the LV obtained by segmentation using the LS algorithm were smoother than the shapes obtained using the RG algorithm. The segmentation based on the level sets algorithm attained higher values for the Dice score with respect to the region growing algorithm. The average Dice score for the levels set algorithm was almost three points higher than the average Dice score for the región growing algorithm. In contrast, the region growing algorithm has a lower computational cost with respect to the level sets algorithm. For instance, the levels set algorithm required 310.2 seconds for segmenting a LV shape while the region growing algorithm required only 30 seconds.