1 Introduction

Synthetic Aperture Radar images (SAR images) of the earth are an important tool for many scientific applications such as high resolution remote sensing for mapping, surface surveillance, search-and-rescue, mine detection, navigation position location, and Automatic Target Recognition. For most of the applications mentioned, besides the importance of segmentation quality in the subsequent analysis for target detection and recognition, the time computing can play a key role. This constraint is difficult to reach because of the large size of SAR images. Amongst several existing works especially designed for SAR data [11, 19], one finds the contribution of the unsupervised statistical approach [17, 20].

Tree main statistical approaches have been suggested depending on the way that neighborhood influences the classification of a given pixel. The local or blind approaches consider that pixels are spatially independent [18]. Contextual methods [14] take into account a neighborhood of limit extent. Global approaches [17] assume that all pixels in an image influence the classification of the pixel of interest. Different estimation algorithms have been applied to the problem of unsupervised segmentation. The most popular one is the expectation-maximization (EM) algorithm [1, 6], with the maximum likelihood (ML) estimation. Many variants of the EM have been presented such as Stochastic EM (SEM) [14], Gibsien EM (GEM) [5], Iterative Conditional Estimation (ICE) [2, 4] and Zhang algorithms [21]. Because of the inherent speckle noise, it is now accepted that the statistics of SAR images can be well modeled by the family of probability distributions known as the Pearson system [8, 12].

We base our work in the field of statistical approaches; the mixture identification will be performed with the generalized mixture expectation maximization algorithm GMEM where the generalized aspect comes from the use of the Pearson system.

In statistical segmentation, the time increases with the size of training data set. Because of the large size of SAR images, the size of the training data is very large. As a result, the time required to segmentation could be prohibitively large which constrain its use in the real-world applications.

In this context, we propose a fast segmentation algorithm based on the principle of Bootstrap sampling. The Bootstrap [9] has not seen much applications in image analysis and classification. The first attempt to apply the bootstrap to image analysis was by Ghorbel and Banga [10] who introduced a Bootstrap scheme in the context of Bayesian image segmentation under the Gaussian assumption. M’hiri et al. [15, 16] extended the technique later and applied them to the segmentation of brain images.

In this paper, we propose to combine the bootstrap technique with the generalized mixture expectation maximization algorithm GMEM for unsupervised SAR image classification. The Bootstrap technique consists on selecting a small representative set a pixels from the original image and mixture identification will be done on bootstrap samples instead of the correlated pixels in the real image. The use of such resampling procedure enables to reduce considerably the computation times while preserving the estimation equality. The reminder of this paper is organized as follow. In Sect. 2, we describe the unsupervised Bayesian segmentation by Bootstrapping GMEM algorithm, Sect. 3 is devoted to experimental results of SAR image classification. Conclusions and future prospects are presented in Sect. 4.

2 Unsupervised Image Segmentation by the BGMEM Algorithm

2.1 Bootsrap Sampling for Image Classification

Efron [9] introduces the Bootstrap term to designate the set of the random re-sampling procedures of the data observed intended to be approached by simulating the statistics of the underlying distribution. The Bootstrap theory is based on the convergence of the empirical law of the sample towards the underlying unknown law when the sample size is sufficiently large.

Given a random sample \(\chi _N=({X_1,X_2,..,X_N}\)) of size N from a population with distribution \(F_X\). The bootstrap approximation was to estimate the distribution of a given unknown statistics \(R(\chi _N,F_X)\) by the bootstrap distribution noted \(F_n^*\) corresponding to the sample \(\chi _N^*={X_1^*,X_2^*,..,X_n^*},(n<N)\) where \( X_1^*,X_2^*,..,X_N^*\) are randomly selected from \(\chi _N\). Since the empirical distribution converges almost surely to the underlying distribution, one can hope that the bootstrap distribution would converge to the true unknown distribution. Details and applications of bootstrap technique can be found in Ref. [3, 22].

Ghorbel and Banga [10] introduce a bootstrap sampling scheme for gray levels image analysis. In this case, the whole image distribution is presented only by a small bootstrap sample size but should be representative of the entire image.

Suppose we have a two-dimensional image of \(N_r*N_c=N\) pixel resolution with \(N_r\) rows and \(N_c\) columns. In statistical segmentation, we suppose that the image is a finished population of N observations. It is then noted \(Y=({y_1, y_2,...,y_N})\) that is to say a sample from an unknown distribution. From this original image, we construct the bootstrap sample \(Y^*=(y_1^*,...,y_n^*)\), \((n<N)\) by randomly selecting n pixels \(y_i^*=(k_i,l_i)^t\) for \(i\in \left\{ 1,..,N\right\} \). \(k_i\) and \(l_i\) are obtained by making independent uniform random trials under the set \(\left\{ 1,..,N_r\right\} \) and the set \(\left\{ 1,..,N_c\right\} \) respectively. The \(Y^*=(y_1^*,...,y_n^*)\) is a resample of size n chosen with replacement from Y.

