Abstract
This paper presents a new unsupervised image segmentation method based on particle swarm optimization and scaled active contours with shape prior. The proposed method uses particle swarm optimization over a polar coordinate system to perform the segmentation task, increasing the searching capability on medical images with respect to different interactive segmentation techniques. This method is used to segment the human heart and ventricular areas from datasets of computed tomography and magnetic resonance images, where the shape prior is acquired by cardiologists, and it is utilized as the initial active contour. Moreover, to assess the performance of the cardiac medical image segmentations obtained by the proposed method and by the interactive techniques regarding the regions delineated by experts, a set of validation metrics has been adopted. The experimental results are promising and suggest that the proposed method is capable of segmenting human heart and ventricular areas accurately, which can significantly help cardiologists in clinical decision support.
1. Introduction
In clinical practice, magnetic resonance imaging (MRI) and computed tomography (CT) scanning are effective and widely used methods for the diagnosis and monitoring of cardiac disease. The process performed by a cardiologist on medical images consists of a visual examination followed by a manual delineation of the human organ, which can be subjective, susceptible to errors, and time consuming. Accordingly, the application of computational methods in order to obtain an efficient and accurate image segmentation for clinical decision support plays an essential role.
Automatic segmentation of human organs is an important and challenging task in medical image analysis. In recent years, several approaches have been reported for this purpose such as rule optimization with region growing in pelvic injuries [1], suppressed fuzzy c-means in brain magnetic resonance images [2], adaptive local multiatlas in human heart [3], graph cut in multiple human organs [4, 5] active contour models (ACM) in lungs from magnetic resonance images of the torso [6], intravascular ultrasound images [7], and human prostate [8].
ACM is an energy-minimizing spline curve consisting of a set of discrete control points known as snaxels. This spline is attracted towards features such as edges or lines through the evaluation of internal and external forces. ACM is prone to stagnate in local minima and is also highly sensitive to initialization because it requires to be close to the target object, otherwise failure of convergence will occur.
Since ACM was proposed by [9], many variations have been developed to improve these shortcomings through the introduction of some prior knowledge such as active shape models [10], ACM based on level set method [11], shape prior in left ventricle [12] and cerebellum [13], registered active shape model (RASM) in traumatic pelvic CT images [14], and ACM with particle swarm optimization (PSO) for womb fibroma in supersonic images [15], where static searching windows (size of 10 × 10 or 15 × 15 pixels) were dynamically generated depending on the initial position of the interactive control points. The performance of PSO approach is robust in local minima problem, and according to the tests, it provides an accurate medical image segmentation within an appropriate executing time.
PSO is a stochastic and population-based method inspired by the cognitive and social behavior of bird flocking to solve optimization problems in continuous spaces [16, 17]. This computational intelligence technique consists of a set of potential solutions known as swarm, where each potential solution is referred to as particle. In the PSO strategy, all the particles are guided by the best particle of the swarm, and each particle keeps track of its best solution found through iterations. Since being not computationally expensive and highly efficient to solve optimization problems, PSO has been used in a wide range of medical applications such as hepatitis disease diagnosis [18], tumor classification [19], and 3D medical registration [20].
Although the aforementioned algorithms provide satisfactory results regarding accuracy segmentation and noise sensitivity, more efforts are necessary to develop clinical decision support methods. In this paper, we introduce a new unsupervised image segmentation method based on particle swarm optimization and scaled active contours with shape prior. The proposed method uses PSO over a polar coordinate system to perform the segmentation task by increasing the exploration and exploitation capabilities with respect to the traditional ACM. In addition, this method utilizes the alignment process proposed in [21] to construct a target object template, which is used to determine the initial positioning of the polar coordinate system. Then, the template is scaled to a different size in order to generate potential solutions and assuming that the target object is confined within them. This proposed method addresses the segmentation problem of human heart and ventricular areas from CT and MR images, respectively. Furthermore, to evaluate the segmentation results with respect to regions outlined by experts and by different computational methods, a set of validation metrics is presented.
The remainder of this paper is organized as follows. In Section 2, we introduce the fundamentals of particle swarm optimization and active contour model with shape prior, along with the proposed segmentation method. The experimental results are discussed in Section 3, and conclusions are given in Section 4.
2. Materials and Methods
2.1. Overview of Particle Swarm Optimization (PSO)
PSO is an artificial intelligence algorithm proposed by [16] and modified by [17] to solve optimization problems. PSO consists of a swarm of potential solutions called particles. Each particle is viewed as a point in an -dimensional space . During each iteration, every particle moves through hyperspace to a new position according to the following velocity equation: where is the current particle in the time step and represents its velocity, is the inertia weight, where is a uniform distribution, represents the learning factor, is the best solution found by the current particle, and is the best solution found by the best particle of the whole swarm. On the other hand, assuming that the new velocity of the current particle has been updated, (2) is used to compute its new position within the search space: According to the above description, the PSO algorithm can be implemented through the following procedure. (1)Establish the swarm size and randomly initialize the position and velocity of each particle. (2)Evaluate each particle in the fitness function in order to update its , if the new fitness is better. (3)Find the best particle in the whole swarm and update , if the fitness value found is better. (4)Stop if the convergence criterion is satisfied (e.g., stability or number of iterations). (5)Update velocity and position of all the particles using (1) and (2), respectively, then repeat steps (2)–(5).
2.2. Active Contour Model with Shape Prior
The conventional active contour model is a parametric curve that is driven by internal and external forces to minimize its energy function [9]. In order to incorporate a shape prior constraint within the active contour, the Chan and Vese method [11] is used, and it is defined by using the following equation: where is the total energy composed of the energies , and and the weighting factors of the shape energy , and . The first energy is the active contour, which can be represented as where is the image domain, is the Heaviside function, is the image intensity, is the gradient operator, and are the weighting parameters of the length and area energies of the contour, is a signed distance function, and and are the mean intensity of the object and background, and they are given by the following equations: is the shape energy defined by the difference between the active curve and the shape template. This energy will be minimized by optimizing the transformation parameters, and it is expressed as follows: where represents the deformed template and is defined in a transformation matrix consisting of translation , scaling , and rotation parameters as follows: where is a scaling factor, is the rotation angle parameter, and are translation parameters in the horizontal and vertical axes. Finally, the third energy term is the image-based force, which is the difference energy, and it is computed as
The three energies of the active contour are iteratively evaluated, and the contour evolution terminates when the difference between the previous and current segmented area becomes stable.
2.3. Proposed Image Segmentation Method
The proposed image segmentation method consists of three main steps as shown in Figure 1, and they are described below.
2.3.1. Shape Representation and Construction of the Aligned Template
In order to generate a template of the target object, a set of selected reference images is aligned, which leads to differences in position, direction, and scale. The aim of this step is to obtain a shape template through the alignment of a set of manually segmented images, which contain the human heart and ventricular areas. In Figure 2, to analyze the alignment procedure, a training set consisting of 8 human hearts is presented.
(a)
(b)
By using the technique developed in [21], we compute the shape alignment by estimating the parameters through (10), where is the and translation, is the scale parameter, and is the rotation angle:
The product of the translation matrix , scaling matrix , and the in-plane rotation matrix maps the coordinates to coordinates . Following the alignment process, the gradient descent method is used to minimize the following energy function: where is the transformed image based on the shape parameters, and is the image domain. The gradient of in (11) is iteratively evaluated until convergence. The alignment results presented in Figure 3(a) are slightly different from the training binary shapes, and in Figure 3(b), the resulting maximum contour after alignment presents a significant variation regarding the maximum contour shown in Figure 2(b).
(a)
(b)
After adjusting each shape parameter, the final aligned template is obtained by superimposing all transformed images, and then, the maximum shape boundary is acquired.
2.3.2. Multiswarm Initialization and Numerical Optimization
The initialization stage is performed on the resulting distance map, where the origin of the polar coordinate system is automatically determined by using the maximum mutual information between the template and the current medical image. The coordinate system divides the target object through , where denotes the number of degrees of each constrained polar section . In our approach, parameters play an important role, which are described as follows number of scaled contours; this parameter has to be considered assuming that the target object is confined within the region of the initial contours. Number of control points determines the number of polar sections in which the target object is divided. Number of iterations is used to obtain proper segmentation results through the evaluation of the fitness function. Inertia weight controls the exploration and exploitation abilities of the swarm by weighing the previous velocity value on the new velocity and the learning factor parameter to scale the step size of the search. The numerical optimization is performed after the scaled contours are produced and the control points (particles) are generated for each constrained polar section , in which one edge sectional solution and one swarm of particles must exist. PSO strategy is applied in each polar section separately in order to be placed on its corresponding edge sectional solution. All the particles are evaluated according to the fitness function derived from the distance map, where through iterations the best particle of each swarm is updated only if a best particle is found in its constrained search space. When the optimization process for each swarm is finished, the resulting segmented object is acquired by connecting the particle of each swarm to each other.
The proposed method presents the following advantages: (1) the initial contours are automatically initialized according to the shape template obtained from the alignment process; (2) the number of discrete points can be adjusted directly by modifying the number of polar sections through the parameter. Due to these advantages, the proposed method can be extended to work with sequential CT and MR images by just applying the maximum mutual information to reproduce the coordinate system on the set of images.
The procedure of the proposed image segmentation method is described as follows. (1)Align predefined shapes according to [21] and obtain template after alignment. (2)Perform maximum mutual information to initialize coordinates of the polar coordinate system. (3)Initialize degrees and number of snakes. (4)Initialize the PSO parameters: number of iterations, inertia weight, and learning factor. (5)Generate one swarm of particles for each polar section assigning the current control points as particles. (6)For each swarm , initialize velocities and assign the initial and . (a)Apply restriction of the search space to ignore improper particles.(b)Evaluate each particle in fitness function derived from the distance map.(c)Update and if better particles are found.(d)Apply (1) and (2), respectively.(e)Stop if the convergence criterion is satisfied (e.g., stability or number of iterations), otherwise go to step (a).
2.4. Evaluation Metrics
In order to assess the performance of the proposed method, the maximum cardinality similarity metric, Hausdorff distance, Jaccard, Dice, and correlation indices have been adopted to compare the segmentation results between the regions outlined by experts and the regions obtained by computational methods:
The Jaccard and Dice indices are computed using (12) and (13), and they are used as measures of similarity between reference () and automatic () segmented objects [13]. If regions and are entirely overlapping, the obtained result is 1, and it is 0 when these two regions are completely different. In addition, the correlation index (14) measures the linear relationship between the reference and automatic segmentations, and it is defined in the range [−1,1]. A correlation of 1 means perfect positive linear relationship, and −1 means negative linear relationship.
On the other hand, the Hausdorff distance is a metric for shape matching in medical image segmentation, which measures the similarity between two superimposed sets via (15), where and are the edge pixels in sets and , respectively, and is the Euclidean distance:
The maximum cardinality similarity metric is used in template-matching applications with a low computational complexity [22]. It is applied to the edge pixels of segmented objects using (16), where TP is the number of true positives, FN is the false negatives, and is the total number of pixels in the image:
3. Experimental Results
In this section, we evaluate the performance of the proposed image segmentation method for segmenting the human heart and ventricular areas on CT and MR images, respectively. The computational implementations were performed using the GNU Compiler Collection (C++) version 4.4.5 running on Debian GNU/Linux 6.0, Intel Core i3 with 2.13 Ghz and 4 Gb of memory.
In Figure 4, the segmentation results on a subset of CT images containing the human heart are presented. The whole dataset consists of 144 CT images of size 512 × 512 pixels obtained from different patients. In Figure 4(a), the human heart outlined by cardiologists is presented. Figure 4(b) illustrates the segmentation results obtained through the interactive Graph Cut method [5], which were computed with an average executing time of 0.185 s per image. In Graph Cut method, the experts defined the human heart area and the background seeds in an interactive way. In Figure 4(c) the segmentation results by using the traditional ACM are presented, in which the noise sensitivity and fitting problem are shown. The ACM parameters were determined according to the multiple objects test presented by [15] as 45 control points, , , and , obtaining in our test an average executing time of 0.157 s per image. Figure 4(d) shows the human heart segmentations obtained with the interactive Tseng method. The parameters of this method were determined according to [15], and they were set as 45 control points, window size as 30 × 30 pixels, and 9 particles for each swarm, given an average executing time of 0.176 s per image. Finally, in Figure 4(e), the segmented images by using the proposed method reveal an appropriate human heart segmentation avoiding the local minima problem. In this simulation, the parameters were set as number of scaled contours = 9, number of snaxels = 45, iterations = 10, inertia weight = 0.5, and learning factor = 0.9, obtaining an average executing time of 0.198 s per image. In our approach, the learning factor and inertia weight were tuned experimentally taking into account the following two considerations: firstly, the number of different potential solutions generated through the iterations and secondly, by considering the number of improper solutions in order to perform local exploitation instead exploration [23, 24].
(a)
(b)
(c)
(d)
(e)
From the dataset of CT images described above, in Table 1, a subset of the segmentation results is presented. These results were obtained by comparing the regions delineated by experts with the segmentation results of the Graph Cut, ACM, Tseng method, and our proposed method. Based on the similarity analysis, the human heart segmentation results suggest that the proposed method can lead to more efficiency in the presence of concavities and noise than the comparative interactive methods.
Figure 5 presents the segmentation results on a subset of MR images containing the human left ventricle. These images have been extracted from a previously segmented dataset with 23 MR images of size 512 × 512 pixels. In Figure 5(a), the human left ventricle delineated by cardiologists in different test images is presented. Figure 5(b) shows the segmentation results obtained in an average executing time of 0.166 s by using the interactive Graph Cut method, in which experts defined the human left ventricle area and the background image interactively. These segmentation results cannot be adjusted to the correct left ventricle boundary, which is improved by the next three methods described below. Figure 5(c) illustrates the segmentation results through the traditional ACM, where the local minima problem is clearly shown. The ACM parameters were adjusted according to the medical tests presented in [15], and they were set as 35 control points, , , and , obtaining an average executing time of 0.108 s per image. In Figure 5(d), the human heart segmentation results by applying the interactive Tseng method are illustrated. The parameters of this simulation were chosen according to supersonic image test presented by [15] as 35 control points, window size as 30 × 30 pixels, and 15 particles for each swarm, given an average executing time of 0.209 s per image. Moreover, Figure 5(e) presents the segmentation results acquired by the proposed method. The segmented images show in a suitable way the boundary of the human left ventricle avoiding the noise sensitivity problem. The parameters of this experiment were set as number of scaled contours = 7, number of snaxels = 35, iterations = 10, inertia weight = 0.4, and learning factor = 0.7, obtaining an average executing time of 0.169 s per image.
(a)
(b)
(c)
(d)
(e)
According to the left ventricle dataset of MR images previously described, in Table 2, the average of the segmentation results obtained by Graph Cut, ACM, Tseng method, and our proposed method is compared with the regions outlined by experts. This similarity analysis shows that the proposed methods perform the human left ventricle segmentation accurately with respect to the other computational methods.
In Figure 6, the segmentation results on a subset of coronal MR images containing ventricular areas are presented. These images have been extracted from a dataset of 19 images with size 256 × 256 pixels. In Figure 6(a), the manual delineations made by experts on the ventricular area are illustrated. Figure 6(b) presents the segmentation results obtained with the interactive Graph Cut, which were acquired with an average executing time of 0.192 s per image. This method fails to fit the ventricular area boundary because of the presence of noise. In Figure 6(c), the traditional ACM is applied, where the segmented ventricular areas cannot fit the true boundary accurately. The ACM parameters were set according to the medical image test introduced by [15] as 25 control points, , , and , obtaining in our tests an average executing time of 0.095s per image. Consequently, in Figure 6(d), the ventricular area segmentation results through the interactive Tseng method are shown. The parameters of this simulation were tuned according to the medical image test presented by [15] as 25 control points, window size as 30 × 30 pixels, and 15 particles for each swarm, given an average executing time of 0.143 s per image. Figure 6(e) shows the segmentation results obtained by the proposed method, where the segmented images fit the true ventricular area boundary more accurately than the previous computational methods. The parameters of this simulation were set as number of scaled contours = 7, number of snaxels = 25, iterations = 10, inertia weight = 0.4, and learning factor = 0.7, achieving an average executing time of 0.154 s per image.
(a)
(b)
(c)
(d)
(e)
In Table 3, the average of the previously segmented ventricular areas through the Graph Cut, ACM, Tseng method, and our proposed method is compared to the regions outlined by experts. The performed similarity analysis suggests that the proposed method is more robust in ventricular area segmentation than the other computational methods according to the Jaccard and Dice indices with 85% and 92%, respectively.
Segmentation results were compared with the regions outlined by cardiologists on CT and MR images in the aforementioned datasets. As shown in similarity analysis, the proposed method is able to detect human cardiac organs with a high accuracy and effectiveness. Moreover, it can help cardiologists to better analyze the medical images and increase their monitoring abilities.
4. Conclusions
We have presented a new unsupervised image segmentation method based on particle swarm optimization and scaled active contours with shape prior. This method generates different scaled contours according to the shape template obtained from an alignment process. Then, particle swarm optimization is used to perform the segmentation task over constrained polar sections. This novel approach was applied to segment the human heart and ventricular areas from CT and MR images. In order to evaluate the segmentation results acquired through the proposed method, correlation and Jaccard and Dice indices were used. According to the experimental results, the proposed method proves to be capable of segmenting human heart and ventricular areas accurately compared to the segmentations acquired by Grap Cut, traditional ACM, and interactive Tseng method. In addition, the experimental results have also shown that the exploitation capability of the proposed method is highly appropriate for cardiac medical image applications.
Acknowledgments
This research has been supported by the National Council of Science and Technology of México (CONACYT) under Grant no. 241224-218157. The authors thank the Auckland MRI Research Group, University of Auckland, and Cardiac Atlas Website for the valuable collaboration providing us the sources of magnetic resonance imaging. The authors wish to thank the Cardiology Department of the Mexican Social Security Institute, UMAE T1, León, for their clinical advice.