COUPLING REGULAR TESSELLATION WITH RJMCMC ALGORITHM TO SEGMENT SAR IMAGE WITH UNKNOWN NUMBER OF CLASSES

This paper presents a Synthetic Aperture Radar (SAR) image segmentation approach with unknown number of classes, which is based on regular tessellation and Reversible Jump Markov Chain Monte Carlo (RJMCMC`) algorithm. First of all, an image domain is portioned into a set of blocks by regular tessellation. The image is modeled on the assumption that intensities of its pixels in each homogeneous region satisfy an identical and independent Gamma distribution. By Bayesian paradigm, the posterior distribution is obtained to build the region-based image segmentation model. Then, a RJMCMC algorithm is designed to simulate from the segmentation model to determine the number of homogeneous regions and segment the image. In order to further improve the segmentation accuracy, a refined operation is performed. To illustrate the feasibility and effectiveness of the proposed approach, two real SAR image is tested.


INTRODUCTION
Image segmentation is a hot topic in Synthetic Aperture Radar (SAR) image processing and involves two tasks: determining the number of homogeneous regions and segmenting them.By far, there are numerous algorithms for SAR image segmentation (Alkasassbeh et al., 2012;Ayed et al., 2005;Cao et al., 2005;Li et al., 2010;Schmitt et al., 2014;Xu et al., 2011;Xue et al., 2010).The statistics-and region-based approaches has been proved that they are much most effective for SAR image segmentation.To this end, geometry tessellation based segmentation approaches have been proposed, in which image domain is partitioned into a set of sub-regions and image model is built based on it.The commonly used geometric tessellations include regular tessellation (Askari et al., 2013), Voronoi tessellation (Lucarini, 2009;Zhao et al., 2013), Poisson tessellation (Schneider, 2009;Schneider, 2010), and so on.Li et al. (2010) proposed a region-based approach for SAR image segmentation based on Voronoi tessellation, Bayesian inference and Reversible Jump Markov Chain Monte Carlo (RJMCMC) algorithm.The proposed approach can segment SAR image well, but the number of homogeneous region need to be indicated by users.In order to automatically determine the number of homogeneous regions, Askari et al. (2013) presented an approach to SAR image segmentation.The approach is based on regular tessellation, Bayesian inference, a maximum likelihood Gamma distribution parameter estimator and RJMCMC algorithm.As the move types are designed unreasonably, the proposed approach can automatically determine the number of homogeneous regions, the segmentation accuracy is low and the proposed algorithm is time-consuming.The proposed approach couples regular tessellation with RJMCMC algorithm to automatically determine the number of homogeneous regions and segment the image.In the RJMCMC algorithm, five move types are designed, involving splitting or merging real classes, sampling parameter vector, sampling label field, birth or death of an empty class and splitting or merging blocks.These moves overcome the instability problem of segmentation optimization.This paper is organized as follows.In Section 2 we present the proposed algorithm.We then in Section 3 detail and discuss the results of real SAR images.Finally, Section 4 contains conclusions and perspectives for further research.

Image segmentation model
A SAR image, represented by a vector Y of pixel intensities, is defined on the image domain D, that is, Y = {Y d , dD}, where dD is a pixel location on D, Y d is a random variable defined on d.In this paper, the image domain D is partitioned into a set of blocks D = {P i , i = 1, …, m} by regular tessellation, where m is the unknown number of blocks and possess a prior distribution with the probability density function (pdf) p(m), the number of pixels in i th block P i is the multiple of 4 and the minimum block is 22. Figure 1 shows an example which D is partitioned into 6 rectangles by regular tessellation, that is, D = {P i , i =1, …, 6}.
A random label field Z = {Z i , i = 1, …, m},which is formed with the set of labels for all blocks, where each random variable Z i is associated with the block P i and takes its value in {1, …, k}, where k is unknown and has a prior distribution p(k).A realization of Z corresponds to a segmentation of the image.The improved stationary Potts model (Besag, 1986) is used to define the label field, which models the correlation of labels among blocks.In a given block P i , NP i ={ P i , P i  P i }denotes the block set of sites neighbouring block P i , i, i  {1, …, m} and i  i, P i  P i are neighbours if and only if they have a mutual boundary.The conditional pdf of Z i given Z i , P i  NP i , can be expressed as where t = 1 if Z i = Z i , otherwise t = 0, and  is a constant to control the neighbourhood dependences between a pair of neighbour blocks (Dryden et al., 2003).On the assumption that pdfs for Z i , i = 1, …, m, are conditionally independent, then the joint pdf of Z can be represented as It can be seen that the segmentation of Y can be fully described by the label field Z.
For the ith block P j , the pixels intensities Y i = {Y d , d P i } are modelled by a identical and independent Gamma distribution conditional on its label Z i = s, the joint pdf of Y i becomes where  () is Gamma function,  s = { s ,  s }.Assume that above joint pdfs are also independent, the joint pdf of Y = {Y i , i = 1, …, m} given Z, , k, m can be written as where parameter vector  = { s , s = 1, …, k} = { 1 , ,k, 1, , k}.Assume that the shape (scale) parameter s ( s ) satisfies identical and independent Gaussian distributions with mean   (  ) and standard deviations   (  ) (Yang et al.,  2006).Accordingly, the pdf of  can be expressed as The number of homogeneous regions k can be considered to have a Poisson distribution with mean  k (Green, 1995).k has a maximum equal to m (the maximum occurring if and only if when each block has a different class).Its prior distribution as follows The number of partitioned blocks in the image domain m is considered as a random variable, and has a prior truncated Poisson distribution with mean  m (truncated with minimum m = k, and the maximum m max = d if each block exactly has one pixel).Its prior distribution satisfies (Green, 1995) Segmentation is performed according to RJMCMC to simulate the posterior distribution defined in Eq. ( 8).Let  = (Z,, k, m).
In each iteration, a new candidate  * for  is drawn from a proposal distribution with density q(,  * ) (assume that the dimension of  * is higher than that of ).The random vector u is defined for accomplishing a transition from (, u) to  * with dimension matching, that is, || + |u| = |*|.The probability of accepting the candidate  * can be computed as (Green, 1995) where q(u) is u's density function and r() (r( * )) is a given move type probability in  ( * ).The Jacobian | ∂( * ) /∂(, u)| is the change of variable from (, u) to  * .The move types are designed as follow.
(1) Splitting or merging real classes In the proposed algorithm, the real class represents that at least one block is labeled the class and its number is k r .For splitting move, a set of blocks with a same class label randomly chosen in the current state k are split to two sub-sets with different class labels.As a result, the new label field is Z where  = ( Z, , k, m ) ,  * = ( Z * ,  * , k+1, m) and where k f r = f k /{k r (2 B -2)}, f k is the probability of choosing a split proposal, B is the number of blocks with label s, A merging is designed in tandem with a splitting so that they form a reversible pair (Green, 1995).So the acceptance probability of the merge move can be generalized as (2) Sampling parameter vector The parameter vector for Gamma distribution  is updated by sequentially updating the elements of  s for each label s  {1, …, k }.First of all, we need judge whether s is an empty.If s is a real class, the acceptance probability can be expressed as If l is an empty class, the second form of R p is used.
(3) Sampling label field For a block i , a new real label s * is randomly chosen from the set of current real classes and s * ≠ s (Z i = s).The acceptance probability can be written as  A death is designed in tandem with a birth (Green, 1995), so the acceptance probability can be written as (5) Splitting or merging blocks A move of splitting a block involves splitting a block to two blocks with different class labels.The operation includes: (i) Randomly choosing P i in the current partition of image domain D = {P 1 , …, P i , …, P m }, and its corresponding label Z i = s; (ii) If the pixel number of P i is more than 4 and the number of its row or column is the multiple of 2, P i can be split.Under the constraint of minimum blocks, splitting types can be determined.In order to choose the splitting type, the splitting types are labeled clockwise.As is shown in Figure 2, the block's row number is 6 and the column number is 10, they are both the multiple of 2, so there are 6 (= (6+10)/2-2) splitting types in the selected P i , the numbers of splitting types are labeled as 1, …., 6, respectively; (iii) One splitting is randomly chosen and splits P i into P i * and P m+1 * , and their corresponding labels are re-allocated s and s * , s  s * , respectively.The new partition of image domain becomes D * = { P 1 , …, P i * , …, P m , P m+1 * }, the acceptance probability of the splitting move can be expressed as, where s m or m m+1 is the probabilities of proposing a split or merge operation.The Jacobian term in Eq. ( 20) is equal to 1.A merging is designed in tandem with a splitting, so the acceptance probability of the merging move can be expressed as,

EXPERIMENTAL RESULTS AND DISCUSSION
The proposed approach is tested on real SAR images with sizes of 128128 as given in Figure 3, where Figure 3a presents a RADARSAT-I image with VV polarization and spatial resolution of 25m, and Figure 3b shows a RADARSAT-II standard mode image with HV polarization and spatial resolution of 25m.Visually, both of them include three homogeneous regions.
Experiment results demonstrate that (1) k r  [2, k max ] before a jumping period and after that k r attains to its stable value.In the experiments, k max =7.As is shown in Figure 4a and b, the iterations are 20 and 16, respectively; (2) during the jumping period, there is no class remaining unchanged with consecutively 10 iterations.Consequently, if k r is unchanged with 30 times continuously in iteration, it is assumed that be the number of the real class.As is shown in Figure 4, their real classes are equal to 3.After k r is determined, only move (2)-( 4) are preformed.Figure 5 shows real SAR image segmentation results.Figure 5a1 and b1 present the regular tessellation of coarse segmentation, and Figure 5a2-b2 shows the corresponding coarse segmentation.From Figure 5a1-b1 and a2-b2, it can be found that the accuracy of region segmentation is very high, but the boundary segmentation accuracy is poorer than it.In order to improve the accuracy of image boundary segmentation, a refined operation is performed.First of all, extracting the outlines of segmentation results and buffer zones with two pixels in width are generated centered on them.Then the pixels in the buffer zones are reclassified by the criterion defined in Eq. ( 1), new boundaries of homogeneous regions are tracked.Figure 5a3 and b3 show the refined segmentation result, respectively.From Figure 4 and 5, it is easy to find that the proposed approach can not only determine the number of classes but also segment the image precisely.
In order to qualitatively evaluate the proposed algorithm, the outlines of segmentation results (red) and refined outlines (green) are overlaid on the original image (see Figure 3) shown in Figure 6.Visually, the extracted outlines of segmentation results poorly match the real outlines, but the refined outlines match better.
In order to asses the segmented results shown in Figure 5a3-b3 quantitatively, some common measures, including Producer's accuracy, User's accuracy, Overall accuracy and Kappa coefficient, are computer based on the simulated real SAR images in Figure 7.
, and the new parameter vector is  * = { s * , s = 1, …, k, k+1}.It is evident that the splitting move does not affect m.The acceptance probability of the splitting move can be written as

*
Birth or death of an empty class In the proposed algorithm, the empty class means that no block is associated with the class and its number is denoted as k e .Given the current number of classes k, Let  = (Z,, k, m), b k or d k+1 is the probabilities of proposing a birth or death operation.For a birth move, first of all, a new empty class k +1 is proposed, then a new set  k+1 by identically and independently.So  * = (Z * , * , k+1, m), where  * = { s , s = 1, …, k, k +1}.The acceptance probability of a birth can be expressed as

Table 1 .
Table 1 lists them.As is shown in Table1, it can be seen that only three accuracies are below 90%.In addition, the kappa coefficients are up to 0.915.It can be illustrate that the proposed approach is feasible and effective.Quantitative evaluation