The optimal size of the bootstrap sample is determined by two representative criteria for the use of this technique in image segmentation.

The sample \(Y^*= (y_1^*,...,y_n^*)\) is representative of an image Y when each image gray level (GL) appears at least one time in the samples. It’s explicitly shown by the C1 and C2 equations.

$$\begin{aligned} C1: n_0 > 4D \,\,\,with \,\,\,D: \,\,\,the \,\,\,variation\,\,\, of \,\,\,the\,\,\, gray\,\,\, level \end{aligned}$$
(1)
$$\begin{aligned} C2: n_0 \,\,\,is\,\,\, representative\,\,\, \Longleftrightarrow B(n_0)=\displaystyle {\sum _{j=1}^D}\frac{\pi _j e^{-n_0 \pi _j}}{1-e^{-n_0 \pi _j}}<\epsilon \quad \\\nonumber \end{aligned}$$
(2)

D is the number of the different gray levels values in the image, \(B(n_0)\) is the sampling characteristic function and \(\epsilon \) is a fixed small value.

The advantage of the C1 criterion is to be able to depart from a minimum sample size \(n_0\). The sampling characteristic function \(B(n_0)\) defined by the criterion C2 takes into account both image pixel distribution and bootstrap sample size . The \(n_0\) is initially computed by the C1 criterion then \(n_0\) is progressively increased. The corresponding bootstrap sample is constructed for each size \(n_0\) by random uniform selecting set of pixels from the image and the C2 criterion is verified . The bootstrap sample is qualify to be representative when the sampling characteristic function value is lower than the expected precision \(\epsilon =10^{-2}\).

An empirical convergence study based on minimizing the mean integrated square error (MISE) between original density and its estimate based on the bootstrap sample show that the fixed precision \(\epsilon =10^{-2}\) ensure the representativity of the bootstrap sample.

2.2 An Overview of the Pearson System of Distributions

The Pearson system [13] is made up of mainly height families of distributions including Gaussian, Gamma and Beta ones and offers a large variety of shapes (symmetrical or not, with finite or semi-finite support, etc.). Each law can be uniquely defined by its mean \(\mu _1\) and its first three centered moments \((\mu _2, \mu _3, \mu _4)\). All of them can be represented in the so-called Pearson diagram (Fig. 1) in which axes \(\beta _1\) and \(\beta _2\) are given by \(\beta _1=\frac{(\mu _3)^2}{(\mu _2)^3}\) and \(\beta _1=\frac{(\mu _4)}{(\mu _2)^3}\).

Gaussian distributions are located at \((\beta _1=0,\beta _2=3)\). Gamma distributions on the straight line \(\beta _2=1.5 \beta _1+3\) and inverse gamma distribution on the curve with the equation \(\beta _2=\frac{3}{\beta _1-32}(-13\beta _1-16-2(\beta _1+4)^\frac{3}{2})\). First kind Beta distributions are located between the lower limit and the gamma line, second kind Beta distributions are located between the gamma and the inverse Gamma distributions, and Type 4 distributions are located between the inverse Gamma distributions and the upper limit. Then it is possible to estimate empirical moments of a distribution from a sample and to assess the family of distributions from coordinates \((\beta _1,\beta _2)\) and determine the parameters that precisely characterize the probability density function within its family.

Fig. 1.
figure 1

Eight families of the Pearson system in the \((\beta _1,\beta _2)\)-diagram.

Full size image

2.3 Description of the Bootstrapped Generalized Expectation Maximization Algorithm

Before proceeding with mixture identification by the BGMEM algorithm, a pre-processing step is needed to determine the size of the representative bootstrap sample because the performance of the mixture identification depends on the pre-processing step. The sample size drawn from the image fulfills the criteria proposed by Ghorbel and Banga in Sect. 2.1. Let given an available bootstrap sample \(Y^*=(y_1^*,...,y_n^*)\) from the original image Y.

The mixture identification by BGMEM algorithm is performed using only the representative Bootstrap sample instead of using the entire image. After initializing parameters from the histogram, the two following steps are iterated.

  • Expectation Step: For each class, the distribution \(f_k\) will be selected from the Pearson system distribution according to the skewness and kurtosis values [7] then the a posterior probability for a pixel \(y_i^*\) to belong to class k at the iteration q is given by:

    $$\begin{aligned} \forall k\in \{1,...,K\} \,\,\,\,\,\,\ P(x_k|y_i^*,\theta ^q)= \frac{\displaystyle {\pi _k^{q-1}f_k(y_i^*|\theta _k^{(q-1)})}}{\displaystyle \sum _{l=1}^K \pi _l^{q-1}f_l(y_i^*|\theta _l^{q-1})} \ \end{aligned}$$
    (3)

  • Maximization Step: the parameters of the mixture are constructed in the following way:

