Irregular Tessellation and Statistical Modeling Based Regionalized Segmentation for SAR Intensity Image

: This paper presents a regionalized segmentation method for synthetic aperture radar (SAR) intensity images based on tessellation with irregular polygons. In the proposed method, the image domain is partitioned into a collection of irregular polygons, which are constructed using sets of nodes and are used to fit homogeneous regions with arbitrary shapes. Each partitioned polygon is taken as the basic processing unit. Assuming the intensities of the pixels in the polygon follow an independent and identical gamma distribution, the likelihood of the image intensities is modeled. After defining the prior distributions of the tessellation and the parameters for the likelihood model, a posterior probability model can be built based on the Bayes theorem as a segmentation model. To obtain optimal segmentation, a reversible-jump Markov chain Monte Carlo (RJMCMC) algorithm is designed to simulate from the segmentation model, where the move operations include updating the gamma distribution parameter, updating labels, moving nodes, merging polygons, splitting polygons, adding nodes, and deleting nodes. Experiments were carried out on synthetic and real SAR intensity images using the proposed method while the regular and Voronoi tessellation-based methods were also preformed for comparison. Our results show the proposed method overcomes some intrinsic limitations of current segmentation methods and is able to generate good results for homogeneous regions with different shapes.


Introduction
The synthetic aperture radar (SAR) system has been applied to numerous fields, providing high-resolution imagery, having a strong penetration capability, and acquiring images both day and night and in all weather conditions [1,2]. In SAR image processing, segmentation is the most important task, because its result is critical for subsequent high-level image processing tasks, such as feature extraction, object recognition, and classification [3,4]. Currently, many SAR image segmentation methods have been proposed, including ones based on region [5,6], edge [7,8], clustering [9,10], and statistical models [11,12]. However, the unique imaging mechanism of the SAR system produces images that contain substantial inherent speckle noise. The results of SAR image segmentation are highly susceptible to the noise. Many studies have shown that the speckle noise inherent in SAR imagery can be characterized as a statistical distribution; thus, the statistical model-based segmentation methods are considered to be more efficient [13,14]. Traditionally, statistical segmentation methods always take the pixel as the primary processing unit, resulting in heavy missegmentation [15,16]. In order to solve this issue, region-based statistical segmentation methods have been proposed [17,18]. Their basic principle is that the image domain is first partitioned into a group of sub-regions, and the partitioned sub-regions are viewed as basic processing units to model image segmentation. Wang et al. [19] utilized regular tessellation to partition the image domain into rectangles or squares. Then, several sub-regions were used to fit a homogeneous region. The segmentation method is straightforward and tackles the issue of error segmentation in homogeneous regions to a certain extent; however, the effect of boundary fitting can be vague and irregular. In order to improve boundary fitting, Zhao et al. [20] proposed using Voronoi tessellation instead of regular tessellation. Although the effect of boundary fitting is slightly improved, there are still difficulties in fitting specific shapes due to the convex characteristics of Voronoi polygons. The proposed method is based on an iterative elaboration of the selected polygons, in order to best fit the shapes in the geometry of the intensity image.

Segmentation Model
Assume that a SAR image contains a known number, k, of homogeneous regions, and its domain is partitioned into an unknown number, m, of polygons by irregular tessellation, where m can be considered as a random variable and follows a prior distribution. Each polygon assigns a random variable, Lj, to indicate its membership to a homogeneous region. The label variables for all polygons form a label field, L = {Lj; j = 1, ..., m}. Let l = {lj ∈{1, ..., k}; j = 1, ...., m} denote a realization of L, which corresponds to a segmentation result of the image.
According to the Bayesian theorem [21], a posterior probability as the segmentation model can be written as where p(Z | m, L, P, β) is the likelihood model for the SAR image, characterized by gamma distribution [22], where α and β = {βl; l = 1, ..., k} are the shape and scale parameters of the gamma distribution, and Δl = {Dj: lj = l, j = 1, ..., m} is the collection of polygons with the same label l. p(L| m) is the statistic model of the label field L, which can be defined based on the Markov random field (MRF) [23,24]  where c > 0 is a constant that controls the neighborhood dependence between pairs of neighboring polygons [25] and δ is an indicator function, ( ) where np indicates the total number of images.
For a multi-look SAR intensity image, the shape parameters α of its gamma distribution can be viewed as the number of its looks [26]. p(β) is the joint distribution of {βl; l = 1, ..., k}. Under the assumption that βl, l = 1, ..., k, follow identical independent Gaussian distributions with the same distribution parameter μβ and σβ , p(β) can be expressed as

