Ultrasound in Med. & Biol., Vol. 31, No. 11, pp. 1485–1497, 2005 Copyright © 2005 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/05/$–see front matter
doi:10.1016/j.ultrasmedbio.2005.07.005
● Original Contribution ULTRASOUND IMAGE SEGMENTATION USING SPECTRAL CLUSTERING NECULAI ARCHIP,* ROBERT ROHLING,† PETER COOPERBERG‡ and HAMID TAHMASEBPOUR‡ *Harvard Medical School, Brigham and Women’s Hospital, Boston, MA, USA; †Department of Electrical and Computer Engineering, University of British Columbia, Vancouver, B.C., Canada; and ‡Department of Radiology, University of British Columbia, Vancouver, B.C., Canada (Received 14 Feburary 2005; in final form 7 July 2005)
Abstract—Segmentation of ultrasound images is necessary in a variety of clinical applications, but the development of automatic techniques is still an open problem. Spectral clustering techniques have recently become popular for data and image analysis. In particular, image segmentation has been proposed via the normalized cut (NCut) criterion. This article describes an initial investigation to determine the suitability of such segmentation techniques for ultrasound images. The adaptation of the NCut technique to ultrasound is described first. Segmentation is then performed on simulated ultrasound images. Tests are also performed on abdominal and fetal images with the segmentation results compared to manual segmentation. The success of the segmentation on these test cases warrants further research into NCut-based segmentation of ultrasound images. E-mail: (
[email protected]). © 2005 World Federation for Ultrasound in Medicine & Biology. Key Words: Image segmentation, Spectral clustering, Fiedler eigenvector, Graph partitioning.
example, an active contour model that uses several temporally adjacent images during the extraction of the tongue surface from a sequence of images is presented by Akgul et al. (2000). Another active contour approach (Hamarneh and Gustavsson, 2000) extracts the left ventricle of the heart in echocardiography by including a priori knowledge of ventricle shape to restrict the allowable range of deformations. Another approach uses model-based initialization and a discrete dynamic contour to segment the prostate (Ladak et al., 2000). The algorithm requires the user to select four points on the organ boundary from which an estimate of the prostate shape is interpolated by using cubic functions. The estimated contour is then used as an initialization stage for an active contour to better fit the prostate boundary. An extension of this algorithm, which uses six user selected control points to initialize a 3D active surface, is proposed by Hu et al. (2002). The active surface is then used to segment 3D ultrasound volumes of the prostate. A 3D deformable surface is also used by Ghanei and SoltanianZadeh (2002) for segmenting the prostate in 3D ultrasound. The user must first outline the prostate in several 2D cross-sections of the 3D volume data to initialize the model. Jacob et al. (2001) use snakes with a temporal Kalman filter to extract periodic motions from a sequence of ultrasound images. The algorithm performs an
INTRODUCTION The segmentation of medical images has been a challenging research topic for many years. Successful segmentation techniques can be used in a variety of clinical applications. Some example applications are the measurement of the volume of cancer tumors and cysts (Pathak et al., 2000), the planning of radiotherapy procedures (Archip et al., 2002), and the reconstruction of three-dimensional anatomical models of patients for minimally invasive surgery (Kaus et al., 2001). In clinical practice, ultrasound images are often segmented manually, but manual techniques can be laborious. More sophisticated techniques are needed. Semiautomatic segmentation techniques are a first step towards providing computer assistance, whereby the user initializes or guides the segmentation process. A common semiautomatic technique is based on active contours, such as “snakes” in 2D or “active surfaces” in 3D (Kass et al., 1988). Active contours are usually initialized near an organ boundary, and the process automatically deforms the contour toward the boundary. For Address correspondence to: Dr. Neculai Archip, Harvard Medical School, Brigham and Women’s Hospital, Thorn 329, Dept Radiology, 75 Francis St, Boston, MA, 02115, E-mail: narchip@bwh. harvard.edu 1485
1486
Ultrasound in Medicine and Biology
optimization on all the images in the sequence to extract the motion parameters. A Kalman filter is also used by Abolmaesumi and Sirouspour (2004) to improve the segmentation of organ cavities in ultrasound. Another approach (Ghiachetti, 1998) uses several levels of lowpass filtering, model-based constrained contour evolution by using snake algorithms and optical flow estimation to perform contour extraction in echocardiography. A statistical shape model is reported for aiding prostate segmentation (Shen et al., 2003). The algorithm uses a Gabor filter bank in multiple scales and orientations to characterize the prostate boundary in transrectal ultrasound. The proposed method then uses a deformable model that converges to the prostate boundary in a coarse to fine approach. Active contours are commonly used for image segmentation. The traditional weaknesses are the need for accurate initialization, correct tuning of the deformation parameters, and iteration to a local minima. The above list of papers on active contours shows various methods to improve on initialization, speed and robustness. Yet many of the successful techniques are specialized for a particular clinical area, limiting their usefulness in other applications. A wide variety of other segmentation methods have been proposed. For example, a set of morphology-based operations is used on the entire image to extract endocardial borders, by Choy and Jin (1996). Applying image-processing operators to the image can also help with organ boundary detection. In other words, a combination of image processing techniques can be used for more effective segmentation. For example, Awad et al. (2003) first used the “sticks” filter algorithm, originally reported by Pathak et al. (1998, 2000), to enhance the contrast of prostate images. Recently, an algorithm was proposed by Mudabhushi and Metaxas (2003) that uses low level image filtering, speckle reduction, probabilistic classification of image pixels based on intensity and texture as well as deformable models in order to segment breast lesions from ultrasound images. Other methods, such as texture segmentation (Muzzolini et al., 1993), adaptive temporal filtering (Evans and Nixon, 1996), neural networks (Botros, 1992), radial basis functions (Liu and Tang, 1999), maximum likelihood methods (Strintzis and Kokkindis, 1997) and anisotropic diffusion filters (Pathak et al., 2000) have also been proposed for segmentation of ultrasound. Despite these attempts, few semiautomatic ultrasound segmentation techniques have found widespread acceptance for general-purpose use. Spectral clustering techniques have recently become popular for data and image analysis in nonultrasound applications. In general terms, clustering is the process of grouping similar pixels together so the image is partitioned into separate regions. If the number of
Volume 31, Number 11, 2005
clusters is known a priori, the method is called supervised, otherwise, it is called unsupervised. Once the image is partitioned, the user simply selects the region of interest according to its cluster of pixels. The purpose of this paper is to explore the potential of spectral clustering with ultrasound. To the best of the authors’ knowledge, this approach has not been previously investigated for ultrasound. The spectral clustering techniques are based on Gestalt laws of image perception which is seen as a meaningful organization of objects in a scene. Various factors that can contribute to this process, include grouping cues such as proximity and similarity. A great deal of research in computer vision over the last few decades has sought principled ways to operationalize these ideas. An essential component is the development of grouping methodologies that use these low-level cues to perform image segmentation. A recently developed idea is to cluster pixels (or other image elements) using pairwise affinities in the form of a normalized cut (NCut) (Shi and Malik, 2000). The basic idea is to use a new graph-theoretic criterion for measuring the goodness of an image partition. The maximization of this criterion can be formulated as a generalized eigen value problem. An important drawback of this approach is the extreme computational demand due to the large matrices involved. In this article, NCut is adapted for ultrasound images and tested on a number of images using a simple speed-up technique. The study is structured as follows. First, anisotropic diffusion is proposed in combination with NCut. Next, the NCut algorithm is described and explained using a simple example. NCut is a supervised technique, meaning that the number of segmented clusters must be specified by the user. Modifications to the NCut algorithm to make it unsupervised are described next. A windowing method is then described to improve performances. Simulated ultrasound images are then tested since the true object boundary is known. Next are tests on images of liver cysts, since they are readily seen and are segmented manually. The amniotic fluid-filled regions of fetal images are segmented next. Although measurements of the segmented amniotic fluid volume in fetus images may be clinically useful, the images in this study are chosen mainly to demonstrate the segmentation technique. Comparisons are, therefore, made between the segmented clusters and the manual segmentation. A few additional tests are also included to show segmentation of the bladder and carotid artery. METHODS The NCut technique can be applied directly to the ultrasound image, but good results in MRI are reported for a combination of NCut with an anisotropic diffusion
Ultrasound spectral clustering ● N. ARCHIP et al.
1487
filter (Carballido-Gamio et al., 2004). Montagnat et al. (2003) also use the anisotropic filter successfully during segmentation of echocardiograms. The goal is to facilitate the NCut segmentation while preserving the image features. Anisotropic filtering The objective of anisotropic diffusion is to reduce the largest variations in pixel intensity due to noise or speckle. The goal is to do this without distortion of the edges, such that (at all scales) intra-region smoothing occurs preferentially over inter-region smoothing (Perona and Malik, 1990; Perona et al., 1994), and to satisfy the causality and immediate localization criteria. Given an image I, the anisotropic diffusion equation is: I ⫽ c(x, y, t)ⵜ2I ⫹ ⵜ c · ⵜ I
(1)
where c(x,y,t) is the conduction coefficient, t is the scale variable, ⵜ and ⵜ2 represent the gradient and the Laplacian operators with respect to the space variables. The objective is to assign a value of 1 to the conduction coefficient inside each region, and 0 at the boundaries. Setting the conduction coefficient as a function of the magnitude of the gradient of brightness function, as suggested by Perona and Malik (1990), provides a good estimate of the edge positions: c(x, y, t) ⫽ g共㛳ⵜI(x, y, t)㛳兲
(2)
where g should be smooth and g( · ) is restricted to a subclass of monotoically decreasing functions. In Fig. 1, an example is presented of anisotropic diffusion filter for an ultrasound image. Supervised spectral clustering using the normalized cut Spectral methods for image segmentation use the eigen vectors and eigen values of a matrix derived from the pairwise affinities of pixels. These eigen vectors induce an embedding of the pixels in a low-dimensional subspace wherein a simple central clustering method (such as k-means) can be used to do the final partitioning. Here, the partitioning is based on the normalized cut (NCut) (Shi and Malik, 2000). Constrained spectral clustering algorithms are also presented by Yu and Shi (2004). The mathematical notations are introduced first. Let I be the original image of size N ⫻ N. The 2 symmetric matrix S ⑀ RN ⫻ N2 denotes the weighted adjacency matrix for a graph G ⫽ (V, E) with nodes V representing pixels and edges E whose weights capture the pairwise affinities between pixels. Let A and B represent a bipartition of V, i.e., A 艛 B ⫽ V and A 艚 B ⫽ A. Let cut(A, B) denote the sum of the weights between A and B: cut(A, B) ⫽ ⌺i⑀A,j⑀B Sij. The degree of the ith node is defined as di ⫽ ⌺j Sij. The total connection from nodes in
Fig. 1. Ultrasound image of a carotid artery (a) before anisotropic filtering and (b) after anisotropic filtering.
the set A to all nodes in the graph is denoted by assoc(A, V) ⫽ ⌺i⑀A,j⑀V Sij. The normalized cut between sets A and B is then given by: NCut(A, B) ⫽
cut(A, B) cut(A, B) ⫹ . assoc(A, V) assoc(B, V)
(3)
The problem is to find A and B such that NCut(A, B) is minimized. Using elements of spectral graph theory (Chung, 1997), it is shown (Shi and Malik, 2000) that an approximate solution may be obtained by thresholding the eigen vector (called the Fiedler vector) corresponding to the second eigen value 2 of the generalized eigenvalue problem: (D ⫺ S)y ⫽ Dy where D is a diagonal matrix with entries Dii ⫽ di.
(4)
1488
Ultrasound in Medicine and Biology
Volume 31, Number 11, 2005
Image segmentation is reduced to the problem of partitioning the set V into disjoint sets V1,.., Vn, such that similarity among nodes in Vi is high and similarity across Vi and Vj is low. Here, we define the similarity function between two pixels i and j as:
Sij ⫽
冦
㛳F(i) ⫺ F(j)㛳2 maxF ⫺ minF if 㛳X(i) ⫺ X(j)㛳2 ⬍ r 㛳X(i) ⫺ X(j)㛳2
1⫺ 0
otherwise
冧
(5)
where X(i) is the spatial location of node i and F(i) is a feature vector, based on intensity. The values maxF and minF are, respectively, the maximum and the minimum values of F for the all the pixels in the image. When the difference between two feature vectors of two pixels i and j is small, the similarity is close to 1. Otherwise, it has a value close to 0. For two pixels separated by a large distance (greater than r) the similarity is very close to zero, so it is set to zero explicitly for computational efficiency. Finally, the similarity matrix S is normalized by S ⫽ S/max(S) as proposed by Meila and Shi (2000). Alternative ways to define the similarity matrix based on machine learning techniques are presented by Bach and Jordan (2004). A simple numerical example is provided for illustration. Consider an image I (see Fig. 2a), consisting of 3 ⫻ 3 pixels, and represented as the following matrix:
冢 冣 2 1 1
I⫽ 2 1 3 . 3 3 3 Then, the similarity matrix between the pixels, will be a 9 ⫻ 9 matrix with the following values:
S⫽
冢
0
0.7 0.3
1
0.7
0
1
0.5
1
0
0
0
0
1
0
0.3
1
0
0
0
0
0.5 0.3
0
冣
0.5 0.3 0.3 0.3 0.2
0.3 1
Fig. 2. Graph of a sample image. (a) Initial image with three groups of pixels, each having the intensities value of 1, 2, 3. (b) The graph based on the image. The nodes correspond to the pixels, and the weighted edges are the similarities between pixels. The optimal partition is shown as dashed ellipses that provide a clustering of the image.
0.7 0.3 0.7 0.5 0.3
0.5
1
1
0.7
0
0
0
0
0
0.3
0
0
0.3
0
0
1
1
1
0.3
0
0
0.7
0
1
0
1
1
0.3
0
0
0.5
0
1
1
0
1
0.2
0
0
0.3
0
1
1
1
0
.
The graph resulting from the matrices I and S is presented in Fig. 2b. The optimal partition of the graph provides a solution to the image clustering problem. Recently, the Nyström method was proposed to efficiently approximate the eigen vectors (Fowkles et al.,
2004). For example, the method was recently used to segment the spine in magnetic resonance images (Carballido-Gamio et al., 2004). A requirement of supervised clustering is the number of clusters. Prior knowledge about the image content (i.e., anatomy in this case) can be used to decide the appropriate number of clusters. Examples of such approaches are presented in Kannan et al. (2004), Meila and Shi (2000), Ng et al. (2002). Choosing too few clusters will not partition the objects separately in the image. Choosing a very large number of clusters results in a single object represented as multiple clusters. A single object should be represented as a single cluster, but the correct number depends on the image content.
Ultrasound spectral clustering ● N. ARCHIP et al.
1489
Fig. 3. The Fiedler vector method applied to a noisy image having three clusters. (a) Initial image. (b) A graphic representation of the Fiedler vector illustrating the effect of noise. (c) An ideal Fiedler vector created by line fitting to make it easier to analyze. (d) Final segmentation.
One of the common problems of the standard NCut technique relates to its recursive nature. In other words, the image is initially bi-clustered and afterwards each of the two partitions are further clustered. It is difficult to come up with a robust strategy to decide in each step which cluster should be split further. Therefore, the clusters obtained with NCut tend to have the same size. Overcoming this problem is a challenging task. Unsupervised clustering, based solely on Fiedler eigenvector is a potential solution. In the unsupervised clustering, the user does not need to explicitly specify the number of clusters. As will be shown, an additional parameter is introduced to determine when there is a significant difference between clusters to allow automatic separation of clusters. More details follow.
Algorithm 1. Unsupervised clustering of a matrix Require: image I having the size N ⫻ N Ensure: the number of clusters K and the clusters C1, . . ., CK 1: Build the similarity matrix S, and the diagonal matrix of node degrees D 2: S ⫽ S/max(S) 3: Solve the generalized eigenvector problem (D ⫺ S)y ⫽ Dy 4: Consider 2 the second smallest eigenvalue and v2 its corresponding vector 5: Consider the permutation 㜷 ⫽ (i1, . ., iN2) sorting 2 6: Apply 㜷 to the image pixels (I1, . ., IN2). Obtain I⬘ ⫽ (I⬘1, . ., IN⬘ 2) 7: K ⫽ 1; CK ⫽ I⬘1 8: i ⫽ 2; 9: while i ⱕ N2 do 10: if I⬘i cannot be added to CK then 11: K ⫽ K ⫹ 1 12: end if 13: Add I⬘i to CK 14: i ⫽ i ⫹ 1 15: end while
1490
Ultrasound in Medicine and Biology
Algorithm 2. Segment the image I by windowing Require: Image I Ensure: the number of clusters K and the clusters C1, . . ., CK 1: Split the image I into the smaller windows W1, . . ., Wm 2: Apply the Algorithm 1 for W1 and obtain K and C1, . . ., CK 3: for i ⫽ 2 to m do 4: Apply the Algorithm 1 for Wi and obtain k⬘i and C⬘1, . . ., Ck⬘i⬘ 5: Merge clusters C⬘j and Cl if similar 6: Add clusters C⬘j to the list C1, . . ., CK if C⬘j not similar with any of existent clusters 7: end for
Unsupervised spectral clustering and the Fiedler eigen vector Instead of minimizing the NCut function, as suggested by Shi and Malik (2000), a property of the Fiedler vector is used. Let each pixel be Ij, j ⫽ 1, . . ., N2. Given the image pixels P ⫽ (I1, . ., IN2) and the Fiedler vector V ⫽ (1, . ., N2), we consider the permutation 㜷 ⫽ (i1, . ., iN2) that sorts the vector V: i1 ⱕ . . ⱕ iN2. By applying 㜷 to the pixel vector P, the vector (Ii1, . ., IiN2) is obtained. It satisfies the property that ?j1, . ., jk with @p, jp ⱕ l ⱕ jp ⫹ 1 : S(Il, Il⫺1) ⬍ ⑀, where S is the similarity function defined in eqn 5. This property provides a way to cluster the original image according to the given similarity metric. As an example, consider the matrix I previously defined. It contains three obvious clusters, given by the values 1, 2, 3. The similarity matrix is built and the Fiedler vector is obtained:
2 ⫽ 共0.10 0.41 0.44 0.03 0.39 ⫺0.34 ⫺0.28 ⫺0.34 ⫺0.36兲. Sorting the components of 2 gives a new vector:
2 ⫽ 共⫺0.36 ⫺0.34 ⫺0.34 ⫺0.28 0.03 0.10 0.39 0.41 0.44兲. To see the relevance of this order, the same permu-
Volume 31, Number 11, 2005
tation is applied to the image I, and expressed as a vector I⬘ which lists the image pixels row-by-row: I ⫽ 共3 3 3 3 2 2 1 1 1 兲 . In its new form, the vector components are grouped in compact blocks which represent the clusters. It becomes clear that the problem of clustering the initial matrix data is reduced to the problem of determining the sequences of maximum length having similar values in the vector obtained from the Fiedler eigen vector. Of course, for real images (usually affected by noise), the sequences are not perfectly aligned as shown in the simple example above. An example of a noisy image is presented in Fig. 3a. The noise in the image is apparent in the Fiedler vector shown in Fig. 3b. To accommodate the noise, a line fitting algorithm (Hair et al., 1998) is used to regularize the data. The result of the line fitting is illustrated in Fig. 3c. The three clusters from the regularization step can be used to segment the image, as shown in Fig. 3d. The pseudocode of the unsupervised NCut technique is listed in Algorithm 1. Calculating the eigen vectors for the Laplacian resulting from the whole image is computationally expensive. Therefore, a grid is used to divide the whole image into smaller windows. We then apply the clustering algorithm to each of these smaller windows. A significant improvement in the execution time is obtained, since the size of the matrices is reduced. Let W1, . ., Wm be the windows. First the window W1 is clustered, yielding K clusters. For each new window Wi clustered into k⬘i groups, if any of k⬘i can be combined with an existent cluster, then the two clusters are merged. Otherwise, the new cluster is added and the number of clusters K is updated. The pseudocode of this technique is listed in Algorithm 2. Since the unsupervised version of NCut is used, there is still no need to know or estimate K in advance. The unsupervised NCut method with windowing requires three parameters to be specified. The first pa-
Table 1. Simulated ultrasound images. The errors between the segmented object and the true object are calculated using the Hausdorff metric and the mean differences True object vs NCut
True object vs NCut-Unsupervised
Image
SNR
CNR
Hausdorff (mm)
Mean (mm)
Hausdorff (mm)
Mean (mm)
1 2 3 4 5 6 Average
4.40 4.36 4.27 4.19 3.90 2.20 3.88
1.59 1.54 1.46 1.40 1.20 0.75 1.32
3.5 2.7 3.7 3.6 3.5 5.1 3.68
1.4 1.9 1.2 1.3 2.2 3.3 1.88
3.3 2.2 3.2 3.8 4.0 4.9 3.56
1.2 1.3 1.4 1.5 2.0 3.1 1.76
Ultrasound spectral clustering ● N. ARCHIP et al.
1491
Fig. 5. Simulated ultrasound image with low-contrast and noisy. (a) Original image of 8 circular objects. (b) Segmentation using NCUT-Unsupervised. (c) The contour of the selected cluster corresponding to the circular objects.
information about the image content. However, for our experiments, we used the same value for all examples with good results. The second parameter (DifClusters) is the threshold value when two clusters are grouped together during windowing. This parameter is embedded in Algorithm 2 on line 5. It can also be customized for a particular application, but is kept constant in the images presented in this paper. The third parameter is the window size and is also kept constant. Tests Simulated ultrasound images are used first since the simulated object boundary is known. In this way, the supervised and unsupervised methods can be easily depicted and compared to the true segmentation. An ultra-
Fig. 4. Simulated ultrasound image. (a) Original image of a circular object. (b) The true contour of the object. (c) Segmentation using NCUT (5 clusters). (d) The contour of the selected cluster corresponding to the circular object. (e) Segmentation using NCUT-Unsupervised. (f) The contour of the selected cluster corresponding to the circular object.
rameter (DifPixels) specifies the maximum difference between two pixels in order to consider them as belonging to the same cluster. This parameter is embedded in Algorithm 1 on line 10. Naturally, this parameter could be determined for specific applications using a priori
Table 2. Liver cyst images. The errors between the two manually segmented contours of the cyst are calculated using the Hausdorff metric and the mean differences. The two manually segmented contours were drawn by two radiologists: Rad.1 and Rad.2
Image
Hausdorf error Rad.1 vs. Rad.2 (mm)
Mean error Rad.1 vs. Rad.2 (mm)
1 2 3 4 5 Average
2.8 2.9 2.4 3.1 2.0 2.6
1.8 1.7 1.4 1.6 1.1 1.5
1492
Ultrasound in Medicine and Biology
Volume 31, Number 11, 2005
Table 3. Liver cyst images. The errors between NCut and a manually segmented contour of the cyst are calculated using the Hausdorff metric and the mean differences. The manually segmented contours were drawn by radiologist Rad.1 or Rad.2 Image
Hausdorff error Rad.1 vs. NCut (mm)
Mean error Rad.1 vs. NCut (mm)
Hausdorff error Rad.2 vs. NCut (mm)
Mean error Rad.2 vs. NCut (mm)
1 2 3 4 5 Average
6.8 2.4 1.8 2.6 1.6 3.1
2.8 1.5 1.0 1.4 0.7 1.5
6.7 1.9 2.3 2.7 1.6 3.1
2.6 1.2 1.1 1.1 0.6 1.3
sound image of a circular object in a homogeneous background is created with the Field II program (Jensen, 1999). Field II models the acoustic field from an ultrasound transducer and simulates the interaction with the object by using the Tupholme-Stepanishen method for pulsed ultrasound. Six simulated images were created with various levels of Gaussian additive noise. Two different metrics are used to compute the error between the boundary of the segmented region and the known object. The Hausdorff metric (Hausdorff, 1962) is a common mathematical measure for comparing two sets of points in terms of their least similar members. Formally, given two finite point sets A ⫽ a1, . . ., ap and B ⫽ b1, . . ., bq, the Hausdorff metric is defined as: H(A, B) ⫽ max兵h(A, B), h(B, A)其 where h(A, B) ⫽ maxa⑀Aminb⑀B 储a ⫺ b储 and 储.储 is the Euclidean norm. Although theoretically attractive, the Hausdorff metric is very sensitive to noise in practice. Therefore, we also list the mean error representing the average difference in Euclidean distance. Abdominal images were also obtained for patients exhibiting cysts in the liver. The cysts have varying degrees of visibility, so are more difficult to segment. The cysts are generally benign, so there is no immediate clinical application for the segmentation; the results are meant to illustrate the operation of the clustering techniques. Five images were used.
The clustering techniques are also used to segment the amniotic fluid present in five images of the fetus. Again the main purpose is to illustrate the operation of the clustering techniques, but clinical applications include the diagnosis and quantification of hydramnios or oligohydramnios. Low contrast images of a liver hemangioma are also used. Finally, some sample images of bladder and artery segmentation are also shown. For the sets of liver and fetus images, a radiologist performed manual segmentation. The Hausdorff and mean errors were calculated between the radiologist and NCut, and the radiologist and NCut-Unsupervised. For the liver images, a second radiologist performed manual segmentation, so it at a comparison can be made between the results with each radiologist. For all tests, the NCut and NCut-Unsupervised were implemented in Matlab (The Mathworks Inc., Natick, MA, USA). No effort was made to optimize the program code for speed. It is expected that an implementation in a compiled computer language (such as C) would improve the performance by an order of magnitude. For NCut, the number of clusters was selected manually for each image to get the best segmentation. For NCut-Unsupervised, the same set of parameters was used. The values used in the tests are DifPixels ⫽ 5, DifClusters ⫽ 20 and window size ⫽ 15 pixels.
Table 4. Liver cyst images. The errors between NCut-Unsupervised and a manually segmented contour of the cyst are calculated using the Hausdorff metric and the mean differences. The manually segmented contours were drawn by radiologist Rad.1 or Rad.2 Image
Hausdorff error Rad.1 vs NCut-Uns. (mm)
Mean error Rad.1 vs. NCut-Uns. (mm)
Hausdorff error Rad.2 vs NCut-Uns. (mm)
Mean error Rad.2 vs NCut-Uns. (mm)
1 2 3 4 5 Average
6.5 2.3 1.7 2.5 1.2 2.8
2.6 1.2 0.8 1.2 0.5 1.3
6.0 2.0 2.5 2.4 1.6 2.9
2.4 1.1 1.1 1.0 0.7 1.3
Ultrasound spectral clustering ● N. ARCHIP et al.
1493
Table 5. Fetus images. The errors between NCut, NCutUnsupervised and a manually segmented contour of the amniotic fluid are calculated using the Hausdorff metric and the mean differences Manual vs. NCut
Manual vs. NCutUnsupervised
Image
Hausdorff
Mean
Hausdorff
Mean
1 2 3 4 5 Average
6.4 5.5 5.1 5.5 4.5 5.4
2.8 2.9 3.2 3.0 2.2 2.8
6.0 5.7 3.8 5.2 4.7 5.3
2.2 3.2 3.1 2.9 1.7 2.6
dowing, numerical problems arise because of the large number of computations required to cope with a matrix of N2 ⫻ N2. In the general case, for a window of m ⫻ m, and an image of n ⫻ n, the complexity of the algorithm is mn
Fig. 6. Liver cyst image. (a) Original image of a liver with a visible cyst. (b) The manually drawn contour of the cyst. (c) Segmentation using NCUT (9 clusters). (d) The contour of the selected cluster corresponding to the cyst. (e) Segmentation using NCUT-Unsupervised. (f) The contour of the selected cluster corresponding to the cyst.
RESULTS The computation time of the Matlab implementation of both NCut and NCut-Unsupervised (with windowing) is approximately 30 seconds on an image of 150 ⫻ 150. The most time-consuming portion of the algorithm is the eigenvector computation (based on the Lanczos method). As indicated previously, the speed of our new method is improved mainly by the windowing method. Without win-
Fig. 7. Fetus images. (a) An original image of a fetus. (b) The manually drawn edge of the amniotic fluid. (c) Segmentation using NCUT (4 clusters). (d) The manual grouping of the two clusters of amniotic fluid from NCut. (e) The contour of the selected clusters from NCut corresponding to the amniotic fluid. (e) Segmentation using NCut-Unsupervised. (f) The contour of the selected cluster from NCut-Unsupervised corresponding to the amniotic fluid.
1494
Ultrasound in Medicine and Biology
Volume 31, Number 11, 2005
Fig. 8. Low contrast image of hepatic hemangiomas. (a) Original image. (b) Segmented with unsupervised-NCUT. (c) The contour of the selected cluster of hepatic hemangiomas.
(O[np] ⫹ O(p[M{n}]), where O(O[np] ⫹ O(p[M{n}]) represents the complexity of the image, p is the maximum number of matrix-vector computations required and M(n) is the cost of computing the similarity matrix. The results in all cases is a partitioning of the image into a set of regions. The clinical user then can select the segmented region of interest and subsequently calculate geometric properties (dimensions, shape, area, volume), build anatomical models or perform other types of analysis. For this work, the boundaries of the segmented regions are compared. The results from the ultrasound simulations are presented in Table 1, and some example images are shown in Fig. 4. Low contrast images are shown in Fig. 5. Information about the image statistics, such as signal-to-noise ratio (SNR) and contrast-to-noise ration (CNR), is also presented. The results from the liver images are presented in Tables 2, 3 and 4. Comparisons are made with both NCut and NCut-Unsupervised using the manual segmentation
results as the reference. To understand the level of variability in the segmentation results, the errors between the results from the two radiologists were also computed and listed in Table 2. The contours of the segmentation results are shown in Fig. 6 for manual segmentation, NCut and NCut-Unsupervised. The results from the fetus images are presented in Table 5, and some example images are shown in Fig. 7. In some cases, the NCut algorithm represented the amniotic fluid as two or more clusters, so these were manually grouped together (simply by clicking with a mouse on both regions). No manual intervention was needed with NCutUnsupervised. The segmentation of low-contrast liver hemangioma is presented in the Fig. 8. The segmentation of the bladder is shown in Fig. 9, and the segmentation of the carotid artery is shown in Fig. 10. DISCUSSION The simulations show a good match between the spectral clusters and the true object. Small differences
Fig. 9. Bladder image. (a) Original image. (b) Segmentation with NCUT (7 clusters). (c) Segmentation with NCUTUnsupervised (7 clusters).
Ultrasound spectral clustering ● N. ARCHIP et al.
1495
Fig. 10. Carotic artery image. (a) Original image. (b) Segmentation with NCUT-Unsupervised. (c) The contour of the selected cluster of the artery.
along the boundary are likely due to the noise and speckle in the images. For the NCut algorithm, the average Hausdorff error is 3.68 mm and the average mean error is 1.88 mm. For the NCut-Unsupervised algorithm, the average Hausdorff error is 3.56 mm and the average mean error is 1.76 mm. There is no significant difference between the NCut and the NCut-Unsupervised techniques (at the 95% confidence level using the paired Student’s t-test). This gives confidence that the changes to the algorithm to make it unsupervised have not caused a drop in accuracy. For the tests on the liver cysts, there was also no significant difference between the following: (1) errors between NCut and manual segmentation and (2) the errors between the two manual segmentations (95% confidence level, using either the Hausdorff error or mean error). Nor was there a significant difference between the following: (1) errors between NCut-Unsupervised and manual segmentation and (2) the errors between the two manual segmentations. In other words, the results from NCut and NCut-unsupervised were within the variation of the two radiologists. For the tests on the amniotic fluid segmentation, there was also no significant difference between the NCut and the NCut-Unsupervised techniques when compared to the manual segmentation results. With NCutUnsupervised, the fluid and the fetus were clearly partitioned into separate clusters. On the other hand, as previously mentioned, the clusters obtained with standard NCut tend to have the same size, so the amniotic fluid was occasionally split in two clusters and manually combined by clicking both regions. The use of windows is transparent to the user. To show the insensitivity of the technique to the size of the window, the amniotic fluid segmented by NCut-Unsupervised with a window size 20 is shown in Fig. 11, compared to a window size of 15 in Fig. 7. There is little difference in the clustering results.
In these test images, the object of interest (simulated disk, liver cyst or amniotic fluid) is represented as a single connected region. It is worth noting that the object does not need to appear as a connected region to be considered the same object (same cluster). For example, the amniotic fluid in a 2D cross-sectional ultrasound image may appear as several unconnected regions surrounding the fetus. But these unconnected regions may still be partitioned as a single cluster. Therefore, it would be possible to select all amniotic fluid regions in a single step by selecting the single cluster representing those regions. Other segmentation techniques, such as active contours or seed-based segmentation, would require each region to be segmented separately before being joined together. All of the tests were done without changing the three parameters of the method. The standard NCut algorithm contains up to six parameters to adjust. The key parameter for the standard NCut is the number of clusters. The number was adjusted for best performance with these test, but the best number is nonintuitive so must be
Fig. 11. The amniotic fluid segmented with NCut-Unsupervised. The window size is 20. There is only small change compared to the use of a window of size 15.
1496
Ultrasound in Medicine and Biology
adjusted interactively. The simulated ultrasound images needed two clusters, the liver images needed 9, the fetus images required 4 and the bladder required 7 to segment the object of interest as a single region. Although good results were obtained in these tests, additional time and effort is required by the operator compared to NCutUnsupervised. The images of the partitioned regions also show the limitations of a spectral clustering approach. Clearly, not all separate tissue types are partitioned into separate regions. The objects of interest in this study were partitioned the easiest because they exhibited the most homogeneous properties. For example, the liver cyst forms a single cluster, but the rest of the liver is formed from four more clusters. Some of the liver clusters are also shared with non-liver tissue. Clearly, the liver could not be completely segmented with the current partition. So far, the method appears to work well for fluid-filled cavities and vessels. To perform other segmentation tasks, a more sophisticated similarity criterion may be needed. Additional features, such as texture, could be added to the similarity criterion. This may slow the creation of the similarity matrix, but the rest of the algorithm will proceed without change. Of course, by using a more specialized similarity criterion, the method may no longer work as a general purpose tool. Similarly, adding a priori information about the object may also help segmentation, but the trade-off between generality and specialization still exists. The key is that spectral clustering is a very active research area and many innovations from the nonmedical field may be incorporated here. The Gestalt basis of the approach provides a solid foundation from which to further the segmentation research. CONCLUSIONS AND FUTURE RESEARCH Spectral clustering techniques are shown to be suitable for ultrasound. The NCut technique has been adapted for use with ultrasound by incorporating an anisotropic diffusion preprocessing step, and by making the algorithm unsupervised. The key is to sort the Fiedler eigen vector. A windowing technique was also shown to be feasible as a method for improving the computational speed of the technique. The results on simulated data, liver cysts and amniotic fluid agree closely with results from manual segmentation. The main conclusion is that unsupervised spectral clustering is feasible. The spectral clustering technique is promising, and future research will focus on the development of new similarity criteria, the addition of other a priori information and further improvements to speed.
Volume 31, Number 11, 2005 Acknowledgments—This work is supported by the Natural Sciences and Engineering Research Council of Canada, the Institute of Robotics and Intelligent Systems, and by the NSF 0426558 and NIH/ NCRR P41 RR13218 grants.
REFERENCES Abolmaesumi P, Sirouspour MR. An interacting multiple model probabilistic data association filter for cavity boundary extraction from ultrasound images. IEEE Trans Med Imag 2004;23(6): 772–784. Akgul YS, Kambhamettu C, Stone M. A task-specific contour tracker for ultrasound. Proc IEEE Workshop on Mathematical Methods in Biomedical Image Analysis 2000;135–142. Archip N, Erard PJ, Egmont-Petersen M, Haefliger JM, Germond J. A knowledge-based approach for automatic detection of spinal cord in ct images. IEEE Trans Med Imag 2002;21(12):1504 –1516. Awad J, Abdel-Galil TK, Salama MMA, Fenster A, Rizkalla K, Downey DB. Prostate’s boundary detection in transrectal ultrasound images using scanning technique. IEEE Canadian Conf on Electrical and Computer Eng 2003:1199 –1202. Bach FR, Jordan MI. Learning spectral clustering. Advances in Neural Information Processing Systems 2004:305–312. Botros MN. A pc-based tissue classification system using artificial neural networks. IEEE Trans Instr Measure 1992;41(5):633– 638. Carballido-Gamio J, Belongie SJ, Majumdar S. Normalized cuts in 3-d for spinal mri segmentation. IEEE Trans Med Imag 2004;23(1): 36 – 44. Choy MM, Jin JS. Morphological image analysis of left-ventricular endocardial borders in 2d echocardiograms. SPIE Proc Medical Imaging 1996:852– 863. Chung FRK. Spectral Graph Theory. Providence, RI: American Mathematical Society, 1997. Evans AN, Nixon MS. Biased motion-adaptive temporal filtering for speckle reduction in echocardiography. IEEE Trans Med Imag 1996;15(1):39 –50. Fowkles C, Belogie S, Chung F, Malik J. Spectral grouping using the nystrom method. IEEE Trans Patt Anal Mach Intel 2004;26(2): 214 –225. Ghanei A, Soltanian-Zadeh H. A discrete curvature-based deformable surface model with application to segmentation of volumetric images. IEEE Trans Info Tech Biomed 2002;6(4):285–295. Ghiachetti A. On line analysis of echocardiographic image sequences. Med Image Anal 1998;2(3):1–25. Hair JF, Tatham RL, Anderson RE, Black W. Multivariate Data Analysis Prentice Hall; 5th edition, 1998. Hamarneh G, Gustavsson T. Combining snakes and active shape models for segmenting the human left ventricle in echocardiographic images. Comput Cardiol 2000;27:115–118. Hausdorff F. Set Theory. New York: Chelsea, 2nd edition, 1962. Hu N, Downey DB, Fenster A, Ladak HM. Prostate surface segmentation from 3d ultrasound images. Proc IEEE Intl Symp Biomed Imag 2002;613– 616. Jacob M, Blu T, Unser M. A unifying approach and interface for spline-based snakes. Proc SPIE Intl Symp Med Imag 2001:340 – 347. Jensen JA. A new calculation procedure for spatial impulse responses in ultrasound. J Acous Soc Amer 1999;105:3266 –3274. Kannan R, Vempala S, Vetta A. On clusterings: Good, bad and spectral. Journal of the Association for Computing Machinery 2004; 51(3):497–515. Kass M, Witkin A, Terzopoulos D. Snakes: active contour models. Intl J Comput Vis 1988;1(4):321–331. Kaus MR, Warfield SK, Nabavi A, Black P, Jolesz FA, Kikinis R. Automated segmentation of mr images of brain tumors. Radiology 2001;218(2):586 –591. Ladak HM, Mao F, Wang Y, Downey DB, Steinman DA, Fenster A. Prostate boundary segmentation from 2D ultrasound images. Med Phys 2000;27:1777–1788.
Ultrasound spectral clustering ● N. ARCHIP et al. Liu J, Tang YY. Adaptive image segmentation with distributed behavior-based agents. IEEE Trans Patt Anal Mach Intell 1999;10(6): 544 –551. Meila M, Shi J. Learning segmentation by random walks. Neural Information Processing Systems 2000:873– 879. Montagnat J, Sermesant M, Delingette H, Malandain G, Ayache N. Anisotropic filtering for model-based segmentation of 4d cylindrical echocardiographic images. Patt Recog Lett 2003;24:815– 828. Mudabhushi A, Metaxas DN. Combining low-, high-level and empirical domain knowledge for automated segmentation of ultrasonic breast lesions. IEEE Trans Med Imag 2003;22(2):155–169. Muzzolini R, Yang YH, Pierson R. Multiresolution texture segmentation with application to diagnostic ultrasound images. IEEE Trans Med Imag 1993;12(1):108 –123. Ng AY, Jordan M, Weiss Y. On spectral clustering: analysis and an algorithm. Advances in Neural Information Processing Systems 2002:849 – 856. Pathak D, Grimm PD, Chalana V, Kim Y. Pubic arch detection in transrectal ultrasound guided prostate cancer therapy. IEEE Trans Med Imag 1998;17(5):762–771.
1497
Pathak D, Haynor DR, Kim Y. Edge-guided boundary delineation in prostate ultrasound images. IEEE Trans Med Imag 2000;19(12): 1211–1219. Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Trans Patt Anal Mach Intell 1990;12(7):629 – 639. Perona P, Shiota T, Malik J. Anisotropic diffusion. In: Romeny BM, ed. Geometry-Driven Diffusion in Computer Vision 1, Norwell, MA, Kluwer Academic Publishers, 1994:73–92. Shen D, Zhan Y, Davatzikos C. Segmentation of the prostate boundaries from ultrasound images using statistical shape model. IEEE Trans Med Imag 2003;22(4):539 –551. Shi J, Malik J. Normalized cuts and image segmentation. IEEE Trans Patt Anal Machine Intell 2000;22(8):888 –905. Strintzis MG, Kokkinidis I. Maximum likelihood motion estimation in ultrasound image sequences. IEEE Signal Proc Lett 1997;4(6): 156 –157. Yu SX, Shi J. Segmentation given partial grouping constraints. IEEE Trans Patt Anal Mach Intell 2004;26(2):173–183.