$$\begin{aligned} \forall k\in \{1,...,K\} \,\,\,\,\,\,\pi _k^{(q)}= \frac{\displaystyle \sum _{i=1}^n P(x_k|y_i^*,\theta ^q)}{\displaystyle n} \end{aligned}$$
(4)
$$\begin{aligned} \forall k\in \{1,...,K\} \,\,\,\,\,\,\mu _{1,k}^{(q)}= \frac{\displaystyle {\sum _{i=1}^n P(x_k|y_i^*,\theta ^q)}y_i^*}{\displaystyle \sum _{i=1}^n P(x_k|y_i^*,\theta ^q)} \end{aligned}$$
(5)

Hight order moments \(\mu _{j}\) with \(j={2,3,4}\) are computed according the equation:

$$\begin{aligned} \forall k\in \{1,...,K\} \,\,\,\,\,\,\mu _{j,k}^{(q)}= \frac{\displaystyle {\sum _{i=1}^n P(x_k|y_i^*,\theta ^q)(y_i^*-\mu _{1,k}^{(q)})^j}}{\displaystyle \sum _{i=1}^n P(x_k|y_i^*,\theta ^q)} \end{aligned}$$
(6)

At the \(q^th\) iteration on the sample the skewness and the kurtosis of the class k are computed as:

$$\begin{aligned} \beta _{1,k}^{(q)}=\frac{ {(\mu _{3,k}^{(q)})}^{2}}{{(\mu _{2,k}^{(q)})}^{3}} \,, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\beta _{2,k}^{(q)}=\frac{ {\mu _{4,k}^{(q)}}}{{(\mu _{2,k}^{(q)})}^{2}} \ \end{aligned}$$
(7)

The algorithm applied to the sample stops when the sample estimated parameters are stagnated.

3 Unsupervised Bayesian Classification of SAR Images

In the unsupervised Bayesian classification, the image is supposed as a realization of mixture distributions and the classification problem consists on mixture identification associated to a Bayesian rule. A large study is done to validate the BGMEM framework on a large data set of both synthetic images and real SAR images. In this section, we will show results only for four SAR images. The original SAR images are \(512*512\) pixels resolutions presented by 256gray levels classified in four classes. Before proceeding with image classification by the BGMEM algorithm, a pre-processing step is needed to determine the size of the representative bootstrap sample according to the representativity criteria presented in Sect. 2.1. The bootstrap sample size must be large enough to ensure a good estimation and small enough to reduce the computation time. The sample size \(n_0=3000\) pixels is validate according to the sampling characteristic function equation \((B(3000)<0.004)\) for the four images.

Table 1. Evaluation of the MSE measure of the mixtures
Full size table
Table 2. Time segmentation in seconds of the SAR images
Full size table

The Fig. 2a shows the fitting of the estimated densities to the image histogram in the case of classical GMEM algorithm for the four images and the Fig. 2b shows the result in the Bootstapped case. The conditional estimated densities fit the image histogram in the two cases of classical and bootstrapped algorithm.

We present in the Table 1 a comparison of the Mean Square Error MSE between the classical GMEM algorithm and the bootstrapped one. The mean square error shows an improvement of estimation when considering the bootstrapped version of the algorithm. The observations are decorrelated by randomly selecting observations which construct the sample. This fact offer best conditions of application of the GMEM algorithm and so better parameter estimation.

Fig. 2.
figure 2

Plot of estimated distributions by (a) GMEM algorithm (b) BGMEM algorithm

Full size image
Fig. 3.
figure 3

(a) Original images (b) segmented image from classical GMEM algorithm (c) segmented image from BGMEM algorithm.

Full size image

The Fig. 3b and c represent the segmentation results we obtained by applying the classical GMEM and the Bootstrapped BGMEM algorithms. We can observe visually a same segmentation quality by the both algorithms which is consequently of good mixture identification. We dont focus our work in the interpretation of the SAR images so the segmentation quality is showed visually but the contribution here was in accelerating the time computing of image classification. The computational time of the Boostrap algorithm and the classical one are given in Table 2. The programs are done with matlab and turned under a core i5 processor time computing may be lower if the algorithms are running on a more powerful processor. On the basis of the bootstrap sample which represents 0.01 of the image size, we have considerable gain in time computing more than a factor of 87 by the boostrapped algorithm with the same accuracy and robustness of image classification.

4 Conclusion

In this paper, we have proposed a Bootstrap model to the unsupervised Bayesian image segmentation. The Bootstrap technique allows an estimation of parameters of the image from a small sized sample instead of the entire image. The principle advantage of the proposed algorithm is the reduction of time computing which make it useful in real-time applications. As future work, we are interest to the study of multivariate Pearson system for the statistical modeling of images acquired by different sensors.