Simulation
In order to obtain the optimal segmentation from Equation (8), a reversible-jump Markov chain Monte Carlo (RJMCMC) [28] algorithm was designed to simulate from the posterior probability and the maximum a posterior (MAP) [29] scheme was adopted. Using MAP, the optimal segmentation result could be expressed as To obtain the optimal segmentation result, seven move operations were designed. Move 1: Updating the gamma distribution parameter. Randomly select a scale parameter from β = {βl; l = 1, ..., k}, for example βl. The proposed βl * is drawn from a Gaussian distribution with the mean βl and standard deviation εβ (i.e., βl *~ N(βl, εβ)). The acceptance probability of changing β to β * = {β1, ..., βl-1, βl * , βl+1, ..., βk} can be calculated as, Move 2: Updating the label. A polygon Dj with the label lj is uniformly drawn with the probability 1/m from {Dj; j = 1, …, m}. To update its label, a proposed label lj * is uniformly drawn from {1, ..., k} and restricted to lj * ≠ lj. The realization l = {l1, ..., lj, ..., lm} of L after updating lj to lj * becomes l * = {l1, ..., lj-1, lj * , lj+1, ..., lm}. The acceptance probability of changing l to l * can be written as Move 3: Moving the node. A node is uniformly selected from mobile nodes Pmm (e.g., (sjg, tjg)) to move to a new proposed candidate point (sjg * , tjg * ). If the node is located inside the image domain, the candidate point is selected randomly within a circle with the center at (sjg, tjg) and the radius r. If the node is on the boundary of the image, it can only move along the side of the boundary that it is on. If the node moves from (sjg, tjg) to (sjg * , tjg * ), the polygon Dj is changed into Dj * . Consequently, the set of the nodes for the polygon j, Pj = {(sj1, tj1), …, (sjg, tjg), …, (sjn j , tjn j )}, become Pj * = {(sj1, tj1), …, (sjg * , tjg * ), …, (sjn j , tjn j )}. The node sets for all polygons become P * = P1∪...∪Pj-1∪Pj * ∪Pj+1∪...∪Pm.  The acceptance probability for the moving node can be written as  The acceptance probability of merging polygons can be written as, Move 5: Splitting polygons. A polygon with more than three nodes is selected randomly, for example Dj. A node (sjg, tjg) is selected from the node set Pj, randomly. Then another node is selected, for example (sjg′, tjg′), which is not adjacent to (sjg, tjg) from the node set Pj. The two selected nodes, (sjg, tjg) and (sjg′, tjg′), are then connected. As a result, Dj is split into two polygons, denoted as Dj *   The acceptance probability of splitting polygons can be written as Move 6: Adding a node. For a randomly selected polygon, for example Dj, with one of its edges inside the image domain with the end points (sjg, tjg) and (sjg+1, tjg+1), a candidate node, denoted as (sjn j +1, tjn j +1), is uniformly drawn from the area from the intersection of the image domain and the circle with the center at the edge and the diameter equal to the length of the edge. Consequently, the set of nodes Pj = {(sjg, tjg) ∈ D; g = 1, …, nj} is changed to Pj * = {(sjg, tjg) ∈ D; g = 1, …, nj+1}, and Pj' = {(sj'g, tj'g) ∈ D; g = 1, …, nj'} is changed to Pj' * = {(sj'g, tj'g) ∈ D; g = 1, …, nj'+1}. Additionally, the polygon Dj and its neighbor Dj', which share the edge with Dj, are changed to Dj * and Dj' * , respectively. Figure  5 shows an example of adding a node, where the edge selected is made up of points (s25, t25) and (s24, t24) from D2, and its neighbor polygon sharing the edge is D3. After the operation, the number of nodes in polygon D2 is transformed from five to six, while the number of nodes in polygon D3 is also changed from five to six. Polygons D2 and D3 are changed into D2 * and D3 * , respectively. The acceptance probability of adding a node can be written as Move 7: Deleting a node. A concave polygon Dj is selected randomly from the current m polygons. Two adjacent edges are selected from Dj. The common node of the two edges is deleted. After deleting, the new edge can be generated by connecting the remaining nodes of the edges.

Experimental Results and Discussion
To verify the feasibility and effectiveness of the proposed algorithm, experiments were performed on both the simulated and intensity SAR images.

