IEEE TRANSACTIONS O N SIGNAL PROCESSING
VOL
90 I
10 NO 1 APKll 1992
An Adaptive Clustering Algorithm for Image Segmentation Thrasyvoulos N . Pappas
Abstract-The problem of segmenting images of objects with smooth surfaces is considered. The algorithm we present is a generalization of the ,K-means clustering algorithm to include spatial constraints and to account for local intensity variations in the image. Spatial constraints a r e included by the use of a Gibbs random field model. Local intensity variations are accounted for in a n iterative procedure involving averaging over a sliding window whose size decreases as the algorithm progresses. Results with an eight-neighbor Gibbs random field model applied to pictures of industrial objects, buildings, aerial photographs, optical characters, and faces, show that the algorithm performs better than the K-means algorithm and its nonadaptive extensions that incorporate spatial constraints by the use of Gibbs random fields. A hierarchical implementation is also presented and results in better performance and faster speed of execution. The segmented images a r e caricatures of the originals which preserve the most significant features, while removing unimportant details. They can be used in image recognition and as crude representations of the image. The caricatures are easy to display or print using a few grey levels and can be coded very efficiently. In particular, segmentation of faces results in binary sketches which preserve the main characteristics of the face, so that it is easily recognizable.
I. INTRODUCTION E present a technique for segmenting a grey-scale image (typically 2.56 levels) into regions of uniform or slowly varying intensity. The segmented image consists of very few levels (typically 2-4). each denoting a different region, as shown in Fig. 1. It is a sketch, or caricature, of the original image which preserves its most significant features, while removing unimportant details. It can thus be the first stage of an image recognition system. However, assuming that the segmented image retains the intelligibility of the original, it can also be used as a crude representation of the image. The caricature has the advantage that it is easy to display or print with very few grey levels. The number of levels is crucial for special display media like paper, cloth, and binary computer screens. Also, the caricature can be coded very efficiently, since we only have to code the transitions between a few grey levels. We develop an algorithm that separates the pixels in the image into clusters based on both their intensity and their Manuscript received April 5 , 1990; revised February 18. 1991. The author is with the Signal Processing Research Department. AT&T Bell Laboratories. Murray Hill. NJ 07974. IEEE Log Number 9106031.
Fig. I . Grey-scale image and its four-level segmentation
relative location. The intensity of each region is assumed to be a slowly varying function plus noise. Thus we assume that we have images of objects with smooth surfaces, no texture. We make use of the spatial information by assuming that the distribution of regions is described by a Gibbs random field. The parameters of the Gibbs random field model the size and the shape of the regions. A well-known clustering procedure is the K-means algorithm [ l ] , [2]. When applied to image segmentation this approach has two problems: it uses no spatial constraints and assumes that each cluster is characterized by a constant intensity. A number of algorithms that were recently introduced by Geman and Geman [ 3 ] , Besag [4], Derin and Elliot [SI, Lakshmanan and Derin [6], and others, use Gibbs random fields and can be regarded as an extension of this algorithm to include spatial constraints. They all assume, however, that the intensity (or some other parameter) of each region is constant. The technique we develop can be regarded as a generalization of the K-means clustering algorithm in two respects: it is adaptive and includes spatial constraints. Like the K-means clustering algorithm, our algorithm is iterative. Each region is characterized by a slowly varying intensity function. Given this intensity function, we define the a posteriori probability density function for the distribution of regions given the observed image. The density has two components, one constrains the region intensity to be close to the data and the other imposes spatial continuity. The algorithm alternates between maximizing the a posteriori probability density and estimating the intensity functions. Initially, the intensity functions are constant in each region and equal to the K-means cluster centers. As the algorithm progresses, the intensities are updated by averaging over a sliding window whose size
1053-587X/92$03 00
C
1992 IEEE
I
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
902
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 40. NO. 4, APRIL 1992
progressively decreases. Thus the algorithm starts with global estimates and slowly adapts to the local characteristics of each region. Our algorithm produces excellent caricatures of a variety of images. They include pictures of industrial objects, buildings, optical characters, faces, and aerial photographs. The binary sketches of faces, in particular, preserve the main characteristics of the face, so that it is easily recognizable. They can thus be used both as representations of the face and as input to a face recognition system. The performance of our technique is clearly superior to K-means clustering for the class of images we are considering, as we show with several examples. It is also better than the nonadaptive extensions of K-means clustering that incorporate spatial constraints by the use of Gibbs random fields. Our experimental results indicate that the performance of the algorithm is also superior to existing region growing techniques (see [7] for a survey). In the two-level segmentation case our technique compares favorably to other adaptive thresholding techniques that can be found in [81, P I . We also compare our adaptive clustering algorithm with the edge detector of [lo], [ 111. Its advantage over the edge detection approach is that it works with regions. In the edge detection case it is often difficult to associate regions with the detected boundaries. In our case regions are always well defined. The comparison becomes particularly interesting when we consider the case of faces. We will show that it is usually easier to recognize a face from a bilevel caricature than from its edges. We conclude our presentation of the adaptive clustering algorithm by considering a hierarchical multigrid implementation. We construct a pyramid of different resolutions of an image by successively filtering and decimating by two, starting from the highest resolution image. The development relies on understanding the behavior of the algorithm at different levels of the resolution pyramid. The hierarchical implementation results in better performance and faster speed of execution. The model for the class of images we are considering is presented in Section 11. The algorithm is presented in Section 111. In Section IV we examine the behavior of the algorithm on a variety of images. In Section V we consider a hierarchical implementation of the algorithm. The paper's conclusions are summarized in Section VI. 11. MODEL We model the grey-scale image as a collection of regions of uniform or slowly varying intensity. The only sharp transitions in grey level occur at the region boundaries. Let the observed grey scale image be y. The intensity of a pixel at location s is denoted by y5, which typically takes values between 0 and 255. A segmentation of the image into regions will be denoted by x , where x , = i means that the pixel at s belongs to region i. The number of different region types (or classes) is K . In this section we develop a model for the a posteriori probability den-
sity function p ( x 1 y). By Bayes' theorem P(X
I Y)
a P( Y I x)p(x>
(1)
where p ( x ) is the a priori density of the region process and p ( y I x ) is the conditional density of the observed image given the distribution of regions. We model the region process by a Markov random field. That is, if N,y is a neighborhood of the pixel at s, then (2) p(x, I xq, all q + s) = P(X, I xq, q E N,). We consider images defined on the Cartesian grid and a neighborhood consisting of the 8 nearest pixels. According to the Hammersley-Clifford theorem [ 121 the density of x is given by a Gibbs density which has the following form: (3) Here 2 is a normalizing constant and the summation is over all cliques C . A clique is a set of points that are neighbors of each other. The clique potentials V , depend only on the pixels that belong to clique C . Our model assumes that the only nonzero potentials are those that correspond to the one- and two-point cliques shown in Fig. 2. The two-point clique potentials are defined as follows: (-B. if x.. = x.. and s. a E C
The parameter /3 is positive, so that two neighboring pixels are more likely to belong to the same class than to different classes. Increasing the value of /3 has the effect of increasing the size of the regions and smoothing their boundaries. The one-point clique potentials are defined as follows:
Vc(x)
= ai,
if x,
=
i and s E C , all i.
(5)
The lower the value of a ' , the more likely that pixel s belongs to class i. Thus, the parameters aireflect our a priori knowledge of the relative likelihood of the different region types. In the remainder of this paper we will assume that all region types are equally likely and set a' = 0 for all i. A more extensive discussion of Gibbs random fields can be found in [3]-[5]. The conditional density is modeled as a white Gaussian process, with mean pi and variance oz.' Each region i is characterized by a different pi which is a slowly varying function of s. Thus, the intensity of each region is modeled as a signal pi plus white Gaussian noise with variance
d. The combined probability density has the form
'Note that Ihe superscript i in pi and a' is an index, not an exponent. The few instances in the text where an exponent is used, as in 0 2 , should be obvious from the context.
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
903
PAPPAS: ADAPTIVE CLUSTERING ALGORITHM
Fig. 2 . Cliques for Gibbs probability density function
We observe that the probability density function has two components. One constrains the region intensity to be close to the data. The other imposes spatial continuity. Note that if the intensity functions pi do not depend on s, then we get a model similar to those used in 131-[5]. If we also set p = 0, then we get the K-means algorithm. Thus our technique can be regarded as a generalization of the K-means clustering algorithm. The Gibbs random field introduces the spatial constraints, and the pixel dependence of the intensity functions makes it adaptive. We found both of these additions to be essential for the proposed algorithm. If we assume that the parameter and the noise variance 2 are known, then we must estimate both the distribution of regions x and the intensity functions p:. This will be done in the following section. Note that the only constraint we impose on the functions p: is that they are slowly varying. The estimation procedure we present in the following section uses this fact and implicitly imposes the constraint. The problem with using a parametric model (e.g., polynomial) for the region intensities is that, since the regions and intensity functions must be estimated concurrently, the parametric model becomes sensitive to inaccurate and fragmented region estimates. Alternative models for the conditional density are possible. For example, Chalmond [13] and Dinten [14] use Besag's autonormal model [4] in the context of image restoration, and Bouman and Liu [15], [16] use an autoregressive model for segmentation of textured images. We would like to point out that the Gibbs field is a good model of the region process only if we have a good model for the conditional density. The Gibbs model by itself is not very useful. It is easy to see that, in the absence of an observation, the optimal segmentation is just one region. In their discussion of the Ising model, Kindermann and h e l l [17] observe that the Gibbs field is ill behaved if there is no strong outside force (observed image with edges, in our case). 111. ALGORITHM In this section we consider an algorithm for estimating the distribution of regions x and the intensity functions p:. Note that the functions pi are defined on the same grid as the original grey-scale image y and the region distribution x. Like the K-means clustering algorithm, our algorithm is iterative. It alternates between estimating x and the intensity functions p:. First, we consider the problem of estimating the local intensity functions. Given the region labels x,we estimate the intensity pi at each pixel s as follows. We average the
I "
..._._.__. .'
Fig. 3 . Local average estimation.
grey levels of all the pixels that belong to region i and are inside a window of width W centered at pixel s, as shown in Fig. 3 . When the number of pixels of type i inside the window centered at s is too small, the estimate pi is not very reliable. At these points the intensity pi is undefined. When this happens, x, cannot be assigned to level i . Thus we have to pick the minimum number of pixels T,,, that are necessary for estimating p:. The higher this threshold the more robust the computation of pf. On the other hand, any small isolated regions with area less than this threshold may disappear. Thus the presence of this threshold strengthens the spatial constraints. A reasonable choice for the minimum number of pixels necessary for estimating p: is Tmln= W , the window width. This value of the threshold guarantees that long one-pixel-wide regions will be preserved. We have to obtain estimates of pi for all region types i and all pixels s. Clearly, the required computation is enormous, especially for large window sizes. Instead, we compute the estimates p: only on a grid of points and use bilinear interpolation to obtain the remaining values. The spacing of the grid points depends on the window size. We chose the spacing equal to half the window size in each dimension (50% overlap). Fig. 4 shows a window and the location of the centers of adjacent windows. The functions pi are smooth and therefore the interpolated values are good approximations of the values we would obtain by an exact calculation. Note that the larger the window the smoother the functions, hence the spacing of the grid points increases with the window size. Note also that, because of the dependence of the grid spacing on the window size, the amount of computation is independent of window size. Fig. 5 shows an original image, its two-level segmentation, and the local intensity functions for the two levels. Second, we consider the problem of estimating the distribution of regions. Given the intensity functions PI,, we
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 40. NO. 4. APRIL 1992
904
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Fig. 4. Window locations for local average estimation.
Fig. 5. Local intensity functions. (a) Original image. (b) Two-level segmentation. (c) Level 0 intensity function. (d) Level 1 intensity function.
want to maximize the a posteriori probability density (6) to obtain the MAP estimate of x. Finding the global maximum of this function requires a lot of computations. It
can be done using the Gibbs sampler and simulated annealing [3]. Instead, we maximize the conditional density at each point x, given the data y and the current x at all
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
905
PAPPAS: ADAPTIVE CLUSTERING ALGORITHM
1
other points:
get segmentation i by K-means
I
window equals whole image
The equality on the left follows from the Markovian property and the whiteness of the noise. The maximization is done at every point in the image and the cycle is repeated until convergence. This is the iterated conditional modes (ICM) approach proposed by Besag [4], sometimes also called the greedy algorithm [ 181. It corresponds to instantaneous freezing in simulated annealing. Bouman and Liu use this greedy algorithm in their multiple resolution approach for image segmentation [ 151, [ 161. The order in which the points are updated within a cycle depends on the computing environment [4]. In our implementation the updating was done in a raster scan. This procedure converges to a local maximum [4]. It stops when the number of pixels that change during a cycle is less than a threshold (typically M / 10 for an M X M image). The convergence is very rapid, usually in less than five cycles, consistent with Besag’s observations. Now we consider the overall adaptive clustering algorithm. We obtain an initial estimate of x by the K-means algorithm [2]. Starting from this segmentation, our algorithm alternates between estimating x and estimating the intensity functions pi. We define an iteration to consist of one update of x and one update of the functions pi. The window size for the intensity function estimation is kept constant until this procedure converges, usually in less than ten iterations (The number of iterations could vary with image content.) Our stopping criterion is that the last iteration converges in one cycle. The whole procedure is then repeated with a new window size. The adaptation is achieved by varying the window size W. Initially, the window for estimating the intensity functions is the whole image and thus the intensity functions of each region are constant. As the algorithm progresses, the window size decreases. The reason for this is that in the early stages of the algorithm, the segmentation is crude, and a large window is necessary for robust estimation of the intensity functions. As the algorithm progresses, the segmentation becomes better, and smaller windows give more reliable and accurate estimates. Thus the algorithm, starting from global estimates, slowly adapts to the local characteristics of each region. A detailed flowchart of our algorithm is given in Fig. 6 . The algorithm stops when the minimum window size is reached. Typically we keep reducing the window size by a factor of two, until a minimum window size of W = 7 pixels. If the window size is kept constant and equal to the whole image, then the result is the same as in the approach of [3]-[6]. This nonadaptive clustering algorithm is only an intermediate result of our adaptive clustering algorithm. It results in a segmentation similar to the K-
get new estimates fis
I
given p: get new segmentation i
I
n=n+l
Fig. 6. Adaptive clustering algorithm
means; only the edges are smoother. The strength of the proposed approach lies in the variation of the window size which allows the intensity functions to adjust to the local characteristics of the image. Note that the estimation procedure imposes the smoothness constraint on the local intensity functions. Also, there is no need to split our regions into a set of connected pixels as is done in [ 191. The variable window size takes care of this problem. We assume that the noise variance 2 is known or can be estimated independently of the algorithm. This is a reasonable assumption in many cases, especially when we have images obtained under similar conditions. From (6) it follows that increasing 2 is equivalent to increasing 0. In fact, the parameter p of the Gibbs random field can be chosen so that the region model is weaker in the early stages (algorithm follows the data), and stronger in the later stages (algorithm follows model) [4], [20]. In our implementations we chose the constant value 0 = 0.5 for all the iterations. Once the algorithm converges at the smallest window size, then p could be increased to obtain smoother segment contours. Thus, given an image to be segmented, the most important parameter to choose is the noise variance. In Section V we show that the performance of the algorithm is reasonable over a wide range of noise variances. In fact, the noise variance controls the amount of detail detected by the algorithm. We also have to choose the number of different regions K. We found that K = 4 works well for most images. Note that, in contrast to the K-means algorithm and some of the other adaptive schemes [7], [8], the choice of K is not critical because the adaptation of the intensity functions allows the same region type to have different intensities in different parts of the image. In the following section we discuss the choice of K more extensively.
Ill
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
I
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 40. NO. 4. APRIL 1992
906
I (C)
Fig. 7 . Example A . (a) Original image. (b) K-means algorithm (K
The amount of computation for the estimation of x depends on the image dimensions, the number of different windows used, the number of classes and, of course, image content. As we pointed out earlier, the amount of computation for each local intensity calculation does not depend on the window size. We should also point out that both the local intensity calculation and each cycle of the probability maximization can be done in parallel. Thus, even though the algorithm is very slow on serial machines, it can be considerably faster on a parallel computer.
IV. EXAMPLES In this section we examine the performance of the algorithm on a variety of images. We discuss the choice of the number of regions K , and compare adaptive clustering to K-means clustering and edge detection. An original image is shown in Fig. 7(a); the resolution is 256 X 256 pixels and the grey-scale range is 0 to 255
=
4 levels). (c) Adaptive clustering ( K = 4 levels)
(8 b). White Gaussian noise with U = 7 grey levels has been added to this image. Fig. 7(b) shows the result of the K-means algorithm for K = 4 levels. The two problems we mentioned earlier are apparent. There are a lot of isolated pixels and areas of uniform intensity (e.g., walls on lower right) are not well separated. Fig. 7(c) shows the result of our adaptive algorithm for K = 4 levels. The superior performance of the adaptive algorithm is apparent. Both problems have been eliminated. We can understand the limitations of the K-means algorithm when we look at a scan line of the image. The dotted line in Fig. 8(a) shows a horizontal scan line in the bottom half of the image of Fig. 7(a), and the solid lines show the four characteristic levels selected by the K-means algorithm. The solid line of Fig. 8(b) shows the resulting segmentation. We observe that while p0 and p' are reasonable models of the local image intensity, p 2 and p3 are not. This is because the selection of the four levels is based on the histogram of the whole image and does not
n
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
907
PAPPAS: ADAPTIVE CLUSTERING ALGORITHM
255
255
image profile
image profile
P:
128
I'
1
64
64
0
0 0
64
128
192
0
255
64
255
128
192
255
(a)
(a)
0 I
255
image profile
image profile
.M3
v: P2
v' 64
.
0 1 0
.. .
I
I
,
64
128
192
0 255
(b) Fig 8 K-means algorithm (a) Image profile and characteristic levels (b) Image profile and model
necessarily match the local characteristics. Fig. 9(a) shows the intensity functions of the adaptive algorithm and Fig. 9(b) shows the corresponding segmentation. The intensity functions have indeed adapted to the local intensities of the image. Another original image is shown in Fig. 10(a); the resolution is 244 X 256 pixels and the grey-scale range is 0 to 255. The standard deviation of the noise was estimated to be (T = 7 grey levels. Fig. lO(b) shows the result of the K-means algorithm. Fig. 1O(c) shows the result of our adaptive algorithm. Again we chose K = 4. The problems we mentioned above are again evident in the K-means result. The adaptive algorithm gives a much better segmentation, even though some problems are still present (e.g., the rings of the wrench in the lower right are not well defined). Overall, we could say that the adaptive algorithm produces a simple representation of the image, while retaining most of the important information. One of the important parameters in the adaptive clustering algorithm is the number of classes K. In most clustering algorithms the choice of the number of classes is important and could have a significant effect on the qual-
0
64
128
192
5
(b)
Fig. 9. Adaptive clustering algorithm. (a) Image profile and local intensity functions ( I S X 15 window). (b) Image profile and model ( I S X 15 window).
ity of segmentation. Thus, many authors consider estimating of the number of classes as a part of their segmentation algorithms, for example [21], [16]. In the Case of the K-means algorithm, choosing the wrong number of classes can be disastrous as we will see in the following examples. However, our adaptive algorithm is quite robust to the choice of K. This is because the characteristic levels of each class adapt to the local characteristics of the image, and thus regions of entirely different intensities can belong to the same class, as long as they are separated in space. Fig. Il(a) shows the result of the K-means algorithm for K = 2. Because of the wide variation of the grey-scale intensity throughout the image, a lot of detail has been lost in the K-means case. Fig. ll(b) shows the result of our adaptive algorithm for K = 2. Here most of the detail in the image has been recovered.2 We can argue that the 'Note that the thin oblique lines of Figs. I I(b) and 7(b) have disappeared in Figs. 1 l ( b ) and 7(b), respectively. This is caused by the combination of the Gibbs constraints and the threshold T,,,, defined in Section 111. The only one-pixel-wide features preserved are those of high contrast.
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
I k k F I R A N S A C T I O N S ON S I G N A L PROCLSSING
9OX
VOL 40 N O 4 A P R I L I Y Y ?
ib)
Fig. I I . Example A. (a) K-means algorithm ( K = 2 levels). (b) Adaptive clustering ( K = 2 levels).
(C)
Fig. 10. Example B . ( a ) Original image. ib) K-nie;in\ algorithm ( K levels). (c) Adaptive clustering ( K = 4 l e v c l s ) .
=
4
segmentation is not perfect; indeed, some edges are missing (e.g., the sky in the upper right) and some spurious edges are present (e.g., around the windows in the lower left). Clearly, K = 4 results in a better segmentation.
In many images K = 4 is the best choice. This appears to be related to the four-color theorem [22], even though a proof would be very difficult. However, for K = 2 the adaptive algorithm results in a very good representation of the image, even though the individual segments don't always correspond to actual edges in the image. The value of K = 2 is especially important because it simplifies both the display (e.g., on paper or on binary computer screens) and coding of the image. As another example, consider the original image shown in Fig. 12(a); the resolution is 512 X 512 pixels and the grey-scale range is 0 to 255. The standard deviation of the noise was assumed to be (T = 4 grey levels. Fig. 12(b) shows the result of the K-means algorithm for K = 2. Fig. 12(c) shows the result of our adaptive algorithm. The advantage of the adaptive algorithm is apparent. For example, an airstrip at the lower left of the image has disappeared in the K-means case, but is nicely preserved by the adaptive algorithm.
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
PAPPAS: ADAPTIVE CLUSTERING ALGORII'HM
(c)
Fig. 12. Example C . (a) Original imagc. ( b ) K-means algorithm ( K
In some cases K is constrained by the statement of the problem. For example, consider the image of Fig. 13(a). Here K = 2 since each point of the image either belongs to the character or to the background. The resolution is 64 X 64 pixels and the grey-scale range is 0 to 255. This image was obtained by blurring the optical character "E:" with a Gaussian filter (the standard deviation is 2.4 pixels) and adding white Gaussian noise with U = 10 grey levels. Fig. 13(b) shows the result of the K-means algorithm. Fig. 13(c) shows the result of our adaptive algorithm. We observe that the adaptive algorithm recovers the hole in the center of the character. This could be crucial in the correct interpretation of the character. Thus our algorithm could be an important component of an optical character recognition system. The performance of the two-level adaptive algorithm is particularly good on grey-scale images of human faces. Consider the original image shown in Fig. 14(a); the resolution is 512 X 512 pixels and the grey-scale range is 0
m
=
2 levels). ( c ) Adaptive clu\tenng ( K = 2 level\)
to 255. The standard deviation of the noise was assumed to be U = 7 grey levels. Fig. 14(b) shows the result of the K-means algorithm and Fig. 14(c) shows the result of our adaptive algorithm. The K-means algorithm produces an image that is too light and, therefore, barely gets some of the facial characteristics. The adaptive algorithm, on the other hand, gives a very good sketch of the face. Note that the nose is implied by a couple of black regions at its base, much like an artist might sketch it with a few brush strokes. In general, the K-means algorithm produces a binary image that may be too light or too dark, depending on the overall image content. The proposed algorithm adapts to the local characteristics and preserves the features of the face. The binary sketches preserve the characteristics of the human faces as we have observed over several examples obtained under different conditions (see Section V for more examples). As we discussed, the sketches can be used for displaying or printing on special media like pa-
--
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
I
910
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 40. NO. 4, APRIL 1992
(a)
(C)
(b)
Fig. 13. Example D. (a) Original image. (b) K-means algorithm ( K
=
2 levels). (c) Adaptive clustering ( K
=
2 levels).
(d)
(C)
Fig. 14. Example E. (4 Original image. (b) K-means algorithm ( K = 2 levels). (c) Adaptive clustering ( K Edge detection.
II
=
2 levels). (d)
91 I
PAPPAS: ADAPTIVE CLUSTERING ALGORITHM
per, cloth and binary computer screens. They can be coded with very few bits of data, since we only have to code the transitions between a few grey levels. In [23] we show that coding of typical face caricatures requires less than 0.1 b/pixel, which is substantially lower than what waveform coding techniques require. The dual to image segmentation is edge detection. It is thus interesting to compare our segmentation results with an edge detector. Here we use the edge detector of Boie and Cox [lo], [ l l ] . Fig. 14(d) shows the result of the edge detector with integrated line detection; their algorithm includes an automatic threshold selection. The algorithm is very successful in detecting the most important edges. We observe, however, that it is not as easy to recognize the person from the edges as it is from the binary sketches. This demonstrates the advantage of representing an image as a set of black and white regions, as opposed to a set of edges. Note the basic tradeoff between region segmentation and edge detection. As we mentioned, not all region transitions correspond to edges in the image. On the other hand, the detected edges do not necessarily define closed regions. As an indication of the amount of computation required by the adaptive clustering algortihm, we look at the running time on a SUN SPARC station 1 +. The computation time for the four-level segmentation of Fig. 7(c) was 1 1 min of user time. The K-means clustering of Fig. 7(b) required only 10 s of user time. The computation time for the two-level segmentation of Fig. l l ( b ) was 7.6 min of user time and for the K-means clustering of Fig. 1 l(a) 4 sec of user time. The resolution of these images is 256 x 256 pixels. When the resolution increases to 512 x 512 pixels, as in Fig. 12, the computation times become 35 min for the adaptive algorithm and 13 s for the K-means algorithm. In our implementations we made no systematic attempt to make the code more efficient, and the choice of termination criteria was conservative.
where I, is the piecewise continuous signal and n, is white Gaussian noise with standard deviation U . The signal at the next resolution has the same form y; =
z:
+ n:
(9)
where I;, the filtered and decimated Z,, is approximately piecewise continuous and n:, the filtered and decimated n , , is approximately white Gaussian with standard deviation U’ = u/2. Thus, the model is still valid at the next resolution and only the noise standard deviation 0 changes by a factor of two. Based on our previous discussion of the role of the parameter p, it is reasonable to assume that it does not change with resolution. Bouman and Liu [15], 1161 adopt a similar approach in their multiple resolution segmentation of textured images. An interesting alternative approach for the choice of the Gibbs parameters is proposed by Gidas [24], who obtains an approximation of p at different resolutions based on a probabilistic model for the whole pyramid. Given our choice of and U , the only remaining parameter is the detection threshold Tmln.We assume that Tmln= W , the window width, is a reasonable choice for the detection threshold T,,,,at the highest level. As we lower the resolution we expect that significant features become smaller and, therefore, we should reduce the detection threshold T,,,,,. We decrease T,,, by a factor of two from one resolution to the next. The hierarchical implementation of the algorithm is described in Fig. 15. We start at the lowest resolution image, and apply our segmentation algorithm as described in the previous sections, starting from the K-means clusters. When the minimum window size is reached, we move to the next level in the pyramid and use the current segmentation, expanded by a factor of two, as a starting point. Since we now have a good starting point, we need not iterate through the whole sequence of windows. We start where we left off at the previous resolution, which is twice the minimum window size of that resolution. This approach reduces the amount of computation sigV . HIERARCHICAL MULTIGRID IMPLEMENTATION nificantly, since most iterations are at the lowest level of In this section we present a hierarchical implementation resolution. The improvement in performance is also sigof our segmentation algorithm. The goal is to improve nificant as the following example indicates. In our impleboth its performance and speed of execution. The develmentations we used a four-level pyramid. opment of such a hierarchical approach relies on underAn original image is shown in Fig. 16(a); the resolution standing the behavior of the algorithm at different levels is 512 x 512 pixels. We assume that the noise level of of the resolution pyramid. In particular, we would like to this image is U = 8 grey levels. Fig. 16(b) shows the know how the various parameters of the algorithm change result of the adaptive algorithm for K = 2 levels, applied with the resolution of the image. to the full-scale image. Fig. 16(c) shows the result of the We construct a pyramid of images at different resoluhierarchical implementation of the two-level adaptive altions in the following way. Starting from the highest resgorithm. Both result in very good segmentations. We can olution image, we obtain the next image in the pyramid argue, however, that the hierarchical segmentation has by filtering and decimating by a factor of two. The lowbetter adjusted to the local image characteristics, as can pass decimation filters we use have a sharp cutoff so that be seen by the difference in the lips and the area above the sharpness of the edges is preserved at the lower resthe hat. olutions. As an indication of the potential computational savings The image model we described in Section I1 can be of the hierarchical implementation we look at the running written as follows: times of the two algorithms on a SUN SPARC station 1 . YS = 4 + n, (8) The full-scale implementation of Fig. 16(b) required 36
+
111
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
1
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 40. NO. 4, APRIL 1992
912
lowest resolution image get segmentation i by K-means
window equals whole image
double window size window size
iterate to get new estimates
highest
Fig. 15. Hierarchical adaptlve clustering algorithm
min of user time, while the hierarchical implementation of Fig. 16(c) required only 8.2 min of user time. The first (lowest resolution) stage required 10 s of user time, the second 16 s, the third 118 s, and the last 348 s. Again, no systematic attempt was made in our implementations to improve the efficiency of the code. We now look at the performance of the algorithm for different noise-level estimates. An original image is shown in Fig. 17(a); the resolution is 512 x 512 pixels. Figs. 17(b)-(e), show the result of the hierarchical implementation of the adaptive algorithm (K = 2) for noise level estimates U = 4, 8, 16, and 32 grey levels, respectively. The original image is the same in all cases. Observe how the amount of detail detected by the algorithm decreases as the noise variance increases. This example illustrates that the algorithm’s performance is reasonable over a wide range of noise variances and, moreover, that we can control the amount of detail by choosing the noise variance. Thus, our knowledge of the level of noise determines the amount of detail in the segmented image.
VI. CONCLUSION We presented an algorithm for segmenting images of objects with smooth surfaces. We used Gibbs random fields to model the region process. The parameters of the Gibbs random field model the size and the shape of the regions. The intensity of each region is modeled as a slowly varying function plus white Gaussian noise. The adaptive technique we developed can be regarded as a generalization of the K-means clustering algorithm. We applied our algorithm to pictures of industrial objects, buildings, aerial photographs, faces, and Optical c h a m ters. The segmented images are useful for image display,
lhl
(C)
F~~ 16 Example F (a) Original image (b) Adaptive clustering ( K (c) Hierarchical adaptive clustering ( K = 2)
II
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
=
2)
913
PAPPAS: ADAPTIVE CLUSTERING ALGORITHM
(d) Fig. 17. Example G : Hierarchical adaptive clustering ( K
(e) =
2 ) . (a) Original image. (b)
U =
32.
U
=
4 . (c)
U
=
8 . (d)
U
=
16. ( e )
914
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 40, NO. 4, APRIL 1992
coding, and recognition. We believe that our approach can be extended to the Of images that regions with slowly varying textures, as well as moving images.
ACKNOWLEDGMENT The author like to thank I. ‘Ox for providing the code for the Boie-Cox edge detector.
REFERENCES [ I ] J. T . Tou and R. C . Gonzalez, Pattern Recognition Principles. Reading MA: Addison-Wesley, 1974. 121 R. M. Gray and Y. Linde, “Vector quantizers and predictive quantizers for Gauss-Markov sources,” IEEE Trans. Commun., vol. COM-30, no. 2, pp. 381-389, Feb. 1982. [3] S. Geman and D . Geman, “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images,” IEEE Trans. Putt. Anal. Machinelntell., vol. PAMI-6, no. 6, pp. 721-741, Nov. 1984. [4] J. Besag, “On the statistical analysis of ditty pictures,” J . Roy. Stat. Soc. B , vol. 48, no. 3, pp. 259-302, 1986. [5] H. Derin and H. Elliot, “Modeling and segmentation of noisy and textured images using Gibbs random fields,” IEEE Trans. Patt. Anal. Machine Intell., vol. PAMI-9, no. 1, pp. 39-55. Jan. 1987. [6] S . Lakshmanan and H . Derin, “Simultaneous parameter estimation and segmentation of Gibbs random fields using simulated annealing,” IEEE Trans. Patt. Anal. Machine Intell., vol. 1 I , no. 8, pp. 799813, Aug. 1989. [7] R . M. Haralick and L. G . Shapiro, “Image segmentation techniques,” Comput. Vision, Graph., Image Processing, vol. 29, pp. 100-132, 1985. [8] P. K. Sahoo, S . Soltani, and A . K . C. Wong, “A survey of thresholding techniques,” Comput. Vision, Graph., Image Processing, vol. 41, pp. 233-260, 1988. [9] J . S . Weszka, “A survey of threshold selection techniques,” Comput. Graph. Image Processing, vol. 7 , pp. 259-265, 1978. [lo] R. A. Boie and I. J. Cox, “Two-dimensional optimum edge recognition using matched and Wiener filters for machine vision,” in Proc. IEEE First Int. Con$ Compur. Vision (London, England), June 8- 1 I , 1987, pp. 450-456. [ I l l I . J. Cox, R. A . Boie. and D. A. Wallach. “Line recognitlon,” in Proc. Tenth Int. Cant Putt. Recogn. (Atlantic City, NJ), June 1621, 1990, pp. 639-645. [ 121 J . Besag, “Spatial interaction and the statistical analysis of lattice systems,’’ J . Roy. Stat. Soc. B , vol. 26, no. 2, pp. 192-236, 1974.
[I31 B. Chalmond, “Image restoration using an estimated Markov model,” Signal Processing, vol. 15, no. 2, pp. 115-129, Sept. 1988. [I41 J. M. Dinten, “Tomographic reconstruction of axially symmetric objects: Regularization by a Markovian modelization.” presented at the 10th 1nt.Congr. Patt. Recog., 1990. [I51 C . Bouman and B. Liu, “Segmentation of textured images using a multiple resolution approach,” in Proc. ICASSP-88 (New York), Apr. 1988,pp. 1124-1127. [16] C. Bouman and B. Liu, “Multiple resolution segmentation of textured images,” IEEE Trans. Putt. Anal. Machine Intell.. vol. 13, no. 2, pp. 991113, Feb. 1991. [I71 R. Kindermann and J. L. Snell, Markov Random Fields and rheir Applications. American Mathematical Society, 1980. [18] H. Derin, “Experimentations with simulated annealing in allocation optimization problems,’’ in Proc. CISS-86, pp. 165-171. [I91 T. Aach, U. Franke, and R. Mester, “Top-down image segmentation using object detection and contour relaxation,” in Proc. ICASSP-89 (Glasgow, U.K.), May 1989. pp. 1703-1706. [20] G. Wolberg and T. Pavlidis, “Restoration of binary images using stochastic relaxation with annealing,” Parr. Recogn. Lett., vol. 3. no. 6, pp. 375-387, Dec. 1985. [21] J. Zhang and J. W. Modestino, “A model-fitting approach to cluster validation with application to stochastic model-based image segmentation,” IEEE Trans. Pat:. Anal. Machines Intell., vol. 12, no. I O , pp. 1009-1017, Oct. 1990. [22] K. Appel and W . Haken, “Every planar map is four colorable: Parts I and 11,” 111. J . Math., vol. 21, no. 3, pp. 429-567, Sept. 1977. [23] T. N. Pappas, “Adaptive thresholding and sketch-coding of grey level images,” Proc. SPlE Int. Soc. Opt. Eng., vol. 1199, Nov. 1989. 1241 B. Gidas, “A renormalization group approach to image processing problems,” IEEE Trans. Putt. Anal. Machine Int., vol. 11, no. 2, pp. 164-180, Feb. 1989.
Thrasyvoulos N. Pappas received the S.B., S.M., and Ph.D. degrees in electrical engineering and computer science from the Massachusetts Institute of Technology, Cambridge, in 1979, 1982, and 1987, respectively. He is currently with the Signal Processing Research Department at AT&T Bell Laboratories, Murray Hill, NJ. His research interests are in image processing and computer vision.
Authorized licensed use limited to: Georgia State University. Downloaded on March 19, 2009 at 15:16 from IEEE Xplore. Restrictions apply.
