Next Article in Journal
Low-Cost Nested-MIMO Array for Large-Scale Wireless Sensor Applications
Next Article in Special Issue
Modeling the Footprint and Equivalent Radiance Transfer Path Length for Tower-Based Hemispherical Observations of Chlorophyll Fluorescence
Previous Article in Journal
Development of Data Registration and Fusion Methods for Measurement of Ultra-Precision Freeform Surfaces
Previous Article in Special Issue
Remote Sensing for Crop Water Management: From ET Modelling to Services for the End Users
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multilook SAR Image Segmentation with an Unknown Number of Clusters Using a Gamma Mixture Model and Hierarchical Clustering

Institute for Remote Sensing Science and Application, School of Geomatics, Liaoning Technical University, Fuxin 123000, China
*
Author to whom correspondence should be addressed.
Sensors 2017, 17(5), 1114; https://doi.org/10.3390/s17051114
Submission received: 29 March 2017 / Revised: 30 April 2017 / Accepted: 10 May 2017 / Published: 12 May 2017

Abstract

:
This paper presents a novel multilook SAR image segmentation algorithm with an unknown number of clusters. Firstly, the marginal probability distribution for a given SAR image is defined by a Gamma mixture model (GaMM), in which the number of components corresponds to the number of homogeneous regions needed to segment and the spatial relationship among neighboring pixels is characterized by a Markov Random Field (MRF) defined by the weighting coefficients of components in GaMM. During the algorithm iteration procedure, the number of clusters is gradually reduced by merging two components until they are equal to one. For each fixed number of clusters, the parameters of GaMM are estimated and the optimal segmentation result corresponding to the number is obtained by maximizing the marginal probability. Finally, the number of clusters with minimum global energy defined as the negative logarithm of marginal probability is indicated as the expected number of clusters with the homogeneous regions needed to be segmented, and the corresponding segmentation result is considered as the final optimal one. The experimental results from the proposed and comparing algorithms for simulated and real multilook SAR images show that the proposed algorithm can find the real number of clusters and obtain more accurate segmentation results simultaneously.

1. Introduction

Synthetic Aperture Radar (SAR) is widely used in civil and military fields because of its characteristics such as high resolution, strong penetration, full time and full weather [1]. However, its specific imaging mechanism leads to speckle noises in SAR images. Although the multilook technique can partly reduce the type of noises, there are still obvious noises remaining in SAR images [2]. Image segmentation, as a crucial step for image analysis and interpretation, involves separating pixels into clusters, so that the pixels in the same cluster are homogeneous and the pixels in different clusters are heterogenous [3]. The main tasks for image segmentation are to determine a suitable number of clusters and obtain an accurate segmentation result [4]. Due to the inherent speckle noise existing in SAR images, these tasks are difficult to fulfill. Therefore, the objective of the proposed algorithm is to obtain the real number of clusters and the accurate segmentation result simultaneously for speckled SAR images. In this paper, the real number is regarded as the number selected by human visual judgment.
Many approaches have been studied to determine the number of clusters, such as clustering methods [5,6,7,8,9] and statistic methods [10,11,12,13,14]. The most popular clustering method is Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA) [5], in which the number of clusters is changed by alternately splitting and merging clusters. However, it requires a lot of parameters given by humans, and the K-means dominating optimization in ISODATA algorithm is not suitable for segmenting SAR images because of the statistical characteristics of the speckle noises. In order to reduce the human influence, another clustering method, the hierarchical clustering included agglomerative procedure (bottom-up) and divisive procedure (top-down) [7], was considered. The agglomerative produce is the most commonly used one and the representative algorithm is Agglomerative Nesting (AGNES) [8], which considers each element as a cluster and the two clusters which are closest to each other with minimum distance are merged at each step. This procedure continues until only one cluster exists. In order to overcome the sensitivity of noises to Euclidean distance, Legendre et al. [9] used the average distance to measure the closeness between two clusters, but the smaller cluster is easily ignored. Generally, statistical methods are considered to be a more appropriate strategy in SAR image segmentation [10,11]. Green [12] proposed the Reversible Jump Markov Chain Monte Carlo (RJMCMC) algorithm which has been successfully applied to image segmentation with an unknown number of clusters. Based on the work of Green, Kato [13] combined a mixture model and RJMCMC to complete the segmentation. It regards the number of clusters as a random variable contained in model parameters. The number jumps by sampling operation such as birth or death class, and the optimal number is obtained by maximizing or minimizing the objective function. However, the statistical methods for image segmentation with unknown number of clusters are too cumbersome because the dimension of parameters always changes along with the changing of the number of clusters, and it easily falls into the local optimum. To sum up, hierarchical clustering is convenient for obtaining the number of clusters, and the statistical model can describe images well [14], so hierarchical clustering combining a statistical model is regarded as a more appropriate approach in SAR image segmentation in this paper.
There are many approaches of hierarchical theory with statistical models. Mignotte et al. [15] proposed an unsupervised segmentation method for sonar images based on a hierarchical Markov Random Field (MRF) model; it defines a causal energy function with a Gibbs distribution form. Although it can avoid the ignoring of clusters with a few pixels, it cannot be applied to SAR image segmentation in the absence of the capability to manage the statistical property of the speckle noises. Upgrading the MRF model, Wang et al. [16] proposed a new unsupervised method based on a hierarchical Triplet Markov Field (TMF) model with Gaussian distribution to segment SAR images. It combines the advantages of both hierarchical and TMF in the discrete wavelet domain to deal with the problem of insufficient local statistical information. However, the unimodal characteristic model is not as flexible as a mixture model [17]. Garcia et al. [18] proposed a hierarchical Gaussian Mixture Model (GMM) algorithm, which is able to automatically learn the optimal number of components for the simplified GMM and successfully segment natural images, but the GMM cannot describe the data with non-Gaussian statistical character [19], such as SAR images. Oliver et al. [20] has proved that the intensities of pixels satisfy a Gamma distribution according to the speckle noise existing in SAR images, Dong et al. [21] fully validated that the Gamma distribution with MRF is more suitable for SAR image segmentation than a Gaussian distribution with MRF, so the hierarchical clustering based on a Gamma mixture model (GaMM) is proposed for SAR image segmentation.
In this paper, an over-segmentation result with the number of clusters which is much larger than the real number is given as the initial status. The number of clusters is reduced gradually by a cluster merging operation until the optimal number of clusters is obtained. Firstly, the original image is divided into many subsets in gray level by a span parameter. Each subset is regarded as a cluster, and the two subsets with minimum global energy are merged until the number of clusters is equal to one. In addition, the segmentation model with current number of clusters is established by marginal probability distribution based on the GaMM. The negative logarithm of marginal probability is described as the global energy function, and the smaller the value is, the better the segmentation is. Furthermore, in order to consider the influence of neighboring pixels, the prior probability is defined by the conditional probability of the membership based on MRF. Then, the segmentation results with the current number of clusters can be obtained by posterior probability according to the Bayesian theory [22], the number of clusters reduces by one after each iteration, and a segmentation result and a global energy are obtained under the current number of clusters. Finally, the optimal number of clusters is the one corresponding to the minimum global energy, and the result corresponding to the optimal number is the final optimal segmentation result.