Experiment Setups
The proposed method was implemented by MATLAB programming and run on an Intel(R) Core (TM) i7-4790 CPU@ 3.6GHz 4G computer. The running time was about 30 min. Table 1 lists the constants used in the proposed algorithm for the experiment. As the number of polygons changed by splitting and merging polygon operations, the value of the initial m had no significant impact on the results. In order to facilitate the following move operation, the initial tessellation partitioned the image domain D into 16 squares, such that the initial m = 16. The constant c was the coefficient used for characterizing the dependence of the neighboring polygons in the improved Potts model as defined in Equation (3), such that the conditional probability of the label for a polygon was a monotonic function. Based on a series of experiments, c ∈ [0.5, 1.5] was recommended [30]. In the experiment, the constant c was set to one. For multi-look SAR images, the shape parameter α of the gamma distribution was equal to the number of its looks. In the experiment, four-look SAR images were used for testing the proposed algorithm, such that α = 4. The initial scale parameters of gamma distribution were also drawn from their distributions, i.e., βl ~ N(μβ, σβ), l = 1, ..., k. The values of the mean μβ and standard deviation σβ were set as 32 and 8, respectively. The constant εβ was the proposed variance for β, which affects the convergence rate of the algorithm under the RJMCMC scheme. When εβ is set to be too large, the parameter estimation is inaccurate, but if εβ is set to be too small, the algorithm converges slowly. According to [31], the proper parameters can be obtained when εβ is set to one. In the experiment, the number of iterations was set to 10,000. Generally, 10,000 iterations are sufficient to ensure convergence of the algorithm.

Simulated SAR Image Segmentation
To create a simulated image, a template of size 128 × 128 was constructed, as shown in Figure  7a. The number of homogeneous regions was three. For the simulated SAR image, each homogeneous region followed a gamma distribution, where the shape parameter was equal to four, and the scale parameter is listed in Table 2. The simulated image is shown in Figure 7b. In order to verify the capability of fitting homogeneous regions with sharp angles in the proposed method, the simulated image was designed to have a sharp angle, as shown in Figure 7b highlighted by a blue rectangle. Figure 7c shows the zoomed image of the sharp angle area.   Figure 8 displays the segmentation results of the simulated image. Figure 8a shows the partition result from the irregular tessellation, and Figure 8b provides the optimal segmentation result. As shown in Figure 8a, each homogeneous region was fitted by one irregular polygon. To evaluate the proposed segmentation method qualitatively, the outlines of the segmented homogeneous regions (marked in red) were overlaid on the original image, as shown in Figure 8c. In order to verify the effectiveness of the proposed irregular partition method, the regular partition method [19] and Voronoi partition method [31] were used as comparison method. The results from the regular tessellation-based and Voronoi tessellation-based segmentation methods are illustrated in Figures 9 and 10. The partition results of the regular and Voronoi tessellations are presented in Figures 9a and Figure 10a, the optimal segmentation results are provided in Figures 9b  and 10b, and the visual evaluation results are presented in Figures 9c and 10c.  Figure 9, regular tessellation solved the problem of error segmentation caused by speckle noise in homogeneous regions to a certain extent, but boundary fitting results were coarse and uneven. From Figure 10, although the Voronoi polygons fitted most of the segments along the boundaries well, the results show that the method had difficulty matching the region with sharp angles (see Figure 10f).  In order to verify the effectiveness of the proposed method for SAR image segmentation, the most classic multi-threshold method [32] was used as the contrast method to carry out experiments on the simulated SAR image, and the experimental results are shown in Figure 11. Figure 11a is the segmentation result. It can be seen that this method segmented region I better, but there were still many noises in regions II and III. In order to evaluate the experimental results qualitatively from the visual vision, the contour lines of the segmentation results were overlapped with the experimental image, as shown in the Figure 11b. It can be clearly seen that there was a lot of noise in regions II and III. To evaluate quantitatively the effectiveness of all the segmentation methods, the confusion matrixes for their segmentation results (Figure 8b, Figure 9b, Figure 10b, and Figure 11a) were calculated, taking the image from Figure 7a as a reference. Based on the matrixes, the producer's accuracy, user's accuracy, the overall accuracy, and the Kappa coefficient were calculated and listed in Table 3 [33]. From Table 3, the overall precision of the proposed segmentation method was 98.57%, which demonstrates that the proposed method achieved excellent segmentation and was comparatively better than the regular, Voronoi tessellation-based segmentation methods and the multi-threshold segmentation method. Table 3. Quantitative evaluation of accuracy for the proposed method, the regular tessellation-based method, the Voronoi tessellation-based method, and the multi-threshold segmentation method.  Figure 12 shows the changes of scale parameters throughout the 10,000 iterations for the different segmentation methods. Figure 12a-c shows the parameter variations of the proposed, the regular, and the Voronoi tessellation-based segmentation methods, respectively. From Figure 12a,c, the parameter for each homogeneous region converged to some steady value quickly (after about 1000 iterations). Figure 12b shows that the parameters were convergent but were not steady.   Table 4 provides the estimation values of the scale parameters and their percentage errors. It can be seen from Table 4 that the estimated parameters of the proposed and the Voronoi tessellation-based segmentation methods were more accurate than that of the regular tessellation-based segmentation method, particularly for region I. Figure 13 illustrates the histograms of the different homogeneous regions in the simulated image and the gamma distributions with real parameters. The histograms and the distribution function curves for the various homogeneous regions all matched accordingly.  Figure 14 presents four multi-look SAR intensity images with a size of 128 × 128 pixels, a spatial resolution of 30 m, and the number of looks equal to four. Figure 14a is a RADARSAT-1 image with HH polarization of a coastal scene, where the darker region indicates water (C1), and the brighter region represents the coast (C2). Figure 14b is a RADARSAT-II estuary image with HV polarization, where the darker regions indicate riverbanks (C1), and the middle portion marks ice (C2). Figure 14c is a RADARSAT-I image of sea ice with HH polarization. It contains three sea-ice structures, which from light to dark represent sea ice (C2), partially melted ice (C1), and water (C3), respectively.

Real SAR Intensity Images
(c) Figure 14d shows a RADARSAT-I standard-mode image with HH polarization, which shows four types of sea-ice structures, representative of different degrees of ice (light) and water (dark). The proposed segmentation method was applied to the four SAR intensity images shown in Figure 14 using the same constants listed in Table 1. The segmentation results are shown in Figure  15. Figure 15a shows the result of dividing into irregular polygons. Figure 15b shows the segmentation results of the test images with varying numbers of homogeneous regions. In general, the number of irregular polygons should be the same as the number of homogeneous regions of the image. However, for non-adjacent homogeneous regions, the number of polygons was different from the number of homogeneous regions. As shown in the results of Figure 14b, the number of irregular polygons was three, while the number of homogeneous regions was two. This means that there were two polygons with the same label. In Figure 15c, the outlines of the segmented homogeneous regions were extracted and overlaid over the real images. From a visual examination, the proposed segmentation method was able to attain the optimal segmentation results.  Figure 16 shows the variations of the scale parameters throughout the iteration process. The abscissa and ordinate represent the number of iterations and the scale parameter value, respectively. It can be seen from Figure 14 that the proposed segmentation method made the parameter value converge quickly and get stable parameter values. To compare the results of the proposed segmentation method, the real images shown in Figure  14 were tested using the regular and Voronoi tessellation-based segmentation methods; their segmentation results are shown in Figures 17 and 18. There were some noises in the segmentation results of the regular tessellation method, as shown in Figure 17b. In Figure 18, the Voronoi tessellation-based segmentation method did not produce misclassified noise, but the boundary fitting was still not sufficiently accurate. Because the Voronoi polygons were convex, there were some constraints when using them to fit homogeneous regions. Figure 19 shows the experimental results of the multi-threshold segmentation method, where Figure 19a is the segmentation result and Figure 19b is the result of the overlap of the contour of the segmentation result and experimental image. It can be seen from the experimental results that this method had a good segmentation effect on the small difference in the pixel gray value (as shown in Figure 14b), but it was difficult to achieve a better segmentation result on the large difference in the pixel gray value.

Conclusion
This paper presents a new regional-based segmentation method for SAR intensity images. In this method, the image domain is partitioned into a collection of irregular polygons, which are constructed with a set of nodes and are used to fit homogeneous regions with arbitrary shapes. The commonly used regional-based segmentation method is based on either regular or Voronoi tessellations. The regular and Voronoi polygons can only fit in convex homogeneous regions but correspond poorly in areas with sharp corners. Compared with widely used segmentation methods, the proposed segmentation method provides greater flexibility, which can be used in both convex and concave homogeneous sections. As the nodes are used in segmenting sections, fewer polygons are needed to fit the homogeneous regions, resulting in more well-defined segments compared with other tessellation-based segmentation methods. Combining irregular tessellation and Bayesian inference to segment SAR images, the experiment results show that the proposed segmentation method can provide better segmentation results and, at the same time, the model parameters can be estimated more accurately. The proposed irregular tessellation method provides a new idea for region-based image segmentation and a better solution for region segmentation with different geometric shapes. However, it has no effect on the image category (sea ice, city, and forest). This paper only takes low-resolution SAR sea-ice images as examples for simple verification. In future work, high-resolution urban and forest images can be considered for segmentation experiments to verify the universality of the proposed method.