- Email: [email protected]
A comparison between two robust techniques for segmentation of blood vessels
Computers in Biology and Medicine 38 (2008) 931 – 940 www.intl.elsevierhealth.com/journals/cobm
A comparison between two robust techniques for segmentation of blood vessels Roberto Rodríguez a , Patricio J. Castillo a , Valia Guerra b , Juan Humberto Sossa Azuela c,∗,1 , Ana G. Suáreza a , Ebroul Izquierdo d a Digital Signal Processing Group, Institute of Cybernetics, Mathematics & Physics (ICIMAF), CP 10 400, Havana, Cuba b Numerical Methods Group, Institute of Cybernetics, Mathematics & Physics (ICIMAF), CP 10 400, Havana, Cuba c Center for Computing Research (CIC), National Polytechnic Institute (IPN), Mexico d Department of Electronic Engineering, Queen Mary, University of London, UK
Received 22 August 2007; accepted 18 June 2008
Abstract Image segmentation plays an important role in image analysis. According to several authors, segmentation terminates when the observer’s goal is satisfied. For this reason, a unique method that can be applied to all possible cases does not yet exist. In this paper, we have carried out a comparison between two current segmentation techniques, namely the mean shift method, for which we propose a new algorithm, and the so-called spectral method. In this investigation the important information to be extracted from an image is the number of blood vessels (BV) present in the image. The results obtained by both strategies were compared with the results provided by manual segmentation. We have found that using the mean shift segmentation an error less than 20% for false positives (FP) and 0% for false negatives (FN) was observed, while for the spectral method more than 45% for FP and 0% for FN were obtained. We discuss the advantages and disadvantages of both methods. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Image segmentation; Mean shift; Spectral methods
1. Introduction Nowadays medical imaging provides major aid in many branches of medicine; it enables and facilitates the capture, transmission and analysis of medical images as well as providing assistance towards medical diagnoses. Medical imaging is still on the rise with new imaging modalities being added and continuous improvements of devices’ capabilities. In particular, segmentation and contour extraction are important steps in many medical imaging systems. Segmented images are now routinely used in a multitude of different applications, such as diagnosis, treatment planning, robotics, pathology, study of ∗ Corresponding author. Tel.: +52 55 57 29 60 00x56505; fax: +52 55 57 29 60 00x56607. E-mail address: [email protected] (J.H.S. Azuela). 1 Centro de Investigación en Computación-IPN, Av. Juan de Dios Bátiz s/n, Esquina con Miguel Othón de Mendizábal, Colonia Nueva Industrial Vallejo, C. P. 07700, México, D. F. MÉXICO.
0010-4825/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2008.06.002
anatomical structure and computer-integrated surgery. Image segmentation remains, however, a difficult task due to, in the one hand, the variability of the shape of the objects and, in the other hand, the variation in image quality. Medical images are often corrupted by noise and sampling artifacts, which can cause considerable difficulties when applying rigid methods. Nevertheless, in spite of the most complex algorithms developed until now, segmentation continues to be very dependent on the application, and a single method that can solve all the problems that can be encountered does not exist. Pathological anatomy is a speciality where the use of different digital image processing technique allows improving the accuracy of the diagnosis of many diseases. One of those pathologies is cancer, which is known to be one of the first causes of death in many countries. For this reason, the study of this disease by efficient methods of image processing is of supreme importance. Many segmentation methods have been proposed for medical-image data [1–6]. Unfortunately, segmentation using
932
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
traditional low-level image processing techniques, such as thresholding, region growing and other classical operations requires a considerable amount of interactive guidance in order to attain satisfactory results. Automating these model-free approaches is difficult because of complexity, shadows, and variability within and across individual objects. Furthermore, noise and other image artifacts can cause incorrect regions or boundary discontinuities in objects recovered from these methods. Many segmentation methods have been proposed for the study of blood vessels (BV) [7–11]. The main goal of this paper is to carry out a comparison between two segmentation methods: the mean shift method and the spectral method. In the case of the mean shift method we proposed a new segmentation algorithm. Our interest in comparing these two techniques is due to the good results that both have offered in other segmentation fields. A goal of this research is to determine which of the methods is the most suitable, robust and faster for segmenting images containing BV. In this application the important information to be extracted from images is just the number of BV. Mean shift is a nonparametric and versatile tool for feature analysis and can provide reliable solutions for many vision tasks. Refer for example to [12,13]. The mean shift was proposed in 1975 by Fukunaga and Hostetler [14] and largely forgotten until Cheng’s paper [15] rekindled interest in it. The spectral method was the other technique to be compared. It is based on the calculation of the eigenvectors and eigenvalues of an obtained matrix from the affinities of the pixels within the image. In particular, a spectral method was used which is called of normalized cut [16,17]. This method measures both, the total dissimilarity between the different groups as well as the total similarity within groups. The main idea is to solve the eigenvalue problem for a small subset of pixels (sample) and then this solution is extrapolated to the full set of pixels. The theoretical aspects of the used methods are exposed later on. The remaining of the paper is organized as follows. In Section 2, the characteristics of the images under study are described. In Section 3, the details of the most significant theoretical aspects of each of the two methods are presented. In Section 4, the steps in the experimentation are explained. Here, we carry out an analysis and discussion of both methods. In addition, the
method of evaluation is also given. Finally, in Section 5 the most important conclusions of this research are given. 2. Method and features of the studied images As biomedical images, we decided to use biopsies, which represent histological cuts in malignant tumors. These biopsies were included in paraffin by using the immunohistoquimic technique with the complex method of avidina biotina. Finally, monoclonal CD34 was contrasted with green methyl in order to accentuate new formation of BV. These biopsies were obtained from soft parts of human bodies. This analysis was carried out for more than 80 patients. In Fig. 1 typical images are given. These images were captured via a Morphometrical Analysis by Digital Image Processing (MADIP) system with a resolution of 512 × 512 × 8 bit/pixels. MADIP was a system created to help physicians in their research work using images [18]. MADIP is a software for the morphometrical characterization of macro and micro lesions in biomedical images. In general, it may be used in any branch of science and technology, where not only quantitative techniques of high precision are required, but it is also necessary to do a big number of measurements with high accuracy and rapidly. In Fig. 2 one can be observe a horizontal profile through the center of a vessel; that is, a plot of the intensities of pixels along a single row. By carefully observing this kind of images, we have found several characteristics of these images, which are common to typical images that we encounter in the tissues of biopsies: 1. The intensity was slightly darker within the blood vessel than in the local surrounding background. It was emphasized that this observation held only within the local surroundings. 2. High local variation of intensity was observed both within the blood vessel and the background. However, the local variation of intensity was higher within the blood vessel than in background regions (see Fig. 2(b)). 3. The variability of BV both, in size and shape can be observed (see Fig. 1). It is important to point out that due to the adopted acquisition protocol the images were corrupted with a lot of noise.
Fig. 1. These images represent the blood vessels in malignant tumors. Some blood vessels are marked with arrows.
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
933
Fig. 2. A profile of intensity through the center of a vessel. Profile is indicated by a line in Fig. 2(a).
The quantity n x /n(h d cd ) is the kernel density estimate fˆU (x) (where U means the uniform kernel) computed with the hypersphere Sh (x), and thus expression (2) can be written as
3. Theoretical considerations 3.1. About the mean shift method The iterative procedure to compute the mean shift is introduced as a normalized density estimate of the gradient [13]. By employing a differentiable kernel, an estimate of the density gradient can be defined as the gradient of the kernel density estimate; that is ∇ˆ f (x) = ∇ fˆ(x) =
n 1
nhd
∇K
i=1
x − xi h
1/2cd−1 (d + 2)(1 − x2 ) if x < 1 0 otherwise
K E (x) =
The density gradient estimate becomes ∇ˆ f E (x) =
1 d +2 · n(h d cd ) h2
(xi − x)
xi ∈Sh (x)
⎛ nx d +2 ⎝ 1 = · n(h d cd ) h2 nx
⎞ (xi − x)⎠
(2)
xi ∈Sh (x)
where the region Sh (x) is a hypersphere of radius h having volume h d cd , centered at x, and containing n x data points, that is, the uniform kernel. In addition, in this case d = 3, due to vector x is a three-dimensional vector; it is of dimension two for the spatial domain and one for the range domain (gray levels). The last term in expression (2) is called the sample mean shift, Mh,U (x) =
1 nx
(xi − x) =
xi ∈Sh (x)
1 nx
xi ∈Sh (x)
(4)
which yields Mh,U (x) =
h 2 ∇ˆ f E (x) d + 2 fˆU (x)
(5)
(1)
Conditions on the kernel K (·) and the window radio h are derived in [14] to guarantee asymptotic unbiasedness, meansquare consistency, and uniform consistency of the estimate in expression (1). For example, for Epanechikov kernel
d +2 Mh,U (x) ∇ˆ f E (x) = fˆU (x) · h2
xi − x
(3)
Eq. (5) shows that an estimate of the normalized gradient can be obtained by computing the sample mean shift in a uniform kernel centered on x. In addition, the mean shift has the direction of the gradient of the density estimate at x when this estimate is obtained with the Epanechnikov kernel. Since the mean shift vector always points towards the direction of the maximum increase in the density, it can define a path leading to a local density maximum; that is, to a mode of the density. In addition, expression (5) shows that the mean is shifted towards the region in which the majority of the points reside. Since the mean shift is proportional to the local gradient estimate, it can define a path leading to a stationary point of the estimated density, where these stationary points are the modes. Moreover, as it was pointed out, the normalized gradient in expression (5) introduces a desirable adaptive behavior, since the mean shift step is large for low-density regions corresponding to valleys, and decreases as x approaches a mode [13]. In [13] it was proven that the mean shift procedure, obtained by successive: • computing the mean shift vector Mh (x), • translating the window Sh (x) by Mh (x), guarantees convergence. It is well known that a digital image can be represented as a two-dimensional array of p-dimensional vectors (pixels), where p = 1 in the gray level case, three for color images, and p > 3 in the multispectral case.
934
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
As was pointed in [13] that when the location and range vectors are concatenated in the joint spatial-range domain of dimension d= p+2, their different nature has to be compensated by proper normalization of parameters h s and h r . Thus, the multivariable kernel is defined as the product of two radially symmetric kernels and the Euclidean metric allows a single bandwidth for each domain, that is
r 2 x s 2 C k y (6) K h s ,h r (x) = 2 p k h h hs hr s r where x s is the spatial part, y r is the range part of a feature vector, k(x) the common profile used in both domains, h s and h r , the employed kernel bandwidths, and C the corresponding normalization constant. Basically, the algorithm that we used in this work is performed into two steps, a filtering step and a segmentation step. The filtering algorithm comprises the following steps [13]: Let {X i }i and {Z i }i , i = 1, . . . , n be the input and filtered images in the joint spatial-range domain. For each pixel p ∈ X i , p = (x , y , z ) ∈ R3 , where (x , y ) ∈ R2 and z ∈ [0, 2 − 1], being the quantity of bits/pixel in the image. Then, for eachi = 1, . . . , n . 1. Initialize j = 1 and y i ,1 = x i . 2. Compute through the mean shift y i ,j +1 (see expression (3)), the mode where the pixel converges; that is, the calculation of the mean shift is carried out until convergence, y = y i ,c . 3. Store at Z i the component of the gray level of calculated value: Z i = (x si , y ri ,c ), where x si is the spatial component and y ri ,c is the range component. For the proposed segmentation algorithm, these steps are the following: 1. Run the mean shift filtering procedure for the image according to former steps, and store all the information about the d-dimensional convergence point in Zi . 2. Define the regions, which are found in the spatial domain (h s ) with intensities smaller or equal than h r /2. 3. For each region, which was defined in the step 2, look for all the pixels belonging the region whose mean intensity values are assigned at Zi . 4. Build the region graph through an adjacent list as follows: for each region, look for all adjacent regions to the right or under the region. 5. While there are nodes in the graph, which have been not visited, run a variant of the “Depth-First Search” (DFS), which allows
concatenating adjacent regions; as parameter we have to provide a node which has been not yet visited. To do this the next steps are performed: 5.1. Mark the mentioned node as visited. 5.2. While there is children not yet visited: children = current child if |Zparent − Zchild | < = hr , then, The region of the child is fused with the region of the father, where the same label is assigned and it is marked as visited. 5.2. To the parent is assigned the children of the visited child from the first position and the children which have been not visited are remained. 5.3. Look for the following child and come back to step 2 while stability it is not obtained (i.e., no more pixel values are modified). 6. Eliminate spatial regions containing less than M pixels, since those regions are considered irrelevant. 7. (Binary): To all the pixels belonging to background, assign the white color to the objects and the black color to the background. It is important to point out that the proposed segmentation algorithm for the mean shift is different to the one reported in [12,13]. In other words, observe that in [13] the steps of the segmentation algorithm are totally different to the one reported in this work. Note that while in [13], the authors work with basins of attraction of the corresponding convergence points; in this paper, we make use of children and parents. Independently, of the goal of the comparison, this was the principal novelty of the use of the mean shift in this paper. 3.2. Spectral methods As have been pointed, spectral methods are based on the calculation of the eigenvectors and eigenvalues of a matrix derived from the affinities of the pixels in the image. Normalized cut is a spectral method that measures both, total dissimilarity between the different groups as well as the total similarity within groups. The main difficulty of this method is the high dimension of the eigenproblem to be solved for the case of natural images. It does not allow using standard techniques. A possible approach in this context is considering the image as a graph only locally connected imposing null values for the affinity between unclosed pixels. It guarantees a sparse structure for the affinity matrix and a considerable reduction of the time and memory computational cost of the problem. However, it is not clear how the neglected affinities affect the image segmentation and techniques that involve dense affinity matrix promise better results. A spectral technique where
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
all affinities are retained was proposed in [16]. The method is an extension of the Nystrom approximation for the integral eigenvalue problem. The main idea is to solve the eigenvalue problem for a small random subset of pixels (sample) and then extrapolating this solution to the full set of pixels. 3.2.1. Normalized cut and Nystrom method Let G be a set of pixels and let W be a weighted affinity matrix, then the method of normalized cut divides G into two subsets G 1 and G 2 such that the normalized cut between them is minimum. The normalized cut is defined as 2cut(G 1 , G 2 ) Ncut(G 1 , G 2 ) = (7) vol(G 1 )vol(G 2 ) where cut(G 1 , G 2 ) is the sum of the weights between G 1 and G 2 , vol(G 1 ) and vol(G 2 ) is the sum of the degrees within the sets G 1 and G 2 , respectively, and the symbol “” denotes the harmonic mean. The solution of this minimization problem can be approached by solving the standard eigenvalue problem (D −1/2 WD−1/2 )V = V
(10)
where the dimensions of the matrices A and B are n × n and (n × (N − n)), respectively, n>N . If A is a positive-definite matrix, the second largest eigenvector of Wˆ can be estimated by calculating the second column of the V matrix, which is defined in the following way:
A −1/2 V= A−1/2 US S (11) BT where US and S contain the spectral decomposition of the matrix, S = A + A−1/2 B B T A−1/2
Note that the eigenproblem dimension is n. It implies that if n is taken sufficiently small, the problem can now be solved in a suitable computation time. It is clear to see that an interesting aspect of the algorithm is to know when the second largest eigenvector of Wˆ is a good approximation of the corresponding eigenvector of W. Although both vectors are theoretically different, the structure of the affinity matrix for natural images permits that generally B T A−1 B is a good approximation for matrix C, and the obtained segmentation by the method gives suitable results. The numerical experimentation carried out in [17] using the Corel database, has roughly shown that by choosing at random 100 samples were sufficient to capture the salient groups in typical natural images. In this case, the elements of the used affinity matrix have the form, 2 2 Wi j = e−Ii −I j / I ∗ e−X i −X j | / X
(13)
where X i is the spatial localization and Ii is the gray level of the ith node.
(8)
where D is a diagonal matrix, Dii = j Wi j and D −1/2 denotes the positive symmetric square root of the inverse of D. Basically, the ith component of the second largest eigenvector determines the embedding of the ith pixel. This process partitions the image in two disjointed sets. Since the dimension N of W is the square of the G dimension, the solution of (8) involves a serious numerical difficulty when the problem is solved by a conventional method. The Nystrom extension is an efficient technique that uses the information of a subset of pixels in order to obtain an approximate solution of the eigenvalue problem. The extension of Nystrom is an efficient technique that uses the information of a subset of pixels in order to obtain an approximate solution of the eigenproblem. For the sake of simplicity, the sample is defined as the first n pixels of the image matrix. The strategy consists of dividing the W matrix of the following way:
A B (9) W= BT C and approximating the matrix W by a low-rank matrix,
A B Wˆ = B T B T A−1 B
935
(12)
4. Experimental results, analysis and discussion 4.1. The method of evaluation In general, a priori, physicians know which structures they want to isolate. For that reason and particularly in this case, manual segmentation is known to provide yet the best and most reliable results. However, this task is very tedious and time consuming for the user, and thus it does not serve the needs of daily clinical use well. Up to now and due to the lack of ground truth, the quantitative evaluation of a segmentation method is difficult to achieve. An alternative is to use manual segmentation results as the ground truth. In order to evaluate the performance of both techniques, we have calculated the percentage of false negatives (FN, objects, which are not found by the strategy) and the false positives (FP, noise, which is classified as objects). These were defined according to the following expressions FP =
fp ∗ 100 Vp + f p
(14a)
FN =
fn ∗ 100 Vp + f n
(14b)
where Vp is the actual number quantity of objects identified by the physician, f n is the quantity of objects, which was not marked by the strategy and f p is the number of spurious regions, which was marked as objects. 4.2. Experimental results Now, we shall expose the results of the blood vessel segmentation obtained by using the mean shift and the spectral method. Calculations for both methods were carried out by means of MATLAB system (version 7.0). In the case of the
936
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
Fig. 3. (a) Original image, (b) binary image by using the mean shift, (c) result obtained via the spectral method.
Fig. 4. (a) Original image, (b) segmentation via the mean shift procedure, (c) obtained segmentation via the spectral method.
Fig. 5. (a) original image, (b) segmentation through the mean shift, (c) obtained segmentation with the spectral method.
spectral method, all computational methods were programmed by the authors and we have included the modifications suggested in [19]. The size of the used sample was n = 30. Fig. 3 represents the first example of this comparison. From Fig. 3, one can observe that the two methods were robust in segmenting a histological image (in this case, BV; the methods were also tested with other biomedical images). According to opinion of the physicians both techniques were able to discriminate the BV, which has a lot of importance from the point of view of the diagnosis. However, one can note that the obtained results by the mean shift were less noisy
(see Fig. 3(b)). This is due to the performance of the used algorithm in the mean shift procedure. In further work, we will implement the filtering of the mean shift in order to improve performance of the spectral method. These investigations will be published in future works. In Figs. 4–10 other results are shown that were obtained for the BV segmentation by using the mean shift and the spectral method. In the presented examples, the reader can appreciate that both methods were able to effectively well isolate the BV. However, one also can see that the segmented images via the mean shift
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
937
Fig. 6. (a) original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method. The circles indicate the differences obtained between the two methods.
Fig. 7. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
Fig. 8. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
algorithm were less noisy than those obtained via the spectral method. In addition, the edges were more robust and better defined in the images segmented with the mean shift (see Figs. 6–9). In order to obtain less noisy images with the spectral method it would be necessary to carry out an additional step of filtering with the aim of eliminating the noise that arises during the segmentation process. This procedure could be carried out through a morphological technique, which can be an opening or a majority filter. Another way could be to use the mean shift procedure as a pre-processing in order to improve
performance of the spectral method. This will be future work. It is necessary to point out that this behavior was the same for all the images that were processed. It is important to mention that all results obtained with both methods were superior than those reported in [11]. 4.3. Quantitative comparison between the two techniques From the theoretical point of view the spectral methods have bigger complexity. The mean shift is simply an average and a subtraction with regard to the point x. However, for the case
938
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
Fig. 9. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
Fig. 10. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
of the spectral methods it was necessary to obtain an approximate solution of the eigenproblem; that is, it was necessary to calculate the eigenvectors and eigenvalues of a matrix derived from the affinities of the pixels in the image. So, the spectral method was more demanding from the numeric point of view. From the point of view of the algorithm, the mean shift was less complex than the spectral method and the segmentation process was carried out in a much faster way. Indeed, in a Pentium IV of 2.6 GHz, segmentation by the spectral method took 20 min, while by means of the mean shift only 2 min was required. Since much of the motivation for biomedical image segmentation is to automate all or part of the manual segmentation, it is more important to compare the results of the segmentation with those obtained by manual segmentation. In Fig. 11 three examples of the computation of the XOR between the manual segmentation and the segmentation by using both methods are shown. Manual segmentation of an image was carried out by several pathologists for they know which objects they want to isolate. In practice this is a very tedious task. In Fig. 11 one can observe that the noise is bigger in the images segmented by means of the spectral method. This noise is associated with the FP; that is, noise classified as objects. In order to evaluate the performance of both methods, the percentage of FP and FN was calculated between the result
from the segmentation and from the manual segmentation. Numerical results of the comparison, by using expression (14), were summarized in Table 1 (obtained segmentation by the mean shift) and Table 2 (obtained segmentation by the spectral method). As an example, we presented only three of the analyzed images. Nevertheless, we carried out this comparison with all images. In Tables 1 and 2, one can observe that the percentage of error for FN was equal to 0%; that is, all regions belonging to the BV were detected. This denoted the correct performance of both methods. This behavior was the same for the 80 segmented images. Nevertheless, in Figs. 11(b) and (c) one can notice that the segmentation process was not completely correct. It can be seen, according to the judgment of physicians that some FP arose. The FP may arise due to transition of illumination or a change of contrast. However, it was evident that the number of FP in the segmented images by the spectral method was practically four times with regard to the number obtained via the mean shift method. In practice, these problems are very difficult to completely eliminate in real images. Then, the result from the full automatic segmentation steps was presented to the physicians for manual segmentation. With a few mouse clicks on the final result, the FP were completely eliminated. Images segmented by using the mean shift required less manual post-processing than images segmented by the spectral method.
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
939
Fig. 11. (a) Manually segmented images, (b) XOR with the segmented images by using the mean shift, (c) XOR with the segmented images by using the spectral method (see Figs. 3–5).
Table 1 Numerical results of comparison between manual and automatic segmentation (mean shift)
Table 2 Numerical results of comparison between manual and automatic segmentation (spectral method)
Images
Vp
fp
fn
FN (%)
FP (%)
Images
Vp
fp
fn
FN (%)
FP (%)
Figs. 11(a) and (b) (upper) Figs. 11(a) and (b) (middle) Figs. 11(a) and (b) (lower)
19 7 9
2 1 2
0 0 0
0 0 0
9.5 12.5 18.1
Figs. 11(a) and (c) (upper) Figs. 11(a) and (c) (middle) Figs. 11(a) and (c) (lower)
19 7 9
19 20 11
0 0 0
0 0 0
50. 74.0 55
5. Conclusions In this work, we have carried out a comparison between two robust segmentation techniques; namely, the mean shift method and the spectral method. Through several experiments with real images, we have demonstrated that the segmented images obtained by using the mean shift were less noisy than those obtained by means of the spectral method. In addition, the edges obtained by using the mean shift were better defined. The results obtained by both methods were compared with the manual segmentation performed by the practitioners, where
with the mean shift segmentation errors were less than 20% for false positives (FP) and 0% for false negatives (FN) were observed. On the other hand, with the spectral method more than 45% for FP and 0% for FN were obtained. We have found that the algorithm of the mean shift was 10 times faster than the spectral method. Conflict of interest statement None declared.
940
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
Acknowledgments The authors would like to thank Gert Wollny from the Max-Planck Institute of Evolutionary Anthropology (Leipzig, Germany) for discussions about the paper and grammar of English language. Authors thank CONACYT for the economical support under project number 46805. References [1] W. Kenong, D. Gauthier, M.D. Levine, Live cell image segmentation, IEEE Trans. Biomed. Eng. 42 (1) (1995). [2] J. Sijbers, P. Scheunders, M. Verhoye, A. Van der Linden, D. Van Dyck, E. Raman, Watershed-based segmentation of 3D MR data for volume quantization, Magn. Reson. Imaging 15 (6) (1997) 679–688. [3] C. Chin-Hsing, J. Lee, J. Wang, C.W. Mao, Color image segmentation for bladder cancer diagnosis, Math. Comput. Modeling 27 (2) (1998) 103–120. [4] R. Rodríguez, T. Alarcón, R. Wong, L. Sanchez, Color segmentation applied to study of the angiogenesis. Part I, J. Intell. Robotic Syst. 34 (1) (2002). [5] P. Schmid, Segmentation of digitized dermatoscopic images by twodimensional color clustering, IEEE Trans. Med. Imaging 18 (2) (1999). [6] J.E. Koss, F.D. Newman, T.K. Johnson, D.L. Kirch, Abdominal organ segmentation using texture transforms and a Hopfield neural network, IEEE Trans. Med. Imaging 18 (7) (1999). [7] S. Chaudhuri, S. Chatterjee, N. Katz, M. Nelson, M. Goldbaum, Detection of blood vessels in retinal images using two-dimensional matched filters, IEEE Trans. Med. Imaging 8 (3) (1989) 263–269. [8] I. Liu, Y. Sun, Recursive tracking of vascular networks in angiograms based on the detection-deletion scheme, IEEE Trans. Med. Imaging 12 (2) (1993) 335–341. [9] Y. Masutani, T. Schiemann, K. Hohne, Vascular shape segmentation and structure extraction using a shape-based region-growing model, Proc. MICCAI’98 (1998) 1242–1249. [10] A. Hoover, V. Kouznetsova, M. Goldbaum, Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response, IEEE Trans. Med. Imaging 19 (2000) 201–210.
[11] R. Rodríguez, T.E. Alarcón, J.J. Abad, Blood vessel segmentation via neural network in histological images, J. Intell. Robotic Syst. 36 (4) (2003). [12] D.I. Comaniciu, Nonparametric robust method for computer vision, Ph.D. Thesis, The State University of New Jersey, New Brunswick, Rutgers, 2000. [13] D. Comaniciu, P. Meer, Mean shift: a robust approach toward feature space analysis, IEEE Trans. Pattern Anal. Mach. Intell. 24 (5) (2002). [14] K. Fukunaga, L.D. Hostetler, The estimation of the gradient of a density function, IEEE Trans. Inf. Theory 21 (1975) 32–40. [15] Y. Cheng, Mean shift, mode seeking, and clustering, IEEE Trans. Pattern Anal. Mach. Intell. 17 (8) (1995) 790–799. [16] J. Shi, J. Malik, Normalized cuts and image segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 22 (8) (2000) 888–905. [17] C. Fowlkes, S. Belongie, F. Chung, J. Malik, Spectral grouping using the Nystrom method, IEEE Trans. Pattern Anal. Mach. Intell. 26 (2) (2004) 214–225. [18] R. Rodríguez, T. Alarcón, L. Sanchez, MADIP: Morphometrical Analysis by Digital Image Processing, in: Proceedings of the Ninth Spanish Symposium on Pattern Recognition and Image Analysis, vol. I, 2001, pp. 291–298. [19] E. Izquierdo, V. Guerra, Numerical stability of Nystrom extension for image segmentation, in: Proceeding of 2004 IEEE International Workshop on Machine Learning for Signal Processing, Sao Luis, Brazil, 2004. Juan Humberto Sossa Azuela was born in Guadalajara, Jalisco, Mexico in 1956. He received a Master degree in Electrical engineering from CINVESTAV-IPN, Mexico in 1987 and a Ph.D. in Informatics from the INPG, Grenoble, France, in 1992. Since 1993, he has been with the National Polytechnical Institute of Mexico, where he is a Full Time Professor in Computer Science. His research interests focus on image analysis, pattern recognition and visual control of robots. Roberto Rodríguez Morales was born in Havana, Cuba. He received a Diploma in Physics from the Physic Faculty, Havana University in 1978 and a Ph.D. from the Technical University of Havana in 1995. Since 1998, he is the Head of the Digital Signal Processing Group of the Institute of Cybernetics, Mathematics & Physics (ICIMAF). His research interests include segmentation, restoration, mathematical morphology, visual pattern recognition, analysis and interpretation of images, theoretical studies of Gaussian scale-space, new Lexis on color via Santana method and mean shift.
A comparison between two robust techniques for segmentation of blood vessels Roberto Rodríguez a , Patricio J. Castillo a , Valia Guerra b , Juan Humberto Sossa Azuela c,∗,1 , Ana G. Suáreza a , Ebroul Izquierdo d a Digital Signal Processing Group, Institute of Cybernetics, Mathematics & Physics (ICIMAF), CP 10 400, Havana, Cuba b Numerical Methods Group, Institute of Cybernetics, Mathematics & Physics (ICIMAF), CP 10 400, Havana, Cuba c Center for Computing Research (CIC), National Polytechnic Institute (IPN), Mexico d Department of Electronic Engineering, Queen Mary, University of London, UK
Received 22 August 2007; accepted 18 June 2008
Abstract Image segmentation plays an important role in image analysis. According to several authors, segmentation terminates when the observer’s goal is satisfied. For this reason, a unique method that can be applied to all possible cases does not yet exist. In this paper, we have carried out a comparison between two current segmentation techniques, namely the mean shift method, for which we propose a new algorithm, and the so-called spectral method. In this investigation the important information to be extracted from an image is the number of blood vessels (BV) present in the image. The results obtained by both strategies were compared with the results provided by manual segmentation. We have found that using the mean shift segmentation an error less than 20% for false positives (FP) and 0% for false negatives (FN) was observed, while for the spectral method more than 45% for FP and 0% for FN were obtained. We discuss the advantages and disadvantages of both methods. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Image segmentation; Mean shift; Spectral methods
1. Introduction Nowadays medical imaging provides major aid in many branches of medicine; it enables and facilitates the capture, transmission and analysis of medical images as well as providing assistance towards medical diagnoses. Medical imaging is still on the rise with new imaging modalities being added and continuous improvements of devices’ capabilities. In particular, segmentation and contour extraction are important steps in many medical imaging systems. Segmented images are now routinely used in a multitude of different applications, such as diagnosis, treatment planning, robotics, pathology, study of ∗ Corresponding author. Tel.: +52 55 57 29 60 00x56505; fax: +52 55 57 29 60 00x56607. E-mail address: [email protected] (J.H.S. Azuela). 1 Centro de Investigación en Computación-IPN, Av. Juan de Dios Bátiz s/n, Esquina con Miguel Othón de Mendizábal, Colonia Nueva Industrial Vallejo, C. P. 07700, México, D. F. MÉXICO.
0010-4825/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2008.06.002
anatomical structure and computer-integrated surgery. Image segmentation remains, however, a difficult task due to, in the one hand, the variability of the shape of the objects and, in the other hand, the variation in image quality. Medical images are often corrupted by noise and sampling artifacts, which can cause considerable difficulties when applying rigid methods. Nevertheless, in spite of the most complex algorithms developed until now, segmentation continues to be very dependent on the application, and a single method that can solve all the problems that can be encountered does not exist. Pathological anatomy is a speciality where the use of different digital image processing technique allows improving the accuracy of the diagnosis of many diseases. One of those pathologies is cancer, which is known to be one of the first causes of death in many countries. For this reason, the study of this disease by efficient methods of image processing is of supreme importance. Many segmentation methods have been proposed for medical-image data [1–6]. Unfortunately, segmentation using
932
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
traditional low-level image processing techniques, such as thresholding, region growing and other classical operations requires a considerable amount of interactive guidance in order to attain satisfactory results. Automating these model-free approaches is difficult because of complexity, shadows, and variability within and across individual objects. Furthermore, noise and other image artifacts can cause incorrect regions or boundary discontinuities in objects recovered from these methods. Many segmentation methods have been proposed for the study of blood vessels (BV) [7–11]. The main goal of this paper is to carry out a comparison between two segmentation methods: the mean shift method and the spectral method. In the case of the mean shift method we proposed a new segmentation algorithm. Our interest in comparing these two techniques is due to the good results that both have offered in other segmentation fields. A goal of this research is to determine which of the methods is the most suitable, robust and faster for segmenting images containing BV. In this application the important information to be extracted from images is just the number of BV. Mean shift is a nonparametric and versatile tool for feature analysis and can provide reliable solutions for many vision tasks. Refer for example to [12,13]. The mean shift was proposed in 1975 by Fukunaga and Hostetler [14] and largely forgotten until Cheng’s paper [15] rekindled interest in it. The spectral method was the other technique to be compared. It is based on the calculation of the eigenvectors and eigenvalues of an obtained matrix from the affinities of the pixels within the image. In particular, a spectral method was used which is called of normalized cut [16,17]. This method measures both, the total dissimilarity between the different groups as well as the total similarity within groups. The main idea is to solve the eigenvalue problem for a small subset of pixels (sample) and then this solution is extrapolated to the full set of pixels. The theoretical aspects of the used methods are exposed later on. The remaining of the paper is organized as follows. In Section 2, the characteristics of the images under study are described. In Section 3, the details of the most significant theoretical aspects of each of the two methods are presented. In Section 4, the steps in the experimentation are explained. Here, we carry out an analysis and discussion of both methods. In addition, the
method of evaluation is also given. Finally, in Section 5 the most important conclusions of this research are given. 2. Method and features of the studied images As biomedical images, we decided to use biopsies, which represent histological cuts in malignant tumors. These biopsies were included in paraffin by using the immunohistoquimic technique with the complex method of avidina biotina. Finally, monoclonal CD34 was contrasted with green methyl in order to accentuate new formation of BV. These biopsies were obtained from soft parts of human bodies. This analysis was carried out for more than 80 patients. In Fig. 1 typical images are given. These images were captured via a Morphometrical Analysis by Digital Image Processing (MADIP) system with a resolution of 512 × 512 × 8 bit/pixels. MADIP was a system created to help physicians in their research work using images [18]. MADIP is a software for the morphometrical characterization of macro and micro lesions in biomedical images. In general, it may be used in any branch of science and technology, where not only quantitative techniques of high precision are required, but it is also necessary to do a big number of measurements with high accuracy and rapidly. In Fig. 2 one can be observe a horizontal profile through the center of a vessel; that is, a plot of the intensities of pixels along a single row. By carefully observing this kind of images, we have found several characteristics of these images, which are common to typical images that we encounter in the tissues of biopsies: 1. The intensity was slightly darker within the blood vessel than in the local surrounding background. It was emphasized that this observation held only within the local surroundings. 2. High local variation of intensity was observed both within the blood vessel and the background. However, the local variation of intensity was higher within the blood vessel than in background regions (see Fig. 2(b)). 3. The variability of BV both, in size and shape can be observed (see Fig. 1). It is important to point out that due to the adopted acquisition protocol the images were corrupted with a lot of noise.
Fig. 1. These images represent the blood vessels in malignant tumors. Some blood vessels are marked with arrows.
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
933
Fig. 2. A profile of intensity through the center of a vessel. Profile is indicated by a line in Fig. 2(a).
The quantity n x /n(h d cd ) is the kernel density estimate fˆU (x) (where U means the uniform kernel) computed with the hypersphere Sh (x), and thus expression (2) can be written as
3. Theoretical considerations 3.1. About the mean shift method The iterative procedure to compute the mean shift is introduced as a normalized density estimate of the gradient [13]. By employing a differentiable kernel, an estimate of the density gradient can be defined as the gradient of the kernel density estimate; that is ∇ˆ f (x) = ∇ fˆ(x) =
n 1
nhd
∇K
i=1
x − xi h
1/2cd−1 (d + 2)(1 − x2 ) if x < 1 0 otherwise
K E (x) =
The density gradient estimate becomes ∇ˆ f E (x) =
1 d +2 · n(h d cd ) h2
(xi − x)
xi ∈Sh (x)
⎛ nx d +2 ⎝ 1 = · n(h d cd ) h2 nx
⎞ (xi − x)⎠
(2)
xi ∈Sh (x)
where the region Sh (x) is a hypersphere of radius h having volume h d cd , centered at x, and containing n x data points, that is, the uniform kernel. In addition, in this case d = 3, due to vector x is a three-dimensional vector; it is of dimension two for the spatial domain and one for the range domain (gray levels). The last term in expression (2) is called the sample mean shift, Mh,U (x) =
1 nx
(xi − x) =
xi ∈Sh (x)
1 nx
xi ∈Sh (x)
(4)
which yields Mh,U (x) =
h 2 ∇ˆ f E (x) d + 2 fˆU (x)
(5)
(1)
Conditions on the kernel K (·) and the window radio h are derived in [14] to guarantee asymptotic unbiasedness, meansquare consistency, and uniform consistency of the estimate in expression (1). For example, for Epanechikov kernel
d +2 Mh,U (x) ∇ˆ f E (x) = fˆU (x) · h2
xi − x
(3)
Eq. (5) shows that an estimate of the normalized gradient can be obtained by computing the sample mean shift in a uniform kernel centered on x. In addition, the mean shift has the direction of the gradient of the density estimate at x when this estimate is obtained with the Epanechnikov kernel. Since the mean shift vector always points towards the direction of the maximum increase in the density, it can define a path leading to a local density maximum; that is, to a mode of the density. In addition, expression (5) shows that the mean is shifted towards the region in which the majority of the points reside. Since the mean shift is proportional to the local gradient estimate, it can define a path leading to a stationary point of the estimated density, where these stationary points are the modes. Moreover, as it was pointed out, the normalized gradient in expression (5) introduces a desirable adaptive behavior, since the mean shift step is large for low-density regions corresponding to valleys, and decreases as x approaches a mode [13]. In [13] it was proven that the mean shift procedure, obtained by successive: • computing the mean shift vector Mh (x), • translating the window Sh (x) by Mh (x), guarantees convergence. It is well known that a digital image can be represented as a two-dimensional array of p-dimensional vectors (pixels), where p = 1 in the gray level case, three for color images, and p > 3 in the multispectral case.
934
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
As was pointed in [13] that when the location and range vectors are concatenated in the joint spatial-range domain of dimension d= p+2, their different nature has to be compensated by proper normalization of parameters h s and h r . Thus, the multivariable kernel is defined as the product of two radially symmetric kernels and the Euclidean metric allows a single bandwidth for each domain, that is
r 2 x s 2 C k y (6) K h s ,h r (x) = 2 p k h h hs hr s r where x s is the spatial part, y r is the range part of a feature vector, k(x) the common profile used in both domains, h s and h r , the employed kernel bandwidths, and C the corresponding normalization constant. Basically, the algorithm that we used in this work is performed into two steps, a filtering step and a segmentation step. The filtering algorithm comprises the following steps [13]: Let {X i }i and {Z i }i , i = 1, . . . , n be the input and filtered images in the joint spatial-range domain. For each pixel p ∈ X i , p = (x , y , z ) ∈ R3 , where (x , y ) ∈ R2 and z ∈ [0, 2 − 1], being the quantity of bits/pixel in the image. Then, for eachi = 1, . . . , n . 1. Initialize j = 1 and y i ,1 = x i . 2. Compute through the mean shift y i ,j +1 (see expression (3)), the mode where the pixel converges; that is, the calculation of the mean shift is carried out until convergence, y = y i ,c . 3. Store at Z i the component of the gray level of calculated value: Z i = (x si , y ri ,c ), where x si is the spatial component and y ri ,c is the range component. For the proposed segmentation algorithm, these steps are the following: 1. Run the mean shift filtering procedure for the image according to former steps, and store all the information about the d-dimensional convergence point in Zi . 2. Define the regions, which are found in the spatial domain (h s ) with intensities smaller or equal than h r /2. 3. For each region, which was defined in the step 2, look for all the pixels belonging the region whose mean intensity values are assigned at Zi . 4. Build the region graph through an adjacent list as follows: for each region, look for all adjacent regions to the right or under the region. 5. While there are nodes in the graph, which have been not visited, run a variant of the “Depth-First Search” (DFS), which allows
concatenating adjacent regions; as parameter we have to provide a node which has been not yet visited. To do this the next steps are performed: 5.1. Mark the mentioned node as visited. 5.2. While there is children not yet visited: children = current child if |Zparent − Zchild | < = hr , then, The region of the child is fused with the region of the father, where the same label is assigned and it is marked as visited. 5.2. To the parent is assigned the children of the visited child from the first position and the children which have been not visited are remained. 5.3. Look for the following child and come back to step 2 while stability it is not obtained (i.e., no more pixel values are modified). 6. Eliminate spatial regions containing less than M pixels, since those regions are considered irrelevant. 7. (Binary): To all the pixels belonging to background, assign the white color to the objects and the black color to the background. It is important to point out that the proposed segmentation algorithm for the mean shift is different to the one reported in [12,13]. In other words, observe that in [13] the steps of the segmentation algorithm are totally different to the one reported in this work. Note that while in [13], the authors work with basins of attraction of the corresponding convergence points; in this paper, we make use of children and parents. Independently, of the goal of the comparison, this was the principal novelty of the use of the mean shift in this paper. 3.2. Spectral methods As have been pointed, spectral methods are based on the calculation of the eigenvectors and eigenvalues of a matrix derived from the affinities of the pixels in the image. Normalized cut is a spectral method that measures both, total dissimilarity between the different groups as well as the total similarity within groups. The main difficulty of this method is the high dimension of the eigenproblem to be solved for the case of natural images. It does not allow using standard techniques. A possible approach in this context is considering the image as a graph only locally connected imposing null values for the affinity between unclosed pixels. It guarantees a sparse structure for the affinity matrix and a considerable reduction of the time and memory computational cost of the problem. However, it is not clear how the neglected affinities affect the image segmentation and techniques that involve dense affinity matrix promise better results. A spectral technique where
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
all affinities are retained was proposed in [16]. The method is an extension of the Nystrom approximation for the integral eigenvalue problem. The main idea is to solve the eigenvalue problem for a small random subset of pixels (sample) and then extrapolating this solution to the full set of pixels. 3.2.1. Normalized cut and Nystrom method Let G be a set of pixels and let W be a weighted affinity matrix, then the method of normalized cut divides G into two subsets G 1 and G 2 such that the normalized cut between them is minimum. The normalized cut is defined as 2cut(G 1 , G 2 ) Ncut(G 1 , G 2 ) = (7) vol(G 1 )vol(G 2 ) where cut(G 1 , G 2 ) is the sum of the weights between G 1 and G 2 , vol(G 1 ) and vol(G 2 ) is the sum of the degrees within the sets G 1 and G 2 , respectively, and the symbol “” denotes the harmonic mean. The solution of this minimization problem can be approached by solving the standard eigenvalue problem (D −1/2 WD−1/2 )V = V
(10)
where the dimensions of the matrices A and B are n × n and (n × (N − n)), respectively, n>N . If A is a positive-definite matrix, the second largest eigenvector of Wˆ can be estimated by calculating the second column of the V matrix, which is defined in the following way:
A −1/2 V= A−1/2 US S (11) BT where US and S contain the spectral decomposition of the matrix, S = A + A−1/2 B B T A−1/2
Note that the eigenproblem dimension is n. It implies that if n is taken sufficiently small, the problem can now be solved in a suitable computation time. It is clear to see that an interesting aspect of the algorithm is to know when the second largest eigenvector of Wˆ is a good approximation of the corresponding eigenvector of W. Although both vectors are theoretically different, the structure of the affinity matrix for natural images permits that generally B T A−1 B is a good approximation for matrix C, and the obtained segmentation by the method gives suitable results. The numerical experimentation carried out in [17] using the Corel database, has roughly shown that by choosing at random 100 samples were sufficient to capture the salient groups in typical natural images. In this case, the elements of the used affinity matrix have the form, 2 2 Wi j = e−Ii −I j / I ∗ e−X i −X j | / X
(13)
where X i is the spatial localization and Ii is the gray level of the ith node.
(8)
where D is a diagonal matrix, Dii = j Wi j and D −1/2 denotes the positive symmetric square root of the inverse of D. Basically, the ith component of the second largest eigenvector determines the embedding of the ith pixel. This process partitions the image in two disjointed sets. Since the dimension N of W is the square of the G dimension, the solution of (8) involves a serious numerical difficulty when the problem is solved by a conventional method. The Nystrom extension is an efficient technique that uses the information of a subset of pixels in order to obtain an approximate solution of the eigenvalue problem. The extension of Nystrom is an efficient technique that uses the information of a subset of pixels in order to obtain an approximate solution of the eigenproblem. For the sake of simplicity, the sample is defined as the first n pixels of the image matrix. The strategy consists of dividing the W matrix of the following way:
A B (9) W= BT C and approximating the matrix W by a low-rank matrix,
A B Wˆ = B T B T A−1 B
935
(12)
4. Experimental results, analysis and discussion 4.1. The method of evaluation In general, a priori, physicians know which structures they want to isolate. For that reason and particularly in this case, manual segmentation is known to provide yet the best and most reliable results. However, this task is very tedious and time consuming for the user, and thus it does not serve the needs of daily clinical use well. Up to now and due to the lack of ground truth, the quantitative evaluation of a segmentation method is difficult to achieve. An alternative is to use manual segmentation results as the ground truth. In order to evaluate the performance of both techniques, we have calculated the percentage of false negatives (FN, objects, which are not found by the strategy) and the false positives (FP, noise, which is classified as objects). These were defined according to the following expressions FP =
fp ∗ 100 Vp + f p
(14a)
FN =
fn ∗ 100 Vp + f n
(14b)
where Vp is the actual number quantity of objects identified by the physician, f n is the quantity of objects, which was not marked by the strategy and f p is the number of spurious regions, which was marked as objects. 4.2. Experimental results Now, we shall expose the results of the blood vessel segmentation obtained by using the mean shift and the spectral method. Calculations for both methods were carried out by means of MATLAB system (version 7.0). In the case of the
936
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
Fig. 3. (a) Original image, (b) binary image by using the mean shift, (c) result obtained via the spectral method.
Fig. 4. (a) Original image, (b) segmentation via the mean shift procedure, (c) obtained segmentation via the spectral method.
Fig. 5. (a) original image, (b) segmentation through the mean shift, (c) obtained segmentation with the spectral method.
spectral method, all computational methods were programmed by the authors and we have included the modifications suggested in [19]. The size of the used sample was n = 30. Fig. 3 represents the first example of this comparison. From Fig. 3, one can observe that the two methods were robust in segmenting a histological image (in this case, BV; the methods were also tested with other biomedical images). According to opinion of the physicians both techniques were able to discriminate the BV, which has a lot of importance from the point of view of the diagnosis. However, one can note that the obtained results by the mean shift were less noisy
(see Fig. 3(b)). This is due to the performance of the used algorithm in the mean shift procedure. In further work, we will implement the filtering of the mean shift in order to improve performance of the spectral method. These investigations will be published in future works. In Figs. 4–10 other results are shown that were obtained for the BV segmentation by using the mean shift and the spectral method. In the presented examples, the reader can appreciate that both methods were able to effectively well isolate the BV. However, one also can see that the segmented images via the mean shift
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
937
Fig. 6. (a) original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method. The circles indicate the differences obtained between the two methods.
Fig. 7. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
Fig. 8. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
algorithm were less noisy than those obtained via the spectral method. In addition, the edges were more robust and better defined in the images segmented with the mean shift (see Figs. 6–9). In order to obtain less noisy images with the spectral method it would be necessary to carry out an additional step of filtering with the aim of eliminating the noise that arises during the segmentation process. This procedure could be carried out through a morphological technique, which can be an opening or a majority filter. Another way could be to use the mean shift procedure as a pre-processing in order to improve
performance of the spectral method. This will be future work. It is necessary to point out that this behavior was the same for all the images that were processed. It is important to mention that all results obtained with both methods were superior than those reported in [11]. 4.3. Quantitative comparison between the two techniques From the theoretical point of view the spectral methods have bigger complexity. The mean shift is simply an average and a subtraction with regard to the point x. However, for the case
938
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
Fig. 9. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
Fig. 10. (a) Original image, (b) segmentation via the mean shift, (c) segmentation via the spectral method.
of the spectral methods it was necessary to obtain an approximate solution of the eigenproblem; that is, it was necessary to calculate the eigenvectors and eigenvalues of a matrix derived from the affinities of the pixels in the image. So, the spectral method was more demanding from the numeric point of view. From the point of view of the algorithm, the mean shift was less complex than the spectral method and the segmentation process was carried out in a much faster way. Indeed, in a Pentium IV of 2.6 GHz, segmentation by the spectral method took 20 min, while by means of the mean shift only 2 min was required. Since much of the motivation for biomedical image segmentation is to automate all or part of the manual segmentation, it is more important to compare the results of the segmentation with those obtained by manual segmentation. In Fig. 11 three examples of the computation of the XOR between the manual segmentation and the segmentation by using both methods are shown. Manual segmentation of an image was carried out by several pathologists for they know which objects they want to isolate. In practice this is a very tedious task. In Fig. 11 one can observe that the noise is bigger in the images segmented by means of the spectral method. This noise is associated with the FP; that is, noise classified as objects. In order to evaluate the performance of both methods, the percentage of FP and FN was calculated between the result
from the segmentation and from the manual segmentation. Numerical results of the comparison, by using expression (14), were summarized in Table 1 (obtained segmentation by the mean shift) and Table 2 (obtained segmentation by the spectral method). As an example, we presented only three of the analyzed images. Nevertheless, we carried out this comparison with all images. In Tables 1 and 2, one can observe that the percentage of error for FN was equal to 0%; that is, all regions belonging to the BV were detected. This denoted the correct performance of both methods. This behavior was the same for the 80 segmented images. Nevertheless, in Figs. 11(b) and (c) one can notice that the segmentation process was not completely correct. It can be seen, according to the judgment of physicians that some FP arose. The FP may arise due to transition of illumination or a change of contrast. However, it was evident that the number of FP in the segmented images by the spectral method was practically four times with regard to the number obtained via the mean shift method. In practice, these problems are very difficult to completely eliminate in real images. Then, the result from the full automatic segmentation steps was presented to the physicians for manual segmentation. With a few mouse clicks on the final result, the FP were completely eliminated. Images segmented by using the mean shift required less manual post-processing than images segmented by the spectral method.
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
939
Fig. 11. (a) Manually segmented images, (b) XOR with the segmented images by using the mean shift, (c) XOR with the segmented images by using the spectral method (see Figs. 3–5).
Table 1 Numerical results of comparison between manual and automatic segmentation (mean shift)
Table 2 Numerical results of comparison between manual and automatic segmentation (spectral method)
Images
Vp
fp
fn
FN (%)
FP (%)
Images
Vp
fp
fn
FN (%)
FP (%)
Figs. 11(a) and (b) (upper) Figs. 11(a) and (b) (middle) Figs. 11(a) and (b) (lower)
19 7 9
2 1 2
0 0 0
0 0 0
9.5 12.5 18.1
Figs. 11(a) and (c) (upper) Figs. 11(a) and (c) (middle) Figs. 11(a) and (c) (lower)
19 7 9
19 20 11
0 0 0
0 0 0
50. 74.0 55
5. Conclusions In this work, we have carried out a comparison between two robust segmentation techniques; namely, the mean shift method and the spectral method. Through several experiments with real images, we have demonstrated that the segmented images obtained by using the mean shift were less noisy than those obtained by means of the spectral method. In addition, the edges obtained by using the mean shift were better defined. The results obtained by both methods were compared with the manual segmentation performed by the practitioners, where
with the mean shift segmentation errors were less than 20% for false positives (FP) and 0% for false negatives (FN) were observed. On the other hand, with the spectral method more than 45% for FP and 0% for FN were obtained. We have found that the algorithm of the mean shift was 10 times faster than the spectral method. Conflict of interest statement None declared.
940
R. Rodríguez et al. / Computers in Biology and Medicine 38 (2008) 931 – 940
Acknowledgments The authors would like to thank Gert Wollny from the Max-Planck Institute of Evolutionary Anthropology (Leipzig, Germany) for discussions about the paper and grammar of English language. Authors thank CONACYT for the economical support under project number 46805. References [1] W. Kenong, D. Gauthier, M.D. Levine, Live cell image segmentation, IEEE Trans. Biomed. Eng. 42 (1) (1995). [2] J. Sijbers, P. Scheunders, M. Verhoye, A. Van der Linden, D. Van Dyck, E. Raman, Watershed-based segmentation of 3D MR data for volume quantization, Magn. Reson. Imaging 15 (6) (1997) 679–688. [3] C. Chin-Hsing, J. Lee, J. Wang, C.W. Mao, Color image segmentation for bladder cancer diagnosis, Math. Comput. Modeling 27 (2) (1998) 103–120. [4] R. Rodríguez, T. Alarcón, R. Wong, L. Sanchez, Color segmentation applied to study of the angiogenesis. Part I, J. Intell. Robotic Syst. 34 (1) (2002). [5] P. Schmid, Segmentation of digitized dermatoscopic images by twodimensional color clustering, IEEE Trans. Med. Imaging 18 (2) (1999). [6] J.E. Koss, F.D. Newman, T.K. Johnson, D.L. Kirch, Abdominal organ segmentation using texture transforms and a Hopfield neural network, IEEE Trans. Med. Imaging 18 (7) (1999). [7] S. Chaudhuri, S. Chatterjee, N. Katz, M. Nelson, M. Goldbaum, Detection of blood vessels in retinal images using two-dimensional matched filters, IEEE Trans. Med. Imaging 8 (3) (1989) 263–269. [8] I. Liu, Y. Sun, Recursive tracking of vascular networks in angiograms based on the detection-deletion scheme, IEEE Trans. Med. Imaging 12 (2) (1993) 335–341. [9] Y. Masutani, T. Schiemann, K. Hohne, Vascular shape segmentation and structure extraction using a shape-based region-growing model, Proc. MICCAI’98 (1998) 1242–1249. [10] A. Hoover, V. Kouznetsova, M. Goldbaum, Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response, IEEE Trans. Med. Imaging 19 (2000) 201–210.
[11] R. Rodríguez, T.E. Alarcón, J.J. Abad, Blood vessel segmentation via neural network in histological images, J. Intell. Robotic Syst. 36 (4) (2003). [12] D.I. Comaniciu, Nonparametric robust method for computer vision, Ph.D. Thesis, The State University of New Jersey, New Brunswick, Rutgers, 2000. [13] D. Comaniciu, P. Meer, Mean shift: a robust approach toward feature space analysis, IEEE Trans. Pattern Anal. Mach. Intell. 24 (5) (2002). [14] K. Fukunaga, L.D. Hostetler, The estimation of the gradient of a density function, IEEE Trans. Inf. Theory 21 (1975) 32–40. [15] Y. Cheng, Mean shift, mode seeking, and clustering, IEEE Trans. Pattern Anal. Mach. Intell. 17 (8) (1995) 790–799. [16] J. Shi, J. Malik, Normalized cuts and image segmentation, IEEE Trans. Pattern Anal. Mach. Intell. 22 (8) (2000) 888–905. [17] C. Fowlkes, S. Belongie, F. Chung, J. Malik, Spectral grouping using the Nystrom method, IEEE Trans. Pattern Anal. Mach. Intell. 26 (2) (2004) 214–225. [18] R. Rodríguez, T. Alarcón, L. Sanchez, MADIP: Morphometrical Analysis by Digital Image Processing, in: Proceedings of the Ninth Spanish Symposium on Pattern Recognition and Image Analysis, vol. I, 2001, pp. 291–298. [19] E. Izquierdo, V. Guerra, Numerical stability of Nystrom extension for image segmentation, in: Proceeding of 2004 IEEE International Workshop on Machine Learning for Signal Processing, Sao Luis, Brazil, 2004. Juan Humberto Sossa Azuela was born in Guadalajara, Jalisco, Mexico in 1956. He received a Master degree in Electrical engineering from CINVESTAV-IPN, Mexico in 1987 and a Ph.D. in Informatics from the INPG, Grenoble, France, in 1992. Since 1993, he has been with the National Polytechnical Institute of Mexico, where he is a Full Time Professor in Computer Science. His research interests focus on image analysis, pattern recognition and visual control of robots. Roberto Rodríguez Morales was born in Havana, Cuba. He received a Diploma in Physics from the Physic Faculty, Havana University in 1978 and a Ph.D. from the Technical University of Havana in 1995. Since 1998, he is the Head of the Digital Signal Processing Group of the Institute of Cybernetics, Mathematics & Physics (ICIMAF). His research interests include segmentation, restoration, mathematical morphology, visual pattern recognition, analysis and interpretation of images, theoretical studies of Gaussian scale-space, new Lexis on color via Santana method and mean shift.