2. Description of the Proposed Algorithm

Given a multilook SAR image x = {xi, i = 1, ..., n}, where i is the index of pixels, n is the number of pixels, xi is the intensity of pixel i. Assume that the image has m clusters, in order to represent the image segmentation, membership matrix z = [zij]n×m is used, where j is the index of clusters, zij ∈ {0, 1}, j = 1 m z i j = 1 , zij = 1 indicates the ith pixel belongs to jth cluster. Furthermore, z can be re-expressed as z = (zi·, i = 1, …, n)Τ = (z·j, j = 1, ..., m), where zi· = (zij, j = 1, ..., m) and z·j = (zij, i = 1, ..., n)Τ, where zi· and z·j are row and column vectors of membership matrix z, respectively. In the image segmentation, x is observable data, z can be considered as unobservable data, {x, z} is the complete data set. The image segmentation is the process of estimating unobservable data z given x.

2.1. Segmentation Process

Assume that the conditional distribution of ith pixel belonging to jth cluster is Gamma distribution and its conditional probability density function (pdf) is:
p ( x i | z i j = 1 ;   θ j ) = Ga ( x i ;   α j , β j ) = x i α j 1 Γ ( α j ) β j α j exp ( x i β j )
where Ga(·) is the pdf of Gamma distribution, θj = (αj, βj) is distribution parameter vector of jth cluster, αj and βj are shape and scale parameters of Gamma distribution, respectively.
According to the Bayesian theory, the joint pdf of xi and zij is:
p ( x i , z i j ;   θ j ) = p ( z i j = 1 ) p ( x i | z i j = 1 ;   θ j )
The marginal probability of xi is:
p ( x i ;   θ ) = j = 1 m p ( z i j = 1 ) p ( x i | z i j = 1 ;   θ j )
where θ = {θj, j = 1, …, m}.
Assume that the marginal distribution of intensity in the SAR image z is independent of each other, then the marginal probability of image x is:
p ( x ;   θ ) = i = 1 n p ( x i ;   θ ) = i = 1 n j = 1 m p ( z i j = 1 ) p ( x i | z i j = 1 ;   θ j )
The prior probability for zij = 1 is denoted p(zij = 1) = πij, where πij ∈ [0, 1] and satisfies j = 1 m π i j = 1 . In order to consider the spatial interaction from neighborhood pixels, the prior probability is defined by a conditional distribution:
π i j = p ( z i j = 1 | z i j N i ) = exp ( η i N i z i j ) j = 1 m exp ( η i N i z i j )
where Ni is the neighborhood pixels set of ith pixel, η is the coefficient indicates the influence of intensities of neighbor pixels, the higher the coefficient η is, the stronger the neighbor influence is.
Combining Equations (4) and (5), the marginal probability can be reformed as:
p ( x ;   θ ) = i = 1 n j = 1 m π i j Ga ( x i ;   α j , β j )
As described in Equation (6), the marginal distribution of observable data x is GaMM, and the prior probability πij is the model weight. The number of clusters is equal to the number of GaMM components.
Equation (6) is regarded as the segmentation model expressed by GaMM. The higher the marginal probability is, the better the segmentation result is. There are many criteria to determine the optimal segmentation. Among them, the maximum posterior probability is one of the most commonly used ones. According to the Bayesian theory, the posterior probability of ith pixel belonging to jth cluster is:
p ( z i j = 1 | x i ;   θ ) = p ( z i j = 1 ) p ( x i | z i j = 1 ;   θ j ) j = 1 m p ( z i j = 1 ) p ( x i | z i j = 1 ;   θ j ) = π i j Ga ( x i ;   α j , β j ) j = 1 m π i j Ga ( x i ;   α j , β j )
Finally, the optimal segmentation corresponds to the resolution of the maximum posterior probability and can be expressed as:
z i j = { 1   if   p ( z i j = 1 | x i ;   θ j ) = max { p ( z i j = 1 | x i ; θ j ) ,   j = 1 ,   ,   m } 0   otherwise
From Equation (8) one can see that the segmentation is the problem of solving the model parameters, θj = {αj, βj}. For multilook SAR images, the shape parameter αj is usually equal to the number of looks. For scale parameter βj, it is difficult to solve by Equation (6), thus, the global energy function L defined as the negative logarithm of marginal probability of x is used instead of Equation (6) to evaluate the quality of results. The smaller the function value is, the bigger the marginal probability is, and the better the segmentation result is:
L = log   p ( x ;   θ ) = log i = 1 n j = 1 m π i j Ga ( x i ;   β j ) = i = 1 n log { j = 1 m π i j Ga ( x i ;   β j ) }
Because the scale parameter βj is differentiable function of global energy L, the derivative of L for βj is:
L β j = i = 1 n 1 j = 1 m π i j Ga ( x i ;   β j ) π i j Ga ( x i ;   β j ) β j = i = 1 n π i j Ga ( x i ;   β j ) j = 1 m π i j Ga ( x i ;   β j ) ( α β j + x i β j 2 ) = 1 β j i = 1 n p ( z i j = 1 | x i , θ j ) ( α + x i β j )
Let Equation (10) be equal to 0, its solution is equivalent to the solution of Equation (11):
i = 1 n p ( z i j = 1 | x i ,   β j ) ( α + x i β j ) = 0
It is difficult to solve Equation (11), because βj exists in both of the first and second term. In order to obtain the best approximate solution, let βj in the first term equal to the value of the last iteration expressed as t, Equation (11) can be rewritten as:
i = 1 n p ( z i j ( t ) = 1 | x i , β j ( t ) ) ( α + x i β j ( t + 1 ) ) = 0
Then the scale parameter β j ( t + 1 ) can be solved as:
β j ( t + 1 ) = i = 1 n p ( z i j ( t ) = 1 | x i ; β j ( t ) ) x i i = 1 n p ( z i j ( t ) = 1 | x i ; β j ( t ) ) α
The optimization process of the segmentation results under fixed m will be stopped when t reaches the maximum value of T. According to our experiments, it is found that taking T = 20 is suitable for all testing images, so in the proposed algorithm, set T = 20 as the iterative stopping condition.
In summary, if the initial segmentation zij(0) under the current number of clusters m is known, the segmentation process regarded as an inner loop under fixed number of clusters is described as follows:
  • Calculating the prior probability πij(t) by Equation (5);
  • Calculating the scale parameter βj(t) by Equation (13), and calculating conditional probability by Equation (1);
  • Calculating the posterior probability p(zij = 1|xi; θj)(t) by Equation (7);
  • Repeating Step 1–3 until t is equal to T;
  • The optimal segmentation result zij(T) under fixed number of clusters m is obtained by Equation (8).

2.2. Hierarchical Clustering Process

In the above Subsection, the optimal segmentation can be obtained under a fixed number of clusters. In order to obtain the optimal result with a real number of clusters, the number of clusters m is viewed as a variable in the proposed hierarchical clustering algorithm. The number m gradually decreases from an initial number which is much greater than the real number of clusters to one by a cluster merging operation. Thus, one can get the optimal results under each number of clusters by inner loop. At last, the optimal segmentation is determined from the series of optimal results by minimizing global energy defined by Equation (9).
If τ is an iterative indicator in the outer loop for hierarchical clustering process, the key point of the hierarchical clustering process is to obtain the initial parameters at the τ + 1 iteration, such as, z(τ+1, 0), p(zij = 1|xi; θj)(τ+1, 0). Further, the optimal results at τ + 1 iteration z(τ+1, T) can be obtained by inner loop. If z·j = {xi; zij = 1}(all pixels in jth cluster) expresses the collection of jth class, and the segmentation result in τ iterative is z(τ, T) = {z·1(τ, T), ..., z·j−1(τ, T), z·j(τ, T), z·j+1(τ, T), ..., z·j’−1(τ, T), z·j(τ, T), z·j’+1(τ, T), ..., z·m(τ, T)}. Assuming that the merged clusters are j and j’, then, z(jj’)(τ, *) = {z·1(τ, *), ..., z·j−1(τ, *), z·j(τ, *), z·j+1(τ, *), ..., z·j’−1(τ, *), z·j(τ, *), z·j’+1(τ, *), ..., z·m−1(τ, *)} ,where z·j(τ, *) = z·j(τ, T) + z·j(τ, T), others stay the same and {z·j(τ, *), z·j’+1(τ, *), ..., z·m−1(τ, *)} are the rearrange of {z·j’+1(τ, T), ..., z·m(τ, T)}. Further, the new global energy function by merging j and j’ is expressed as Ljj(τ, *), let ΔL(τ) = Ljj(τ, *), j = 1, ..., m − 1, j’ = j + 1, ..., m, choose the clusters j ^ and j ^ which corresponding to the minimum global energy function, that is:
{ j ^ , j ^ } = arg min { Δ L ( τ ) : L j j ( τ ,   * ) , j = 1 , , m 1 , j = j + 1 , , m }
The segmentation result in τ + 1 iteration is z ( τ + 1 ,   0 ) = z ( j ^ j ^ ) ( τ ,   * ) , let m(τ+1) = m(τ) − 1, and the initial posterior probability in τ + 1 iteration p(zij = 1|xi; θj)(τ+1, 0) is:
p ( z i j = 1 | x i ;   θ j ) ( τ + 1 ,   0 ) = p ( z i j ^ = 1 | x i ;   θ j ^ ) ( τ ,   T ) + p ( z i j ^ = 1 | x i ;   θ j ^ ) ( τ ,   T )
where p ( z i j ^ = 1 | x i ;   θ j ^ ) ( τ + 1 ,   0 ) does not exist now.
After merging, z(τ+1, T) can be obtained by inner loop and the global energy function value L(τ+1) is recorded by new parameters πij(τ+1, T), βj(τ+1, T) at each iteration, then the real number of clusters m ^ can be estimated by:
{ m ^ , τ ^ } = arg max { L ( τ ) , τ = 0 ,   ,   m ( 0 ,   0 ) }
The final optimal segmentation is:
z ^ i j = z i j ( τ ^ ,   T )
In summary, an initial segmentation must be given before the process of both segmentation and hierarchical clustering. The initial segmentation z(τ=0, t=0) is given by dividing the image x into many sub-classes using span parameter d. Then the ith pixel in jth cluster can be expressed as:
z i j ( 0 ,   0 ) = x i d
where ⌈·⌉ expresses for rounding up to an integer. The jth cluster will be deleted if there are no pixels in it until all clusters are assigned a number of pixels. Then, m(0, 0) is equal to the column number of z(0, 0), and the initial scale parameter βj(τ=0, t=0) is calculated by:
β j ( 0 ,   0 ) = 1 α n j x i z j ( 0 , 0 ) x i
where nj is the number of pixels in jth cluster.

2.3. General Procedure of the Proposed Algorithm

The proposed algorithm uses the outer loop with hierarchical clustering to merge clusters in order to find the optimal number of clusters and the inner loop with the GaMM as the segmentation model to obtain the optimal segmentation under current number of clusters. The outer loop stops until the class number is reduced to one. The optimal segmentation under the estimated real number of clusters can be obtained by minimizing the global energy defined in Equation (9).
The process of the proposed algorithm can be summarized as follows:
  • Initialization. The iterative indicator of inner loop t = 0 and outer loop τ = 0, the maximum number of inner loop T = 20 according to the experimental experience, shape parameter α which is equal to the looks of the SAR images, intensity of neighbor influence η within the range of 0~1;
  • Coarse segmentation. The initial segmentation zij(0, 0) is obtained by Equation (18), and the initial scale parameter βj(0, 0) is calculated by Equation (19).
  • Inner loop. Calculating the prior probability πij(τ, t) by Equation (5), calculating the posterior probability p(zij = 1|xi; θj)(τ, t) by Equation (7). In the inner iteration, the scale parameter βj(τ, t) is calculated by Equation (13), after T inner iterations, the segmentation zij(τ, T) under the current outer loop is obtained and L(τ) is recorded by Equation (9);
  • Merging clusters. Choosing clusters j and j’ by Equation (14) which correspond to the minimum global energy function. Merging clusters j and j’, and updating the parameters of πij(τ, *), βj(τ, *) by Equations (5) and (13), respectively;
  • Updating parameters. After merging, let m(τ+1) = m(τ) − 1, and βj(τ+1, 0) is obtained by Equation (15), then the segmentation zij(τ+1, T) and global energy L(τ+1) are obtained by inner loop;
  • Repeat Step 4–5, and stop the outer loop when m(τ) = 1;
  • The real number of clusters m ^ and final optimal segmentation result z ^ i j is found by Equations (16) and (17).

3. Experimental Results and Discussion

The proposed algorithm is tested with a simulated SAR image and four real SAR images. In addition, the comparing algorithms, including ISODATA [5], improved AGNES [8], HMRF-FCM [6] and Gamma-MRF [21] algorithms are also used to evaluate the proposed algorithm.

3.1 Simulated SAR Image Segmentation

Figure 1a is a template with size of 128 × 128 pixels for generating a simulated SAR image which covers four homogeneous regions labeled I-IV, respectively. Figure 1b is the simulated image, in which the intensities of pixels in each homogeneous region are drawn from a Gamma distribution with the shape and scale parameters listed in Table 1.
In order to generate initial segmentation of image shown in Figure 1b, we set span parameter d = 30. Then the initial number of clusters m is ⌈255/30⌉ which is equal to 9. In this experiment, α = 4 as the looks of SAR image. The coefficient η, which expresses the influence of neighbor pixels, is always taken from 0 to 1. If η is too small, the results will exist many misclassification pixels although the number of clusters is correct, contrarily, it will be over segmented. According to experiments, η = 0.8 which is suitable for the simulated images. Figure 2a is the initial segmentation with m = 9. Figure 2b–h show the corresponding segmentation results with the gradually decreased number of clusters during outer loop. It shows that the initial segmentation in gray level already can roughly distinguish each homogeneous regions shown in Figure 2a. With the increase of the iteration, the number of clusters reduce one each time, and the optimal segmentation can be obtained under the condition of each number of clusters shown in Figure 2b–h. For the homogeneous region with smaller scale parameters like region I, it is easy to be segmented, and the heavily interwoven regions like regions II and III must be divided by iterations with the help of MRF. In the iterations, the global energy value L and the number of clusters m with iteration τ are shown in Figure 3. Figure 3a,b are the full view and local amplification image. It can be seen that the minimum global energy function value is 8.15465 × 104 corresponding to the iteration 6, and the corresponding number of clusters is 4. Thus, the final optimal segmentation result is the one shown in Figure 2f.
The final optimal segmentation results by the proposed algorithm are displayed in Figure 4. Figure 4a is the final optimal segmentation, Figure 4b is the outline of the homogeneous regions, and Figure 4c shows the result of outlines overlaid on the original image.
One can see from Figure 4 that the proposed algorithm can segment each homogeneous region well. There are almost no misclassified pixels, and the boundary is smooth and accurate. For qualitative analysis of the proposed algorithm, the segmentation results of the comparing algorithms are displayed in Figure 5, where Figure 5a–d are the segmentation results of ISODATA, AGNES, HMRF-FCM and Gamma-MRF, respectively. It shows that the ISODATA and AGNES algorithms can find the right number of clusters but ISODATA need more prior condition to limit the range of results and the result is pretty bad. Although AGNES avoids the help of prior knowledge, its segmentation effect is equivalent to ISODATA, there are many misclassified pixels in Figure 5a,b. The HMRF-FCM algorithm can reduce the influence of noise by introducing MRF to model the relation between neighbor pixels. It still cannot distinguish the homogeneous region well. For the Gamma-MRF algorithm, the Gamma distribution combining MRF theory can accurately describe the distribution of pixel intensities, but there are also many pixel errors. From the visual viewpoint, the proposed algorithm is much better than the compared ones.
In order to demonstrate the superiority of GaMM, the histogram for simulated image and the estimated curve with Gaussian distribution, Gamma distribution and GaMM, respectively, are shown in Figure 6. From Figure 6, one can see that the GaMM with different weights can fit the histogram of the image more accurately.
To evaluate the proposed and comparing algorithms in quantitative, the main parameters, running time, accuracy and Kappa coefficient are shown in Table 2. For the proposed algorithm, the running time is 1139.19 s. The user’s accuracy and producer’s accuracy of each homogeneous region are more than 98%, and the overall accuracy and kappa coefficient are 99.34% and 0.99, respectively. For the comparing algorithms, the producer’s accuracy of region III is only 36.16% and 40.28 % in ISODATA and AGNES, and the overall accuracy and kappa coefficient are about 60% and 0.5, respectively. HMRF-FCM and Gamma-MRF are better than the two algorithms discussed above, but the accuracy is still lower than that for the proposed algorithm. Generally, the running time of statistical model is usually much higher than clustering methods, but compared with Gamma-MRF, the proposed algorithm greatly reduces the time. Thus, Table 2 precisely describes the effectiveness and robustness of the proposed algorithm.
Furthermore, the mean and variance listed in Table 3 are also used to evaluate the quality of the proposed algorithm. One can see that the deviation rates between the estimated and actual values for all homogeneous regions are very small. Table 3 illustrates the accuracy of the algorithm.
In order to analyze the robustness of the method and the segmentation accuracy affected by different span parameter d and neighborhood coefficient η, the relationship between the number of clusters, overall accuracy, running time and span parameter d, neighborhood coefficient η are given in Table 4. When the neighborhood coefficient is small, the number of clusters obtained are less than the real one. On the contrary, it is larger than the real one. In this part, we only consider the overall accuracy under the correct number of clusters. We can see that the accuracy of d = 30 is almost higher than others under the same neighborhood coefficient. In Table 4, when d = 30, η = 0.8, the overall accuracy is the highest.
In order to prove that T = 20 is the optimal result, the change of scale parameter β with iteration is shown in Figure 7. From Figure 7 one can see that scale parameter β for each cluster converges to its stable value after 20 iterations. Consequently, 20 is regard as T’s optimal choice, because it not only ensures the accuracy of segmentation but also saves time. Furthermore, the sensitivity of the value of T is also analyzed through the segmentation accuracy with different T under d = 30, η = 0.8, the results are listed in Table 5. If the value of T is too small, the segmentation result is inaccurate, if the value of T is too large, the running time will be increased.

3.2. Real SAR Images Segmentation

Figure 8 shows the real SAR images with size of 128 × 128 pixel and the real number of clusters are 2, 3, 3 and 4, respectively. Figure 8a is a RADARSAT-I intensity image with HH polarization and the spatial resolution of 30 m, Figure 8b,c are RADARSAT-II images with VV polarization and the spatial resolution of 25 m, Figure 8d is a RADARSAT-I image with HH polarization and the spatial resolution of 25 m. All of them are sea ice images in different periods, from light to dark, representing sea ice with density from high to low. The image in Figure 8a with two clusters is easy to segment for all algorithms. For the images in Figure 8b,c with three clusters, it is difficult to segment the different types of sea ices. For the image in Figure 8d with four clusters, the difficulty is that the water is easy to ignore because its area is too small compared with the background of low density sea ice areas.
Figure 9 is the segmentation results with the proposed algorithm and the comparing algorithms of ISODATA, AGNES, HMRF-FCM and Gamma-MRF. Figure 9(a1–d1) are the results of the ISODATA algorithm, Figure 9(a2–d2) are the results of the AGNES algorithm, Figure 9(a3–d3) are the results of the HMRF-FCM algorithm, Figure 9(a4–d4) are the results of the Gamma-MRF algorithm and Figure 9(a5–d5) are the results of the proposed algorithm.
It can be seen that the ISODATA algorithm cannot always determine the right number of clusters, and as shown as Figure 9(c1), the real number of the original image is 3 while the estimated number is 2, which is a fatal error for image segmentation.
Although the AGNES algorithm can find the right number of clusters, it cannot overcome the noise effect well and ignores the existence of water class, see Figure 9(d2).
Although the HMRF-FCM algorithm combines a fuzzy cluster based on Gaussian distribution and MRF, it also cannot obtain a satisfactory segmentation result.
The Gamma-MRF algorithm adopts the Gamma distribution to describe the data, and its results are not good enough either. Besides, both HMRF-FCM and Gamma-MRF lack a mechanism for changing the number of clusters, which makes the algorithm not flexible enough in applications.
As shown in Figure 9(a5–d5), the segmentation results from the proposed algorithm are all the optimal segmentation results under a real number of clusters. One can see that the sea ice images during different periods are evidently distinguished all the time by the proposed algorithm, and there are almost no mis-segmentations in the homogeneous regions.
In addition, the relation of the given parameter, global energy value and the number of clusters are shown in Figure 10. The minimum global energy value is 8.00387 × 104, 8.93877 × 104, 8.88819 × 104 and 7.20323 × 104 corresponding to the iterations 8, 7, 11, 6, and the corresponding number of clusters are 2, 3, 3, 4, respectively. It illustrates that Figure 9(a5–d5) with the proposed algorithm under the real number of clusters are the perfect results.
The values of estimated parameters and running time for real images segmentation results are shown in Table 6. Combining with the results in Figure 9 and the running time in Table 6, it can be seen that the running time of the proposed algorithm has a high cost performance.
In addition, segmentation experiments for SAR images covering urban and agricultural areas were also carried out. Figure 11(a1,b1) show the HH RADARSAT-I/II intensity images for an urban area with two looks and an agricultural area with three looks, respectively. The scene shown in Figure 11(a1) contains three objects, including water, streets and buildings, and the scene in Figure 11(b1) contains water and two kinds of farmland. Figure 11(a2,b2) are the corresponding segmentation results from the proposed algorithm with d = 20, η = 0.3, α = 2 and d = 20, η = 0.35, α = 3, respectively. One can see that the proposed algorithm can separate the street from water on the basis of segmenting buildings and effectively distinguish different kinds of farmland. The water and streets shown in Figure 11(a2) are not clearly displayed, just because their grayscales are similar. Figure 11 illustrates that the proposed algorithm can segment not only simple SAR images well, but also complex images with different noise levels.

4. Conclusions

In this paper, a novel hierarchical clustering segmentation algorithm for multilook SAR image segmentation is proposed. The GaMM can describe the statistical characteristics of the image more accurately than unimodal character model. The hierarchical clustering with global energy merging criterion can determine the real number of clusters and obtain the accurate segmentation results simultaneously. The qualitative and quantitative analyses of the simulated and real images with the proposed and comparing algorithms fully demonstrate the superiority of the proposed algorithm in multilook SAR image segmentation. Some improvements of this research are still needed. On the one hand, the coarse segmentation method could be replaced with some more efficient approaches, such as, Fuzzy C-means (FCM). On the other hand, the image segmentation model can be improved by considering texture features.

Acknowledgments

This work was supported by the Liaoning Provincial Nature Science Foundation of China (No. 2015020090), Liaoning Technical University Graduate Education Cultivation Project of China (No. YB201605) and National Natural Science Foundation of China (No. 41301479; No. 41271435).

Author Contributions

Quanhua Zhao and Yu Li conceived and designed the experiment; Xiaoli Li performed the experiment; Quanhua Zhao, Yu Li and Xiaoli Li analyzed the data and wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bombrun, L.; Vasile, G.; Gay, M.; Totir, F. Hierarchical segmentation of polarimetric SAR images using heterogeneous clutter models. IEEE Trans. Geosci. Remote Sens. 2011, 49, 726–737. [Google Scholar] [CrossRef]
  2. Lee, J.S.; Schuler, R.H.; Lang, K.J. K-distribution for multi-look processed polarimetric SAR imagery. Proceedings of 1994 IEEE International Geoscience and Remote Sensing Symposium (IGARSS 1994), Pasadena, CA, USA, 8–12 August 1994; pp. 2179–2181. [Google Scholar]
  3. Shang, R.H.; Tian, P.P.; Jiao, L.C.; Stolkin, R.; Feng, J.; Hou, B.; Zhang, X.G. A spatial fuzzy clustering algorithm with kernel metric based on immune clone for SAR image segmentation. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 9, 1640–1652. [Google Scholar] [CrossRef]
  4. Zhao, Q.H.; Li, Y.; Wang, Y. SAR image segmentation with unknown number of classes combined Voronoi tessellation and RJMCMC algorithm. In Proceedings of the ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Prague, Czech Republic, 12–19 July 2016; pp. 119–124. [Google Scholar]
  5. Fergani, B.; Kholladi, M.K.; Bahri, M. Comparison of FCM and FISODATA. Int. J. Comput. Appl. 2012, 56, 35–39. [Google Scholar] [CrossRef]
  6. Chatzis, S.P.; Varvarigou, T.A. A fuzzy clustering approach toward hidden Markov random field models for enhanced spatially constrained image segmentation. IEEE Trans. Fuzzy Syst. 2008, 16, 1351–1361. [Google Scholar] [CrossRef]
  7. Sharma, A.; Boroevich, K.A.; Shigemizu, D.; Kamatani, Y.; Kubo, M.; Tsunoda, T. Hierarchical maximum likelihood clustering approach. IEEE Trans. Biomed. Eng. 2017, 64, 112–122. [Google Scholar] [CrossRef] [PubMed]
  8. Leonard, K.; Peter, J.R. Agglomerative Nesting (Program AGNES). In Finding Groups in Data: An Introduction to Cluster Analysis; John Wiley & Sons: Brussels, Belgium, 2005; Volume 5, pp. 199–252. [Google Scholar]
  9. Legendre, P.; Legendre, L. Cluster analysis. In Numerical Ecology, 2nd ed.; Elsevier: Amsterdam, The Netherlands, 1998; Volume 8, pp. 303–381. [Google Scholar]
  10. Lopes, A.; Laur, H.; Nezry, E. Statistical distribution and texture multilook and complex SAR images. In Proceedings of the 10th Annual International Symposium on Geoscience and Remote Sensing (IGARSS 1990), Baltimore, MD, USA, 20–24 May 1990; pp. 2427–2430. [Google Scholar]
  11. Tison, C.; Nicolas, J.M.; Tupin, F.; Maître, H. A new statistical model for Markovian classification of urban areas in high-resolution SAR images. IEEE Trans. Geosci. Remote Sens. 2004, 42, 2046–2057. [Google Scholar] [CrossRef]
  12. Green, P.J. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 1995, 82, 711–732. [Google Scholar] [CrossRef]
  13. Kato, Z. Segmentation of color images via reversible jump MCMC sampling. Image Vision Comput. 2008, 26, 361–371. [Google Scholar] [CrossRef]
  14. Lee, S.; Crawford, M.M. Unsupervised multistage image classification using hierarchical clustering with a Bayesian similarity measure. IEEE Trans. Image Process. 2005, 14, 312–320. [Google Scholar] [PubMed]
  15. Mignotte, M.; Collet, C.; Perez, P.; Bouthemy, P. Sonar image segmentation using an unsupervised hierarchical MRF model. IEEE Trans. Image Process. 2000, 9, 1216–1231. [Google Scholar] [CrossRef] [PubMed]
  16. Wang, J.J.; Jiao, S.H.; Shen, L.Y.; Sun, J.Y.; Tang, L. Unsupervised SAR image segmentation based on a hierarchical TMF model in the discrete wavelet domain for sea area detection. Discret. Dyn. Nat. Soc. 2014, 2014, 1–9. [Google Scholar] [CrossRef]
  17. Peng, Y.J.; Chen, J.Y.; Xu, X.; Pu, F.L. SAR Images statistical modeling and classification based on the mixture of Alpha-Stable distributions. Remote Sens. 2013, 5, 2145–2163. [Google Scholar] [CrossRef]
  18. Garcia, V.; Nielsen, F.; Nock, R. Hierarchical Gaussian mixture model. In Proceedings of the 2010 IEEE International Conference on Acoustics Speech, and Signal Processing (ICASSP 2010), Dallas, TX, USA, 14–19 March 2010; pp. 4070–4073. [Google Scholar]
  19. AlSaeed, D.H.; Bouridane, A.; ElZaart, A.; Sammouda, R. Two modified Otsu image segmentation methods based on Lognormal and Gamma distribution models. In Proceedings of the 2012 International Conference on Information Technology and e-Services, Sousse, Tunisia, 24–26 March 2012; pp. 1–5. [Google Scholar]
  20. Oliver, C.J.; Quegan, S. Understanding Synthetic Aperture Radar Images; SciTech Publishing Incorporated: Raleigh, NC, USA, 2004. [Google Scholar]
  21. Dong, Y.; Forester, B.; Milne, A. Evaluation of radar image segmentation by Markov random field model with Gaussian distribution and Gamma distribution. In Proceedings of the 1998 IEEE International Geoscience and Remote Sensing, Seattle, WA, USA, 6–10 July 1998; pp. 1617–1619. [Google Scholar]
  22. Nikou, C.; Likas, A.C.; Galatsanos, N.P. A Bayesian framework for image segmentation with spatially varying mixtures. IEEE Trans. Image Process. 2010, 19, 2278–2289. [Google Scholar] [CrossRef] [PubMed]
Figure 1. (a) Template; (b) Simulated image.
Figure 1. (a) Template; (b) Simulated image.
Sensors 17 01114 g001
Figure 2. Segmentation process images in iterations. (a) Initial, m = 9; (b) τ = 1, m = 8; (c) τ = 2, m = 7; (d) τ = 3, m = 6; (e) τ = 4, m = 5; (f) τ = 5, m = 4; (g) τ = 6, m = 3; (h) τ = 7, m = 2.
Figure 2. Segmentation process images in iterations. (a) Initial, m = 9; (b) τ = 1, m = 8; (c) τ = 2, m = 7; (d) τ = 3, m = 6; (e) τ = 4, m = 5; (f) τ = 5, m = 4; (g) τ = 6, m = 3; (h) τ = 7, m = 2.
Sensors 17 01114 g002
Figure 3. Global energy value L and the number of clusters m at each iteration. (a) For the whole m; (b) For part of m amplification.
Figure 3. Global energy value L and the number of clusters m at each iteration. (a) For the whole m; (b) For part of m amplification.
Sensors 17 01114 g003
Figure 4. Segmentation results of the proposed algorithm. (a) Final optimal result; (b) Outline; (c) Superposition.
Figure 4. Segmentation results of the proposed algorithm. (a) Final optimal result; (b) Outline; (c) Superposition.
Sensors 17 01114 g004
Figure 5. Segmentation results of comparing algorithms with (a) ISODATA; (b) AGNES; (c) HMRF-FCM; (d) Gamma-MRF.
Figure 5. Segmentation results of comparing algorithms with (a) ISODATA; (b) AGNES; (c) HMRF-FCM; (d) Gamma-MRF.
Sensors 17 01114 g005
Figure 6. Histogram for simulated image and estimated curve with Gaussian distribution, Gamma distribution and Gamma mixture model.
Figure 6. Histogram for simulated image and estimated curve with Gaussian distribution, Gamma distribution and Gamma mixture model.
Sensors 17 01114 g006
Figure 7. The changing curve of scale parameters in each clusters in inner loop iterations.
Figure 7. The changing curve of scale parameters in each clusters in inner loop iterations.
Sensors 17 01114 g007
Figure 8. Original real SAR images. (ad) are sea ice images with 2, 3, 3, 4 clusters.
Figure 8. Original real SAR images. (ad) are sea ice images with 2, 3, 3, 4 clusters.
Sensors 17 01114 g008
Figure 9. Segmentation results of real SAR image with the proposed and comparing algorithms. (a1d1) ISODATA; (a2d2) AGNES; (a3d3) HMRF-FCM; (a4d4) Gamma-MRF; (a5d5) The proposed algorithm.
Figure 9. Segmentation results of real SAR image with the proposed and comparing algorithms. (a1d1) ISODATA; (a2d2) AGNES; (a3d3) HMRF-FCM; (a4d4) Gamma-MRF; (a5d5) The proposed algorithm.
Sensors 17 01114 g009aSensors 17 01114 g009b
Figure 10. Global energy value L and the number of clusters m at each iteration. (a) Original SAR image 1 with d = 30, η = 0.5; (b) Original SAR image 2 with d = 30, η = 0.4; (c) Original SAR image 3 with d = 20, η = 0.7; (d) Original SAR image 4 with d = 30, η = 0.5.
Figure 10. Global energy value L and the number of clusters m at each iteration. (a) Original SAR image 1 with d = 30, η = 0.5; (b) Original SAR image 2 with d = 30, η = 0.4; (c) Original SAR image 3 with d = 20, η = 0.7; (d) Original SAR image 4 with d = 30, η = 0.5.
Sensors 17 01114 g010aSensors 17 01114 g010b
Figure 11. Segmentation results of urban and agricultural SAR images with the proposed algorithm, where (a1,b1) are the original SAR images, (a2,b2) are the segmentation results from the proposed algorithm.
Figure 11. Segmentation results of urban and agricultural SAR images with the proposed algorithm, where (a1,b1) are the original SAR images, (a2,b2) are the segmentation results from the proposed algorithm.
Sensors 17 01114 g011
Table 1. Gamma distribution parameters for each homogeneous region.
Table 1. Gamma distribution parameters for each homogeneous region.
ParametersHomogeneous Region
IIIIIIIV
α4444
β5203565
Table 2. Main parameters, running time, accuracy and Kappa coefficient of simulated image.
Table 2. Main parameters, running time, accuracy and Kappa coefficient of simulated image.
AlgorithmMain ParameterValueTime (s)AccuracyHomogeneous Region
IIIIIIIV
The proposedspan parameter d301139.19user’s99.7699.4299.5098.67
producer’s10099.5798.1499.70
neighborhood coefficient η0.8overall99.34
Kappa0.99
ISODATAexpected number of clusters45.08user’s69.2952.8350.6777.41
minimum number of pixel2500producer’s99.6154.0136.1663.74
upper bounds of standar d deviation10overall63.34
lower limit of distance10Kappa0.51
AGENESspan parameter d3030.77user’s67.8653.1949.6279.08
producer’s99.9350.4640.2860.64
overall62.74
Kappa0.50
HMRF
FCM
fuzzy coefficient0.532.4449user’s98.3180.3295.1682.16
producer’s99.1996.0360.9994.79
neighborhood coefficient1overall87.76
Kappa0.84
Gamma
MRF
neighborhood coefficient0.25902.28user’s96.9893.0292.0192.06
producer’s98.7590.9688.4596.06
shape parameter4overall93.54
Kappa0.91
Table 3. Comparison of the estimated and actual values of mean and variance.
Table 3. Comparison of the estimated and actual values of mean and variance.
ParametersHomogeneous Region
IIIIIIIV
meanactual/estimated19.90/19.8779.24/79.39137.07/137.02208.95/208.45
deviation rate−0.00150.0019−0.0004−0.0024
varianceactual/estimated9.962/9.89239.092/39.07261.402/61.16257.532/57.922
deviation rate−0.0140−0.0010−0.00780.0136
Table 4. The number of clusters and overall accuracy of segmentation with different d and η.
Table 4. The number of clusters and overall accuracy of segmentation with different d and η.
dηApproximate Time (s)
Number of Clusters/Overall Accuracy
0.10.20.30.40.50.60.70.80.91.0
10334/94.734/97.064/97.334/97.454/97.394/97.994/97.56621,600
20334/93.614/96.054/96.604/97.664/98.185563000
30334/94.624/97.594/98.294/99.044/98.844/99.34661000
40334/93.384/93.714/97.784/98.314/98.724/98.6066800
50334/93.214/96.634/97.524/97.924/98.00555500
Table 5. The running time, number of clusters and overall accuracy of segmentation with different T.
Table 5. The running time, number of clusters and overall accuracy of segmentation with different T.
ParametersT
15203050
number of clusters/accuracy64/99.154/99.344/99.344/99.34
time (s)305.687435.1101139.191150.1071266.915
Table 6. Main parameters and running time of real images.
Table 6. Main parameters and running time of real images.
AlgorithmMain ParameterImage 1Image 2Image 3Image 4
ValueTime (s)ValueTime (s)ValueTime (s)ValueTime (s)
The proposedspan parameter d301387.64301450.60202099.92301192.03
neighborhood coefficient η0.50.40.70.5
inner loop T20202020
shape parameter4444
ISODATAexpected number of clusters22.9634.0133.3745.42
minimum number of pixel5000300030001000
upper bounds of standard deviation10201010
lower limit of distance10201010
AGENESspan parameter d303.48305.12309.17305.02
HMRF
FCM
fuzzy coefficient0.52.840.15.82860.56.170.58.60
neighborhood coefficient0.740.750.95
Gamma
MRF
neighborhood coefficient0.24226.350.210,593.520.410,367.170.349378.35
shape parameter4444

Share and Cite

MDPI and ACS Style

Zhao, Q.; Li, X.; Li, Y. Multilook SAR Image Segmentation with an Unknown Number of Clusters Using a Gamma Mixture Model and Hierarchical Clustering. Sensors 2017, 17, 1114. https://doi.org/10.3390/s17051114

AMA Style

Zhao Q, Li X, Li Y. Multilook SAR Image Segmentation with an Unknown Number of Clusters Using a Gamma Mixture Model and Hierarchical Clustering. Sensors. 2017; 17(5):1114. https://doi.org/10.3390/s17051114

Chicago/Turabian Style

Zhao, Quanhua, Xiaoli Li, and Yu Li. 2017. "Multilook SAR Image Segmentation with an Unknown Number of Clusters Using a Gamma Mixture Model and Hierarchical Clustering" Sensors 17, no. 5: 1114. https://doi.org/10.3390/s17051114

APA Style

Zhao, Q., Li, X., & Li, Y. (2017). Multilook SAR Image Segmentation with an Unknown Number of Clusters Using a Gamma Mixture Model and Hierarchical Clustering. Sensors, 17(5), 1114. https://doi.org/10.3390/s17051114

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop