Gaussian Mixture based image segmentation algorithm

Pergamon Pattern Recognition, Vol. 31, No. 3, pp. 315 325, 1998 © 1997 Pattern Recognition Society. Published by Elsevier Science Ltd Printed in Great Britain. All rights reserved 0031-3203/98517.00 + .00 PII: S0031-3203 (97) 00045-9 A GAUSSIAN-MIXTURE-BASED IMAGE SEGMENTATION ALGORITHM LALIT GUPTA and THOTSAPON SORTRAKUL Department of Electrical Engineering, Southern Illinois University, Carbondale, IL 62901, U.S.A. (Received 12 October 1995; accepted 24 May 1996) Abstract--This paper focuses on the formulation, development, and evaluation of an autonomous segmen- tation algorithm which can segment targets in a wide class of highly degraded images. A segmentation algorithm based on a Gaussian-mixture model of a two-class image is selected because it has the potential for effective segmentation provided that the histogram of the image approximates a Gaussian mixture and the parameters of the model can be estimated accurately. A selective sampling approach based on the Laplacian of the image is developed to transform the histogram of any image into an approximation of a Gaussian mixture and a new estimation method which uses information derived from the tails of the mixture density is formulated to estimate the model parameters. The resulting selective-sampling-Gaus- sian-mixture parameter-estimation segmentation algorithm is tested and evaluated on a set of real degraded target images and the results show that the algorithm is able to accurately segment diverse images. © 1997 Pattern Recognition Society. Published by Elsevier Science Ltd. Segmentation Histogram Thresholding Gaussian mixture 1. INTRODUCTION This paper focuses on developing an autonomous segmentation algorithm which can segment targets in a wide class of degraded images. The autonomous segmentation of targets is one of the most complex problems in automatic target recognition (ATR) as the environment of the target is not controllable. Segmentation in the context of ATR problems is the process of separating a target from the background as accurately as possible. Because the subsequent stages of the ATR system depend upon the information provided by the segmentation stage, reducing segmen- tation error clearly reduces the probability of target misclassification. A multitude of image segmentation algorithms have been developed based on grey level thresholding, 11-13~ edge detection, ~11-16) region grow- ing, m'12'tT-19) region splitting and merging,/11,12,2°-22) re- laxation,~11,12,23 25~ and fuzzy set theory313'26-28) It is well known in image segmentation research that there is no single algorithm that can effectively segment all types of images. Additionally, the application of differ- ent algorithms to the same image generally produce different segmentation results. The selection of an appropriate segmentation algorithm depends largely on the type of images and the application areas. One approach is to select/design an effective algorithm for a class of images generated by a particular sensor. This approach is especially suitable for segmenting images generated in problems such as industrial in- spection where the environment is controllable. The drawback to this approach when applied to segment- ing images in problems such as ATR is that there is no control of the environment. The images, therefore, have varying characteristics and the results of seg- menting some images may be seriously incorrect. Since the performance of the segmentation stage has a major influence on the performance of the ATR system, the goal is to focus on the formulation of an autonomous segmentation algorithm capable of accu- rately segmenting a wide class of highly degraded images. The segmentation algorithm developed in this pa- per is based on a parametric model in which the probability density function of the gray levels in the image is a mixture of two Gaussian density functions. This model has received considerable attention in the development of segmentation algorithms and it has been noted that the performance is influenced by the shape of the image histogram and the accuracy of the estimates of the model parameters. The model-based segmentation algorithms perform poorly if the histo- gram of an image is a poor approximation of a mix- ture of two Gaussian density functions. It is especially poor when the two modes of the mixture density are not well-separated. The application of this model in image segmentation is, therefore, highly limited as the histograms of the images encountered in practical problems are seldom approximations of Gaussian mixtures with well-defined modes. Additionally, no method has been proposed to broaden the application of Gaussian-mixture-based segmentation to images that have histograms which are poor approxima- tions of Gaussian mixtures. In order to broaden the application, a selective sampling method is pro- posed to transform the histogram of any image into 315 316 L. GUPTA and T. SORTRAKUL a histogram which approximates a Gaussian mixture. A new method is also proposed to accurately estimate the model parameters by first estimating the two Gaussian components of the mixture density and then estimating the parameters of the model from the two components. 2. THE GAUSSIAN MIXTURE MODEL FOR IMAGE SEGMENTATION Gaussian-mixture-based segmentation is based on histogram thresholding which is one of the most fre- quently used methods for segmenting images. In his- togram thresholding, it is assumed that an image consists of two regions or classes: target and back- ground, with each region having a uni-modal gray-level distribution. The segmentation problem, therefore, involves determining a suitable threshold to partition the image into the target and background regions. Gaussian-mixture segmentation algorithms assume a parametric model in which the probability density function (PDF) of the gray levels in the image are a mixture of two Gaussian density functions with given means, standard deviations, and proportions. That is, it is assumed that an image with the target and background intensities is combined with Gaus- sian noise. The mixture probability density function of such an image can be written as ~t21 p(z) = (P1)p(z/C O + (Pz)p(z/C2), where p(z/Ci) = (2~0./2) - 1/2exp [ - (z - pi)2/(20.2)], mixture component i, i = 1, 2 pi = mean of intensity level i, i = 1, 2 0.i = standard deviation of intensity level i, i = 1, 2 Pi = prior probability of intensity level i, i = 1, 2; (P1 + P2 = 1). For a given threshold T, the overall probability of error is E(T) = PxEI(T) + PIE2(T), where and Ea(T)=f~: ,p (z /Cz )dz E2(T) = ]~jr p(z/C1)dz. The threshold value for which E(T) is minimum is given by solving AT 2 + BT + C=O, (1) where A = o , 2 - 0 .~, B = 2(pla 2 - p2a~), 20.2 ,20.2 C = t~,2 - - IA1 2 q- 2 r r2a 2 ln(0.2P1/alP2). If the five model parameters P1 (or P2),/~1,1~2, 0.1, and 0"2 are known, the optimum threshold can be deter- mined arlalytically from equation (1), however, if the parameters are unknown, they have to be estimated from the gray-level histogram of the image. Several image segmentation algorithms based on estimating the parameters of the Gaussian mix- ture using curve-fitting techniques, (1'3) truncated histograms,16, 71 and the E-M algorithm I1°1 have been proposed. The algorithms based on curve-fitting techniques use a goodness-of-fit criterion through a conjugate gradient hill-climbing procedure to fit histograms with bimodal Gaussian density curves. The drawback with the curve-fitting methods is that they are computationally involved. The drawback with the iterative methods which estimate the model parameters from truncated histograms is that the esti- mation of parameters is biased even if the optimal threshold is selected. An improvement in estimating the variances when a threshold in an iteration is away from the overlapping regions of the mixture compo- nents is proposed in reference (7), however, the para- meters are otherwise estimated as described in refer- ence (6). The drawback of the E-M-based algorithm is that it tends to give poor estimates for small sample sizes and the time for convergence can be extremely high. Additionally, the E M algorithm fails to con- verge if one or both variances approach zero as in the cases of targets and backgrounds with uniform intensities. The performance of Gaussian-mixture-based seg- mentation algorithms clearly depend on how well the histogram of an image approximates a Gaussian mix- ture and the accuracy of the estimates of the model parameters. The histograms of real images seldom approximate a Gaussian mixture and none of the algorithms described propose methods to broaden the application of Gaussian-mixture-based segmentation to images with histograms that are a poor approxima- tion of a Gaussian mixture. The focus of the segmen- tation algorithm developed in this paper, therefore, is to broaden the application of Gaussian mixture segmentation algorithms through a histogram transformation approach and to accurately estimate the mixture density parameters through a method specifically designed to estimate the mixture density parameters from the transformed histogram. The formulations of the selective sampling algorithm for histogram transformation, the Gaussian-mixture parameter-estimation algorithm, and the selective- sampling-Gaussian-mixture parameter-estimation image segmentation algorithm are described in the following sections. Gaussian-mixture-based image 317 3. THE SELECTIVE SAMPLING (SS) ALGORITHM Several histogram modification methods for image segmentation have been proposed to facilitate the detection of the bottom of the valley between two modes by producing a transformed gray-level histo- gram which is more bimodal, or has a deeper valley between the modes, or converting the valley between the modes into a peak. (5' 12) The approach developed in this paper focuses on selecting a set of pixels which are typical of the target and the background. Only the selected pixels contribute to the transformed histo- gram. The pixels selected are the set of pixels which most likely fall on the two sides of the boundary of a target. The Laplacian L[f(x, y)] of the imagef(x, y) is used to identify and select such pixels. 12 The La- placian is a second-derivative operator which gives high values to pixels on both sides of the edge (modelled as a step or a ramp function) and low- values to all other pixels. Pixels on the bright side of an edge have a negative Laplacian whereas pixels on the dark side of an edge have a positive Laplacian. The histogram of the selected non-zero Laplacian pixels will, therefore, be bi-modal with the 2 modes representing typical target and background gray- levels and the valley between the two peaks would also be well-defined. Additionally, the two peaks will be approximately equal in size because the number of pixels on both sides of the boundary would be ap- proximately equal even if the target and the back- ground have unequal areas. The Laplacian of a continuous function g(x, y) is defined as L[g(x, y)] = [(02g)/Ox 2) q-(a2g)/c~y2)]. For a discrete imagef(x, y), the Laplacian at (x, y) can be approximated by L[f(x, y)] =f(x - 1,y) +f(x + 1, y) +f(x, y - I) +f(x, y + 1) - 4[f(x, y)]. The Laplacian appears to be ideal for selecting pixels to generate a well-defined bi-modal histogram as it produces values with differing signs at each edge, thus, facilitating the identification of target and back- ground pixels. Rules based on the sign and magnitude of the Laplacian can be easily developed to identify the target and background pixels whether an edge is modeled as a step or a ramp. The problem with the Laplacian is that as a sec- ond-order derivative, it is highly sensitive to noise and produces high values not only at the boundary of the target but also in any region where sharp gray-level transitions occur. Therefore, the applicability of this method to highly degraded images is questionable. To overcome this problem, the pixel selection rule de- veloped is based on using the Laplacian to accurately locate edges in noise and using the connectivity of the pixels to exclude pixels that may be noise pixels. Without a loss of generality, it will be assumed that an image is composed of a bright target on a dark back- ground. In the first step, the sign transition in the Laplacian is used to determine which type of edge is encountered when the image rows and columns are scanned. A transition from plus to minus indicates a step edge transition from a background pixel to a target pixel whereas a transition from minus to plus indicates a step edge transition from a target pixel to a background pixel. Ramp edge transitions have sim- ilar sign transitions in the Laplacian along with a zero (or zeros) between the sign changes. Figure l(a) illus- trates the Laplacian computed for step and ramp edge models. The gray level is denoted by f(x), the La- placian by L(x), and the location of the pixels selected are underlined. Note that the average value of the Laplacian is zero at an edge transition. As it is impos- sible to characterize all possible distortions which may occur in edges, only edges which are neither steps nor ramps but are monotonically increasing or de- creasing are considered as distorted edges. A set of rules are developed to compute a more accurate La- placian should such distortions occur in the edges. By first detecting an edge and then combining adjacent Laplacians, the average value L(x) of the Laplacian at an edge transition can be made close to zero as shown by the underlined values in Fig. l(b). The next step involves excluding pixels with low Laplacian magni- tudes by comparing them with an empirically deter- mined threshold. Such pixels typically result from noise and are not edge pixels. In the final step, the connectivity of the remaining pixels is determined by examining the eight-neighborhood of each pixel in order to exclude isolated pixels as the goal is to retain only boundary edge pixels which are generally con- nected to each other. The detailed steps of the selec- tive sampling algorithm are described below: Step 1: Computation of the Laplacian of the image. Compute the Laplacian L[f(x,y)] of the image f(x, y). Let the symbols 0, +, and - represent the zero, positive and negative Laplacian values of L[f(x, y)] and let nO, n +, and n - represent a string of n, (n > 2) continuous zeros, positive, and negative Laplacians. Step 2: Detection of all edges in the image. Scan the rows and columns of L If(x, y)] and use the following rules to detect edges: (a) Noise-free step edge detection A noise-free step edge is characterized by the se- quence: [. . . 0,0,( +, - ), 0,0 ... ] or [ . . .0 ,0 , ( - , +) ,0 ,0 ...3. (b) Noise-free ramp edge detection A noise-free ramp edge is characterized by the sequence: [. . . 0, 0,( +, n0, - ) ,0 ,0 ... ] or [ . . . 0 ,0 , ( - ,n0 , +) ,0 ,0 ...3 318 L. GUPTA and T. SORTRAKUL f(x) L(x) 06 6 -6 _( 0036 6600 030-3 0 -660 (a) Noise-flee step and ramp edges. 6 6 3 0 0 -503 f(x) L(x) Ll(X) 016 1 4-5 50-5 J -L_ 6 0 0 1 6 6 6 1 0 6 6 1 0 0 014 -5 0 -541 0 -5 41 0 0 5 0 -5 0--5 0 _5 0 -5 0 5 f(x) L(x) Ll(x) 0056 0056 6500 665 00 05-4 -1 05-41 -1 -4 50 0 -1 -450 050-5 050-5 -5050 0 -505__0 (b) Distorted edges. Fig. 1. The Laplacian computed for noise-free and distorted edges. (c) Detection of monotonically increasing/decreasing edges A monotonically increasing edge is characterized by the sequence: [... O,O,(n +,kO, m - ),0,0 .... ] and a monotonically decreasing edge is character- ized by [ .... 0,0,(m - ,k0, n +), 0,0 . . . . ]. If a monotonically increasing edge is detected, the Laplacian is updated by summing the n positive and m negative Laplacians and storing them at the first occurrence and last occurrence of positive and nega- tive values, respectively. All other Laplacians are set to zero. Similarly, for the monotonically decreasing edge, the Laplacian is updated by summing the m negative and n positive Laplacians and storing them at the first occurrence and last occurrence of negative and positive values, respectively. It is also possible to have sequences of zeros between a string of positive or negative Laplacians. For this case also, the strings of positive and negative Laplacians are sum- med and stored at the extreme locations. Step 3: Exclusion of small edges. Generate an image .fl (x, y) according to f l (x ,y )= ~1 if ,L[f(x,y)], >Lo, otherwise. Lo is an empirically determined threshold. Step 4: Removal of isolated disconnected pixels. Gen- erate the transformed image f2 (x, y) containing only the selected pixels according to 0 if [ f l (x, y) = 0] or [f l (x, y) = 1 fz(x,y) = and ~f l (x,y) < 2], 8N. f(x, y) otherwise. 8N denotes the set of 8-connected pixels of the pixel at (x, y). Gaussian-mixture-based image 319 Step 5: Histogram off(x, y). Compute the histogram hE (z) off2 (x, y). It is important to note that f2(x, y) is an image which contains the gray levels of the pixels selected from the original imagef(x, y). Given that, typically, only pixels on both sides of an edge transition are selected, the number of pixels used in computing the transformed histogram will be relatively small. The means of the target and background pixels can be expected to be well separated and the proportion of target and background pixels will be approximately equal irrespective of the target and background sizes. t p(z)LITj ~ , ~ p(z)RITj ~.~ P2p(z/C2) PIp(z/CI) hITj @v i Tj ~zlTj ~21Tj z Fig. 2. Illustration of the mixture density, components, and tails. 4. THE GAUSSIAN-MIXTURE PARAMETER ESTIMATION (GMPE) ALGORITHM The parameter-estimation algorithm focuses on the tails of the mixture density to first estimate the two Gaussian components of the mixture density and then estimate the parameters of the model from the two components. For a given threshold T j, the notations used in the formulation are as follows: p(z) = mixture density p(z/Ci) = mixture component i, i = 1, 2 Pip(z/Ci) = weighted mixture component i, i = 1, 2 p(z)L[Tj = left truncated component p(z)R[Tj = right truncated component Li] Tj = left tail of weighted component Pip(z/Ci), i, i = 1, 2 R~[Tj = right tail of weighted component Pip(z/Ci), i, i = 1, 2. The mixture density, the components, and the tails are illustrated in Fig. 2. The parameter estimation is done iteratively, and each iteration consists of three estimation steps. For a given threshold Tj, the para- meters ]21 and ]22 are estimated initially using the truncated components p(Z)L[ Tj and p(Z)R[ T~. The next step involves estimating the weighted compo- nents P1 P(z/C1) and P2 P (z/C2) and the proportions P~ and P2. Finally, the threshold is estimated from the parameters of the estimated mixture components p(z/C1)l Tj and p(z/C2)[Tj. In the formulation of the parameter-estimation algorithm, it is assumed that the contribution of tail L21Tj is negligible in tail L1 [ Tj and that tail R1 [ Tj is negligible in tail R2I Tj. Then, using the symmetrical property of the Gaussian function, the tails R~[Tj and L2[ Tj can be estimated from the tails of the mixture p(z). The detailed steps of the estimation algorithm are described next. Step 0: Selection of initial threshold To. Let the initial estimate of the optimum threshold To be ]20, where ]20 is the mean of p(z). Step 1: Estimation of