Parallel Regional Segmentation Method of High-Resolution Remote Sensing Image Based on Minimum Spanning Tree

: With ﬁner spatial scale, high-resolution images provide complex, spatial, and massive information on the earth’s surface, which brings new challenges to remote sensing segmentation methods. In view of these challenges, ﬁnding a more e ﬀ ective segmentation model and parallel processing method is crucial to improve the segmentation accuracy and process e ﬃ ciency of large-scale high-resolution images. To this end, this study proposed a minimum spanning tree (MST) model integrated into a regional-based parallel segmentation method. First, an image was decomposed into several blocks by regular tessellation. The corresponding homogeneous regions were obtained using the minimum heterogeneity rule (MHR) partitioning technique in a multicore parallel processing mode, and the initial segmentation results were obtained by the parallel block merging method. On this basis, a regionalized fuzzy c-means (FCM) method based on master-slave parallel mode was proposed to achieve fast and optimal segmentation. The proposed segmentation approach was tested on high-resolution images. The results from the qualitative assessment, quantitative evaluation, and parallel analysis veriﬁed the feasibility and validity of the proposed method. the parallel RHMRF-FCM algorithm.


Introduction
Remote sensing image segmentation is a process of partitioning image domain into several meaningful regions and serves as a bridge between remote sensing image and high-level information processing [1,2]. Remote sensing image is a type of data that records spectral information of the earth's surface from a spaceborne sensor. It is a standardized spatial data that depicts the spatial grid (pixel) distribution of spectral information after the surface has reflected solar electromagnetic waves [3,4]. In recent years, with more satellites being sent to space, coupled with improvements in sensor technology, high-resolution remote sensing data has gradually became ubiquitous [5]. With a finer spatial scale, high-resolution images provide more detailed information through abundant geometric, textural, and spectral features, resulting in new challenges to the traditional remote sensing image segmentation methods [6,7]. These challenges include geometric noise caused by peripheral objects, the phenomena of different objects having the same spectral signature, and similar object types shown as having varying spectral signatures in high-resolution images [8]. Improvements in spatial resolution have been accompanied by the rapid growth in data scale. Large-scale high-resolution remote sensed imagery is characterized by the use of massive data, which results in additional difficulty for traditional serial methods to meet speed requirements in data processing [9,10]. Traditional segmentation algorithms, based on spectral feature clustering, are unable to effectively segment where • stands for ceil operation, g × g is the tessellation parameter.

Block Partitioning
For subimage Z t , its image domain F t can be divided into I t S disjoint nonempty regions according to the MHR, such that, F = {F ti ; i = 1, 2, . . . , I t S }, the corresponding sub-image is Z ti = {z a ; a = 1, 2, . . . , N ti , (u a , v a ) ∈ F ti }, N ti is the number of pixels in the i-th region of block t, a is the local pixel index of the i-th region of block t, S is the threshold of MHR, and each homogeneous region F ti corresponds to a certain land use and land cover (LULC) class or part of it. In order to better express the explicit and implicit features of blocks, the MST model is used to express blocks as image minimum spanning tree field (IMSTF). Let the pixels of block be mapped to vertices, the 8-neighborhood of pixels is taken as the adjacency relationship, and the nonsimilarity measure between the spectral measure vectors between the pixels is taken as the edge weight, the block can be expressed as graph Remote Sens. 2020, 12, 783 3 of 30 The proposed method extended the hidden Markov random field (HMRF) from modeling the neighborhood relationship of pixels to the neighborhood relationship of regions, which more effectively obtained model parameter estimation and optimal image segmentation.

Methods
Given a high-resolution image Z = {zn = Z(un, vn); (un, vn) ∈ F, n =1, 2, …, N, un ∈ [1, row], vn ∈ [1, col]} defined in the image domain F, where (un, vn) is the coordinate of the pixel n, n is the pixel index, N is the number of pixels, row and col is the number of rows and columns of the image, respectively, N = row ×col, Z(un, vn) is the spatial sampling function, zn = {znb; b = 1, 2, …, B} represents the spectral measure vector of pixel n, b is the band index, and B is the number of bands. Suppose that Z consists of K homogeneous regions, such that, F = {Dk; k = 1, 2, …, K}, where k is the index of homogeneous region, and Dk is homogeneous region K.
where Tr is the rows of blocks, and Tc is the columns of blocks.
, r c row col T T g g where •    stands for ceil operation, g × g is the tessellation parameter.

Block Partitioning
For subimage Zt, its image domain Ft can be divided into It S disjoint nonempty regions according to the MHR, such that, F = {Fti; i = 1, 2, …, It S }, the corresponding sub-image is Zti = {za; a = 1, 2,…, Nti, (ua, va) ∈ Fti}, Nti is the number of pixels in the i-th region of block t, a is the local pixel index of the ith region of block t, S is the threshold of MHR, and each homogeneous region Fti corresponds to a certain land use and land cover (LULC) class or part of it. In order to better express the explicit and implicit features of blocks, the MST model is used to express blocks as image minimum spanning tree field (IMSTF). Let the pixels of block be mapped to vertices, the 8-neighborhood of pixels is taken as the adjacency relationship, and the nonsimilarity measure between the spectral measure vectors between the pixels is taken as the edge weight, the block can be expressed as graph Ɠ t = (Vt, Et, Ѱt), where Vt = {va = (ua, va, za); a = 1, 2, …, Nt} is the vertex set, Et = {ec = (vc1vc2); c = 1, 2, …, Ct, vc1, vc2 ∈ Vt} is the edge set, and Ѱt = {φc = φ(vc1, vc2); c = 1, 2, …, Ct} is the edge weights set. In order to improve the computational efficiency, each edge in Et has to be different. Therefore, the construction method for the edge set was proposed.
For a vertex a ∈ {row t , 2row t , …, (col t -1)×row t } or a ∈ {(col t -1)×row t +1, (col t -1)×row t +2, …, Nt}, which is in the last row or column but not the lower right corner, its adjacent edge set is Ea = {(vava +row t )} or Ea  , v a , z a ); a = 1, 2, . . . , N t } is the vertex set, E t = {e c = (v c1 v c2 ); c = 1, 2, . . . , C t , v c1 , v c2 ∈ V t } is the edge set, and Ψ t = {ϕ c = ϕ(v c1 , v c2 ); c = 1, 2, . . . , C t } is the edge weights set. In order to improve the computational efficiency, each edge in E t has to be different. Therefore, the construction method for the edge set was proposed.
As shown in Figure 1, for a vertex a ∈ ({1, 2, . . . , N t − row t }/{row t , 2row t , . . . , (col t −1)×row t }), which is not in the last row or column, its adjacent edge set is E a = {(v a v a +1 ), (v a v a +rowt ), (v a v a +rowt+1 ), (v a +rowt v a +1 )}. For a vertex a ∈ {row t , 2row t , . . . , (col t −1)×row t } or a ∈ {(col t −1)×row t +1, (col t −1)×row t +2, . . . , N t }, which is in the last row or column but not the lower right corner, its adjacent edge set is E a = {(v a v a +rowt )} or E a = {(v a v a +1 )}, respectively. The adjacent edge set of a vertex in the lower right corner Remote Sens. 2020, 12, 783 4 of 29 is E a = ∅. The edge set E = {E a ; a = 1, 2, . . . , N t } = {e c = (v c1 v c2 ); c = 1, 2, . . . , C t , v c1 ∈ V t , v c2 ∈ N p (v c1 )}, where N p (v c1 ) is the 8-neighborhood of vertex v c1 , c is the index of edge, C t is the number of edges, and its calculation equation as follows: ϕ c = ϕ(v c1 , v c2 ) is defined as the nonsimilarity measure of two vertices spectral vectors connected by e c : where d c1c2 and θ c1c2 are the euclidean distance and angle of the spectral vector between vertex v c1 and v c2 , respectively. These parameters can be expressed as: where t 1 , t 2 ∈ (0,1) are the coefficients of d c1c2 and θ c1c2 in Equation (4), respectively. After a lot of experiments, we found that the fine results could be obtained by setting t 1 = 0.1 and t 2 = 0.8.
Remote Sens. 2020, 12, 783 4 of 30 which is in the last row or column but not the lower right corner, its adjacent edge set is E a = {(v a v a +row t )} or E a = {(v a v a +1 )}, respectively. The adjacent edge set of a vertex in the lower right corner is E a = . The edge set E = {E a ; a =1, 2, …, N t } = {e c = (v c1 v c2 ); c = 1, 2, …, For the block graph Ɠ t , the MST rule was introduced to simplify Ɠ t , and T t = (V t¯, E t¯, Ψ t¯) was obtained. The vertex set V t¯= V t , edge set E t¯= {e g = (v g1 v g2 ); g = 1, 2, …, N t -1, v g1 , v g2  V t¯} , and Ψ t¯= {φ g = φ(v g1 , v g2 ); g = 1, 2, …, N t -1} is the set of connection weights corresponding to the connected edges. φ g is defined by Equation (4). T t has the following advantages: (1). The operation time complexity based on T t is O(|V t¯| ), which has linear complexity. (2). Every edge of E t¯i s a cut edge, such that deleting or restoring an edge will increase or decrease a region, which greatly simplifies the complexity of the partition operation, and each region corresponds to a tree. (3). The boundary of the target on T t is generated along the edge, which is helpful to express the complex and irregular boundary of the target. The classical algorithms for solving MST are the Prim [36] and Kruskal [37] algorithms, which are both based on the greedy strategy. The Prim algorithm is a vertex-driven mode with time Markov random field (HMRF) from modeling the neighborhood relationship of regions, which more ion and optimal image segmentation.
be divided into It S disjoint nonempty regions according the corresponding sub-image is Zti = {za; a = 1, 2,…, Nti, i-th region of block t, a is the local pixel index of the i-R, and each homogeneous region Fti corresponds to a or part of it. In order to better express the explicit and sed to express blocks as image minimum spanning tree ed to vertices, the 8-neighborhood of pixels is taken as larity measure between the spectral measure vectors ht, the block can be expressed as graph Ɠ t = (Vt, Et, s the vertex set, Et = {ec = (vc1vc2); c = 1, 2, …, Ct, vc1, vc2 ∈ 1, 2, …, Ct} is the edge weights set. In order to improve has to be different. Therefore, the construction method he proposed method extended the hidden Markov random field (HMRF) from modeling the eighborhood relationship of pixels to the neighborhood relationship of regions, which more ffectively obtained model parameter estimation and optimal image segmentation.

Given a high-resolution image
here Tr is the rows of blocks, and Tc is the columns of blocks.
, va) ∈ Fti}, Nti is the number of pixels in the i-th region of block t, a is the local pixel index of the ih region of block t, S is the threshold of MHR, and each homogeneous region Fti corresponds to a ertain land use and land cover (LULC) class or part of it. In order to better express the explicit and mplicit features of blocks, the MST model is used to express blocks as image minimum spanning tree ield (IMSTF). Let the pixels of block be mapped to vertices, the 8-neighborhood of pixels is taken as he adjacency relationship, and the nonsimilarity measure between the spectral measure vectors etween the pixels is taken as the edge weight, the block can be expressed as graph Ɠ t = (Vt, Et, t), where Vt = {va = (ua, va, za); a = 1, 2, …, Nt} is the vertex set, Et = {ec = (vc1vc2); c = 1, 2, …, Ct, vc1, vc2 ∈ t} is the edge set, and Ѱt = {φc = φ(vc1, vc2); c = 1, 2, …, Ct} is the edge weights set. In order to improve he computational efficiency, each edge in Et has to be different. Therefore, the construction method or the edge set was proposed.
As shown in Figure 1, for a vertex a ∈ ({1, 2, …, Ntrow (1). The operation time complexity based on T t is O(|V t¯| ), which has linear complexity. (2). Every edge of E t¯i s a cut edge, such that deleting or restoring an edge will increase or decrease a region, which greatly simplifies the complexity of the partition operation, and each region corresponds to a tree. (3). The boundary of the target on T t is generated along the edge, which is helpful to express the complex and irregular boundary of the target.
The classical algorithms for solving MST are the Prim [36] and Kruskal [37] algorithms, which are both based on the greedy strategy. The Prim algorithm is a vertex-driven mode with time complexity Remote Sens. 2020, 12, 783 5 of 29 of O(|E| + |V|log|V|) based on the Fibonacci heap and is suitable for dense graphs. Kruskal algorithm is an edge-driven algorithm with time complexity of O(|E|log|E|), and is appropriate for sparse graphs. Generally speaking, the graph is called sparse graph when 0 ≤ |E| ≤ |V|log|V|. Otherwise, the graph is dense. From Equation (3), n be divided into It S disjoint nonempty regions according }, the corresponding sub-image is Zti = {za; a = 1, 2,…, Nti, e i-th region of block t, a is the local pixel index of the i-HR, and each homogeneous region Fti corresponds to a s or part of it. In order to better express the explicit and used to express blocks as image minimum spanning tree ped to vertices, the 8-neighborhood of pixels is taken as ilarity measure between the spectral measure vectors ght, the block can be expressed as graph Ɠ t = (Vt, Et, } is the vertex set, Et = {ec = (vc1vc2); c = 1, 2, …, Ct, vc1, vc2 ∈ = 1, 2, …, Ct} is the edge weights set. In order to improve t has to be different. Therefore, the construction method is the 8-t is constructed by 8-neighborhood, which can be considered as a sparse graph. Thus, the solution for the block MST in this study was realized using the Kruskal algorithm.
At the start of the Kruskal algorithm, all vertices are treated as an MST. The merging tree operation by ascending order of edge weights is then defined until all MSTs are merged into one. Each merging tree operation ensures that the merged graph is a tree which does not contain any cycles. The algorithm's greedy strategy is embodied in the order of merging, such that the edge set is rearranged in ascending order of edge weights before iteration, so the edges with small edge weights are merged early in the iteration procedure. Traditional cycle detection algorithms used by the Kruskal are inefficient because of their search strategy, e.g., disjoint set algorithm and rank-and-path compression algorithm take time complexity O(log|V|) [38,39]. For improving the performance, a new cycle detection method based on tree marker vectors was proposed, and its pseudo-code o is presented in Algorithm 1.

Algorithm 1. Proposed image MST generation method based on Kruskal algorithm
Input: graph t is the total number of block t, and rowt and espectively. The number of blocks can be (1) locks.
into It S disjoint nonempty regions according onding sub-image is Zti = {za; a = 1, 2,…, Nti, of block t, a is the local pixel index of the ihomogeneous region Fti corresponds to a t. In order to better express the explicit and ess blocks as image minimum spanning tree es, the 8-neighborhood of pixels is taken as ure between the spectral measure vectors can be expressed as graph Ɠ t = (Vt, Et, set, Et = {ec = (vc1vc2); c = 1, 2, …, Ct, vc1, vc2 ∈ is the edge weights set. In order to improve fferent. Therefore, the construction method The edge set E was rearranged in ascending order of edge weight in step 1 of Algorithm 1. The subsequent cycle judgment process (steps 4, 5) had no order operation in the edge set. Thus, the E t¯o f the MST was strictly in the ascending order, which is an important characteristic of the Kruskal algorithm. The function unique (V t¯) was used for removing duplicate vertices. According step 3 in Algorithm 1, the cycle detection is very simple and takes constant time complexity, which shows the efficiency of proposed method.
The iteration procedure of the proposed image MST generation method based on the Kruskal algorithm is shown in Figure 2. The red and the black points represent the vertices inside and outside of the V t¯, respectively. The red solid lines and the black dotted lines represent the edges inside and outside of the E t¯, respectively. The red numbers were used to show the dynamic change procedure of the tree marker vector. The tree marker vector t was used to mark which tree each vertex belongs to. In the initial stage, each vertex was considered as a tree, as shown in Figure 2a. During the iteration procedure, the merging process for trees was accomplished by changing the tree marker of the vertices. As shown in Figure 2b, both vertices 1 and 2 belong to the tree marked 1 after changing the vertex marker from 2 to 1. The judgment of cycle becomes very simple after the introduction of the tree marker vector t, which only requires the determination of whether the tree markers of the two vertices in the merging edge are the same at each iteration. If they are the same, a cycle will be introduced after the merging operation, and the operation is refused. Oherwise, the operation is accepted. IMSTF reflects the spectral similarity and spatial clustering among pixels and expresses the explicit and implicit features of images. Essentially, it is an extension of the raster structure of image domain and explicitly depicts the connective relationship between vertices, especially its spatial clustering attribute, which has a positive effect on subsequent image processing. On this basis, in order to solve the problems of geometric noise and heterogeneity increase among homogeneous regions in high-resolution images, the similarity of spectral measures and the shape proportionality between regions were taken as the criterion for merging regions based on MHR.
Each vertex was assumed as a homogeneous region. Each edge e g E t¯w as a cut edge and was therefore regarded as the merging edge of its two connected regions. The criterion for merging regions was defined, which used the expression h(F ti , F ti  )  S to determine whether the regions merged or not. In the expression, F ti and F ti  are two regions in the block t connected by e g . Finally, IMSTF reflects the spectral similarity and spatial clustering among pixels and expresses the explicit and implicit features of images. Essentially, it is an extension of the raster structure of image domain and explicitly depicts the connective relationship between vertices, especially its spatial clustering attribute, which has a positive effect on subsequent image processing. On this basis, in order to solve the problems of geometric noise and heterogeneity increase among homogeneous regions in high-resolution images, the similarity of spectral measures and the shape proportionality between regions were taken as the criterion for merging regions based on MHR.
Each vertex was assumed as a homogeneous region. Each edge e g ∈E t¯w as a cut edge and was therefore regarded as the merging edge of its two connected regions. The criterion for merging regions was defined, which used the expression h(F ti , F ti ) ≤ S to determine whether the regions merged or not.
In the expression, F ti and F ti are two regions in the block t connected by e g . Finally, each cut edge was traversed, which was then followed by the merging procedure. The dividing scale S was used to describe the regional heterogeneity.
The merging criterion h(F i , F i ) is a linear combination of spectral similarity and shape proportionality between regions, and is defined as, where: where α ∈ [0, 1] is the coefficient of spectral similarity, and h color is the spectral similarity. σ b,merge , σ b,ti , and σ b,ti are the standard deviations of the spectral measure of the merged region (F ti ∪F ti ), region ti, and ti , respectively. N merge , N ti , and N ti are the pixels numbers of the merged region (F ti ∪F ti ), the region ti, and ti , respectively. 1/σ b 2 is the weight of band b and σ b 2 is the variance of band b.
The shape proportionality h shape is a linear combination of the regional shape compactness h compt and the smoothness h smooth , and γ ∈ [0, 1] is the weight of the regional shape compactness. l merge , l ti , and l ti are the outline lengths of the merged region (F ti ∪F ti ), region ti, and ti , respectively. b merge , b ti, and b ti are the regional outer rectangular perimeters of the merged region (F ti ∪F ti ), region ti, and ti , respectively. The solving procedure for the regional outer rectangular perimeter is shown in Figure 3a, where a green square represents a pixel in the region, and the regional outer rectangular perimeter is the number of red border squares. Figure 3b,c show the solving procedure of the regional outline. First, an augmented matrix was generated according to the coordinate range of the pixels of the region. As shown in Figure 3b, a green square represents an augmented pixel, and the pixel in the region is marked as 1 (see yellow squares). The element values of the four adjacent pixels were then summed up for each yellow square and used as replacement values. In Figure 3c, the yellow squares, with a value of 4, represent the interior pixels of the region, while the blue squares, with a value less than 4, represent the pixels on the regional outline. The length of the regional outline was determined by the number of blue squares.
Remote Sens. 2020, 12, 783 7 of 30 each cut edge was traversed, which was then followed by the merging procedure. The dividing scale S was used to describe the regional heterogeneity. The merging criterion h(F i , F i  ) is a linear combination of spectral similarity and shape proportionality between regions, and is defined as, where: where   [0, 1] is the coefficient of spectral similarity, and h color is the spectral similarity. σ b,merge , σ b,ti , and σ b, ti are the standard deviations of the spectral measure of the merged region (F ti F ti ), region ti, and ti, respectively. N merge , N ti , and N ti  are the pixels numbers of the merged region (F ti F ti  ), the region ti, and ti  , respectively. 1/σ b 2 is the weight of band b and σ b 2 is the variance of band b. The shape proportionality h shape is a linear combination of the regional shape compactness h compt and the smoothness h smooth , and   [0, 1] is the weight of the regional shape compactness. l merge , l ti , and l ti are the outline lengths of the merged region (F ti F ti ), region ti, and ti, respectively. b merge , b ti, and b ti are the regional outer rectangular perimeters of the merged region (F ti  F ti  ), region ti, and ti  , respectively.
The solving procedure for the regional outer rectangular perimeter is shown in Figure 3a, where a green square represents a pixel in the region, and the regional outer rectangular perimeter is the number of red border squares. Figures 3b, c show the solving procedure of the regional outline. First, an augmented matrix was generated according to the coordinate range of the pixels of the region. As shown in Figure 3b, a green square represents an augmented pixel, and the pixel in the region is marked as 1 (see yellow squares). The element values of the four adjacent pixels were then summed up for each yellow square and used as replacement values. In Figure 3c, the yellow squares, with a value of 4, represent the interior pixels of the region, while the blue squares, with a value less than 4, represent the pixels on the regional outline. The length of the regional outline was determined by the number of blue squares.
(a) (c) (b) Figure 3. Procedure for the regional outer rectangular perimeter and regional outline. (a) shows the outer rectangular perimeter of a region, and (b, c) show the procedure of region outline. The merging procedure is presented in Figure 4. Figure 4a,b illustrate the regional outer rectangular perimeter of two regions (b ti and b ti ) and the regional outline (l ti andl ti ) of region F ti and F ti , respectively. Figure 4c shows the regional outer rectangular perimeter b merge and regional outline l merge of the merged region (F ti ∪F ti ). The shape complementarity between regions F ti and F ti was good, so the shape proportionality h shape of merged region (F ti ∪F ti ) was low, which motivates the merging procedure of two regions (F ti and F ti ).
Remote Sens. 2020, 12, 783 8 of 30 The merging procedure is presented in Figure 4. Figures 4a-b illustrate the regional outer rectangular perimeter of two regions (b ti and b ti ) and the regional outline (l ti and l ti ) of region F ti and F ti , respectively. Figure 4c shows the regional outer rectangular perimeter b merge and regional outline l merge of the merged region (F ti  F ti  ). The shape complementarity between regions F ti and F ti  was good, so the shape proportionality h shape of merged region (F ti  F ti  ) was low, which motivates the merging procedure of two regions (F ti and F ti ). Figure 4. Procedure for merging two regions. (a) shows the regional outer rectangular perimeter bti and regional outline l ti of region F ti ; (b) shows the regional outer rectangular perimeter b ti  and regional outline l ti of region F ti ; (c) shows shows the regional outer rectangular perimeter b merge and regional outline lmerge of merged region (FtiFti).
The procedure for the MHR partitioning procedure based on IMSTF is shown in Figure 5. Figure 5a is a simulated-color image, which consists of three homogeneous regions and is marked by three colors: Red, green, and blue. The blue region contains a yellow pixel, which was used to simulate geometric noise. Figure 5b is the MST of the simulated image, which was able to distinctly delineate boundaries. Figure 5c represents the forest obtained by the MHR criterion, with each tree corresponding to a homogeneous region. The MHR describes the regional structure of the image. Based on the simulation results, the MHR based on IMSTF is able to identify the irregular boundaries of the target and deal with the geometric noise effectively. In addition, because E t¯i s arranged in ascending order of edge weights and appears as a disordered hash in space, the smaller edge weight means two regions connected by the edge are more similar in spectral, and the more preferred these two regions will be merged. The procedure of the MHR partitioning method based on the IMSTF can be considered as a procedure for global and progressive region growth. Thus, the problem of the number and spatial distribution of initial growing seed points can be avoided, which is difficult to evade in traditional region growth algorithms. The procedure for the global progressive region growth is similar to the process of simulated annealing. In the initial stage, the regions with the most similar spectral characteristics are preferentially merged to form several pure region growth cores. After more iterations, the (a) shows the regional outer rectangular perimeter b ti and regional outline l ti of region F ti ; (b) shows the regional outer rectangular perimeter b ti and regional outline l ti of region F ti ; (c) shows shows the regional outer rectangular perimeter b merge and regional outline l merge of merged region (F ti ∪F ti ).
The procedure for the MHR partitioning procedure based on IMSTF is shown in Figure 5. Figure 5a is a simulated-color image, which consists of three homogeneous regions and is marked by three colors: Red, green, and blue. The blue region contains a yellow pixel, which was used to simulate geometric noise. Figure 5b is the MST of the simulated image, which was able to distinctly delineate boundaries. Figure 5c represents the forest obtained by the MHR criterion, with each tree corresponding to a homogeneous region. The MHR describes the regional structure of the image. Based on the simulation results, the MHR based on IMSTF is able to identify the irregular boundaries of the target and deal with the geometric noise effectively.
Remote Sens. 2020, 12, 783 8 of 30 The merging procedure is presented in Figure 4. Figures 4a-b illustrate the regional outer rectangular perimeter of two regions (bti and bti′) and the regional outline (lti and lti′) of region Fti and Fti′, respectively. Figure 4c shows the regional outer rectangular perimeter bmerge and regional outline lmerge of the merged region (Fti∪Fti′). The shape complementarity between regions Fti and Fti′ was good, so the shape proportionality hshape of merged region (Fti∪Fti′) was low, which motivates the merging procedure of two regions (Fti and Fti′). (a) shows the regional outer rectangular perimeter bti and regional outline lti of region Fti; (b) shows the regional outer rectangular perimeter bti′ and regional outline lti′ of region Fti′; (c) shows shows the regional outer rectangular perimeter bmerge and regional outline lmerge of merged region (Fti∪Fti′).
The procedure for the MHR partitioning procedure based on IMSTF is shown in Figure 5. Figure  5a is a simulated-color image, which consists of three homogeneous regions and is marked by three colors: Red, green, and blue. The blue region contains a yellow pixel, which was used to simulate geometric noise. Figure 5b is the MST of the simulated image, which was able to distinctly delineate boundaries. Figure 5c represents the forest obtained by the MHR criterion, with each tree corresponding to a homogeneous region. The MHR describes the regional structure of the image. Based on the simulation results, the MHR based on IMSTF is able to identify the irregular boundaries of the target and deal with the geometric noise effectively. In addition, because Et¯ is arranged in ascending order of edge weights and appears as a disordered hash in space, the smaller edge weight means two regions connected by the edge are more similar in spectral, and the more preferred these two regions will be merged. The procedure of the MHR partitioning method based on the IMSTF can be considered as a procedure for global and progressive region growth. Thus, the problem of the number and spatial distribution of initial growing seed points can be avoided, which is difficult to evade in traditional region growth algorithms. The procedure for the global progressive region growth is similar to the process of simulated annealing. In the initial stage, the regions with the most similar spectral characteristics are preferentially merged to form several pure region growth cores. After more iterations, the complex In addition, because E t¯i s arranged in ascending order of edge weights and appears as a disordered hash in space, the smaller edge weight means two regions connected by the edge are more similar in spectral, and the more preferred these two regions will be merged. The procedure of the MHR partitioning method based on the IMSTF can be considered as a procedure for global and progressive region growth. Thus, the problem of the number and spatial distribution of initial growing seed points can be avoided, which is difficult to evade in traditional region growth algorithms. The procedure for the global progressive region growth is similar to the process of simulated annealing. In the initial stage, the regions with the most similar spectral characteristics are preferentially merged to form several pure region growth cores. After more iterations, the complex merging situation (with larger boundary weights) can be judged gradually, and the shape constraint plays a crucial role in making the region grow more compactly.
In terms of block tessellation, assuming each block is independent, each block can be processed in parallel. The parallel task pool Y = {y t = y(Z t , F t ); t = 1, 2, . . . , T} is formed by T blocks, where y(Z t , F t ) is the task of the MHR partitioning method which includes block graph representation, block MST solution, and MHR partitioning, y t is the result of block partitioning, and t is the task index. The parallel computing system was made up of P parallel computing units and used to process Y parallel. Since the amount of data in each block was approximately the same, the runtime for each task can also be considered as roughly similar. Therefore, T tasks can be aggregated into P task groups in a balanced task distribution manner. Each parallel computing unit executes its allocated tasks in sequence. Assume that an image is tessellated into 256 blocks with five parallel execution units, which are represented by five color blocks, as shown in Figure 6. The tasks distributed by each parallel computing unit are shown in the grid of color blocks at the left side of Figure 6. Each color block corresponds to a thread in the execution procedure. Threads with the same color formed thread blocks, which were executed in sequence by the corresponding color computing units. After all the blocks were divided into regions, the global result was reduced according to the block index.
Remote Sens. 2020, 12, 783 9 of 30 complex merging situation (with larger boundary weights) can be judged gradually, and the shape constraint plays a crucial role in making the region grow more compactly.
In terms of block tessellation, assuming each block is independent, each block can be processed in parallel. The parallel task pool Y = {y t = y(Z t , F t ); t = 1, 2, …, T} is formed by T blocks, where y(Z t , F t ) is the task of the MHR partitioning method which includes block graph representation, block MST solution, and MHR partitioning, y t is the result of block partitioning, and t is the task index. The parallel computing system was made up of P parallel computing units and used to process Y parallel. Since the amount of data in each block was approximately the same, the runtime for each task can also be considered as roughly similar. Therefore, T tasks can be aggregated into P task groups in a balanced task distribution manner. Each parallel computing unit executes its allocated tasks in sequence. Assume that an image is tessellated into 256 blocks with five parallel execution units, which are represented by five color blocks, as shown in Figure 6. The tasks distributed by each parallel computing unit are shown in the grid of color blocks at the left side of Figure 6. Each color block corresponds to a thread in the execution procedure. Threads with the same color formed thread blocks, which were executed in sequence by the corresponding color computing units. After all the blocks were divided into regions, the global result was reduced according to the block index.

Block Merging
In order to decrease the damage in landscape caused by regular tessellation, the results from block paritioning had to be merged. For the resulting partition y t of block t, its four neighborhood blocks were y t-1 , y t+1 , y t+Tr , and y t-Tr (see Figure 7a), where t = 1, 2, …, T, and T r was the row number of the block. The purpose of block merging is to stitch the bordered regions by reconstructing the correlation between adjacent blocks. In merging the block y t , the neighborhood blocks y t-1 , y t+1 , y t+Tr , and y t-Tr , were merged pair-wisely to avoid read-write conflicts, resulting in pairs (y t ,y t-1) , (y t , y t+Tr ), (y t , y t+1 ), and (y t , y t-Tr ). Figure 7a shows the association between blocks, where the red, blue, yellow, and green bidirectional arrows indicate the associations between y t and y t-1 , y t+Tr , y t+1 , y t-Tr , respectively. All bidirectional arrows identified by each color form the interblock merge set, which are represented by  t_t-1 ,  t_t+Tr ,  t_t+1 and  t_t-Tr , respectively. Each bidirectional arrow became the merge edge between two adjacent blocks, such that the red and yellow arrows were vertically correlated, while the blue and green arrows were horizontally correlated. In addition, the squares with different color bordered in Figure 7a represent different homogeneous regions, and the squares with different color fillings represent the divided regions of each block. These homogeneous regions were distributed in different blocks due to the influence of image tessellation. In effect, there was a redundancy problem in the merging procedure. For example, for the block y t , its left block y t-Tr and upper block y t-1 were merged, with the center on block y t-Tr-1 . To avoid this, the actual merging procedure only needed to consider the merging of the central block with its right and lower blocks,

Block Merging
In order to decrease the damage in landscape caused by regular tessellation, the results from block paritioning had to be merged. For the resulting partition y t of block t, its four neighborhood blocks were y t−1 , y t+1 , y t+Tr , and y t−Tr (see Figure 7a), where t = 1, 2, . . . , T, and T r was the row number of the block. The purpose of block merging is to stitch the bordered regions by reconstructing the correlation between adjacent blocks. In merging the block y t , the neighborhood blocks y t−1 , y t+1 , y t+Tr , and y t−Tr , were merged pair-wisely to avoid read-write conflicts, resulting in pairs (y t , y t−1) , (y t , y t+Tr ), (y t , y t+1 ), and (y t , y t−Tr ). Figure 7a shows the association between blocks, where the red, blue, yellow, and green bidirectional arrows indicate the associations between y t and y t−1 , y t+Tr , y t+1 ,y t−Tr , respectively. All bidirectional arrows identified by each color form the interblock merge set, which are represented by Λ t_t−1 , Λ t_t+Tr , Λ t_t+1 and Λ t_t−Tr , respectively. Each bidirectional arrow became the merge edge between two adjacent blocks, such that the red and yellow arrows were vertically correlated, while the blue and green arrows were horizontally correlated. In addition, the squares with different color bordered in Figure 7a represent different homogeneous regions, and the squares with different color fillings represent the divided regions of each block. These homogeneous regions were distributed in different blocks due to the influence of image tessellation. In effect, there was a redundancy problem in the merging procedure.
For example, for the block y t , its left block y t−Tr and upper block y t−1 were merged, with the center on block y t−Tr-1 . To avoid this, the actual merging procedure only needed to consider the merging of the central block with its right and lower blocks, as shown in Figure 7b. For the last row or the last column of blocks, only the right and lower blocks had to be merged with the central block (i.e., the lower right block did not need to be merged).
Remote Sens. 2020, 12,783 10 of 30 as shown in Figure 7b. For the last row or the last column of blocks, only the right and lower blocks had to be merged with the central block (i.e., the lower right block did not need to be merged).  Figure 8 shows the dependence relationship of horizontal interblock merging. For the three blocks y t-Tr , y t and y t+Tr , as shown in Figure 8a, the merged result of (y t-Tr , y t ) was considered as input for the merging procedure of (y t , y t+Tr ), which means that the merging result of (y t , y t+Tr ) shown in Figure 8c depends on the merged result of (y t-Tr , y t ) shown in Figure 8b. This dependence relationship created difficulties for the parallel merging procedure, which made it necessary to design a better parallel scheduling scheme to decouple the dependence relationship. The dependence relationship of the vertical merging process was similar to the horizontal procedure and is not discussed here. The merging procedure of blocks can be summarized as follows: For each merged edge in  t_t+1 , if two pixels connected by the edge belong to the same region, continue on with the next merged edge. Otherwise, determine whether or not to merge these two regions based on the equation h(F ti , F t+1 i )  S. If the criterion is satisfied, the two regions are merged. Otherwise, move on to the next merging edge. The flowchart of the algorithm is shown in Figure 9. Using the same principle, the  Figure 8 shows the dependence relationship of horizontal interblock merging. For the three blocks y t−Tr , y t and y t+Tr , as shown in Figure 8a, the merged result of (y t−Tr , y t ) was considered as input for the merging procedure of (y t , y t+Tr ), which means that the merging result of (y t , y t+Tr ) shown in Figure 8c depends on the merged result of (y t−Tr , y t ) shown in Figure 8b. This dependence relationship created difficulties for the parallel merging procedure, which made it necessary to design a better parallel scheduling scheme to decouple the dependence relationship. The dependence relationship of the vertical merging process was similar to the horizontal procedure and is not discussed here.
Remote Sens. 2020, 12,783 10 of 30 as shown in Figure 7b. For the last row or the last column of blocks, only the right and lower blocks had to be merged with the central block (i.e., the lower right block did not need to be merged).  Figure 8 shows the dependence relationship of horizontal interblock merging. For the three blocks y t-Tr , y t and y t+Tr , as shown in Figure 8a, the merged result of (y t-Tr , y t ) was considered as input for the merging procedure of (y t , y t+Tr ), which means that the merging result of (y t , y t+Tr ) shown in Figure 8c depends on the merged result of (y t-Tr , y t ) shown in Figure 8b. This dependence relationship created difficulties for the parallel merging procedure, which made it necessary to design a better parallel scheduling scheme to decouple the dependence relationship. The dependence relationship of the vertical merging process was similar to the horizontal procedure and is not discussed here. The merging procedure of blocks can be summarized as follows: For each merged edge in  t_t+1 , if two pixels connected by the edge belong to the same region, continue on with the next merged edge. Otherwise, determine whether or not to merge these two regions based on the equation h(F ti , F t+1 i )  S. If the criterion is satisfied, the two regions are merged. Otherwise, move on to the next merging edge. The flowchart of the algorithm is shown in Figure 9. Using the same principle, the merging set  t_t+Tr was assessed, then the process was continued until all blocks were merged. The merging procedure of blocks can be summarized as follows: For each merged edge in Λ t_t+1 , if two pixels connected by the edge belong to the same region, continue on with the next merged edge.
Otherwise, determine whether or not to merge these two regions based on the equation h(F ti , F t+1 i ) ≤ S. If the criterion is satisfied, the two regions are merged. Otherwise, move on to the next merging edge. The flowchart of the algorithm is shown in Figure 9. Using the same principle, the merging set Λ t_t+Tr was assessed, then the process was continued until all blocks were merged.
Since the merging procedure was carried out serially, when the number of blocks is large, the required processing time will be enormous. To further improve the speed of merging, a hierarchical parallel scheduling method was proposed. The computational dependency diagram for merging blocks is shown in Figure 10a. The computational dependency diagram can be divided into two parts according to the horizontal and vertical merging modes, as shown in Figures 10b, c, respectively. The red bidirectional arrow in Figure 10b indicates a horizontal first layer merging procedure, and the blue bidirectional arrow denotes a horizontal second layer merging procedure. Similarly, the green bidirectional arrow indicates a vertical first layer merging procedure, and the purple bidirectional arrow denotes a vertical second layer merging procedure.  Since the merging procedure was carried out serially, when the number of blocks is large, the required processing time will be enormous. To further improve the speed of merging, a hierarchical parallel scheduling method was proposed. The computational dependency diagram for merging blocks is shown in Figure 10a. The computational dependency diagram can be divided into two parts according to the horizontal and vertical merging modes, as shown in Figure 10b,c, respectively. The red bidirectional arrow in Figure 10b indicates a horizontal first layer merging procedure, and the blue bidirectional arrow denotes a horizontal second layer merging procedure. Similarly, the green bidirectional arrow indicates a vertical first layer merging procedure, and the purple bidirectional arrow denotes a vertical second layer merging procedure.
Since the merging procedure was carried out serially, when the number of blocks is large, the required processing time will be enormous. To further improve the speed of merging, a hierarchical parallel scheduling method was proposed. The computational dependency diagram for merging blocks is shown in Figure 10a. The computational dependency diagram can be divided into two parts according to the horizontal and vertical merging modes, as shown in Figures 10b, c, respectively. The red bidirectional arrow in Figure 10b indicates a horizontal first layer merging procedure, and the blue bidirectional arrow denotes a horizontal second layer merging procedure. Similarly, the green bidirectional arrow indicates a vertical first layer merging procedure, and the purple bidirectional arrow denotes a vertical second layer merging procedure.  The horizontal and vertical merging tasks of Figure 10b,c can be aggregated into the form presented in Figure 11, where Figure 11a,b are the first and second horizontal merging task sets and Figure 11c,d are the first and second vertically merging task sets, respectively. The index sets MI 1 , MI 2 , MI 3 , and MI 4 are the block pairs of the first horizontal layer, second horizontal layer, first vertical layer, and second vertical layer, respectively. These index sets can be obtained by the following equations: where T r and T c are the row and column numbers of the block, respectively, u = mod(a, b) is the remainder function, a is the dividend, b is the divisor, and u is the remainder. a, b, and u are all integers.   (15) where T r and T c are the row and column numbers of the block, respectively, u = mod(a, b) is the remainder function, a is the dividend, b is the divisor, and u is the remainder. a, b, and u are all integers.
In terms of the aggregated result of the horizontal and vertical merging tasks, the optimized procedure can be realized using a multithread parallel execution mode. Figure 12   In terms of the aggregated result of the horizontal and vertical merging tasks, the optimized procedure can be realized using a multithread parallel execution mode. Figure 12

Regional Hidden Markov Random Field -Fuzzy C-Means (RHMRF-FCM) Segmentation
The parallel MHR partitioning procedure divides the image domain into regions. However, these regions are optimal locally and not necessarily optimal globally. Thus, the general segmentation algorithms need to adopt an additional clustering process to obtain global optimal segmentation results. In this section, the RHMRF model was developed based on the idea of pixel HMRF-FCM image segmentation and using the region as the basic decisionmaking unit.

RHMRF-FCM Segmentation
The image domain F was divided into I S homogeneous regions using a parallel MHR partitioning method, such that F = {F i ; i = 1, 2, …, I S } is the set of pixels corresponding to each region Z i = {z a ; a = 1, 2,…, N i }, and its label set is L = {L i ; i = 1, 2, …, I S }. Then, the objective function of regionalized HMRF-FCM [40] can be written as: (16) where the first term on the right side (Equation (16)) is the linear weighted sum of information and its fuzzy relations with all the regions, which indicates the cost of label configuration and is the main factor for clustering decisionmaking. The second term refers to the regularization of KL divergence information. It introduces the relationship between the fuzzy relation and the local label of the region, which is the expression of the region's local constraints and can restrain the over-segmentation. λ  [0,1] is the fuzzy degree, which is used to coordinate the relationship between two decision items. R = {r ik ; i = 1, 2, …, I S , k = 1,2, …, K} is a fuzzy membership set, where r ik describes the degree which region i belongs to the jth cluster, and η = { ik ; i = 1, 2, …, I S , k = 1,2, …, K} is a prior distribution of HMRF model.  ik is the pointwise prior probability of the RHMRF model states [41], which can be expressed as: , , where NF i is the neighborhood relationship of the region F i , i is the neighborhood region index of F i , k is the labeled index, and β ∈ [0, 1] is a parameter of the clique potentials known as the inverse supercritical temperature.

Regional Hidden Markov Random Field-Fuzzy C-Means (RHMRF-FCM) Segmentation
The parallel MHR partitioning procedure divides the image domain into regions. However, these regions are optimal locally and not necessarily optimal globally. Thus, the general segmentation algorithms need to adopt an additional clustering process to obtain global optimal segmentation results. In this section, the RHMRF model was developed based on the idea of pixel HMRF-FCM image segmentation and using the region as the basic decisionmaking unit.

RHMRF-FCM Segmentation
The image domain F was divided into I S homogeneous regions using a parallel MHR partitioning method, such that F = {F i ; i = 1, 2, . . . , I S } is the set of pixels corresponding to each region Z i = {z a ; a = 1, 2, . . . , N i }, and its label set is L = {L i ; i = 1, 2, . . . , I S }. Then, the objective function of regionalized HMRF-FCM [40] can be written as: where the first term on the right side (Equation (16)) is the linear weighted sum of information and its fuzzy relations with all the regions, which indicates the cost of label configuration and is the main factor for clustering decisionmaking. The second term refers to the regularization of KL divergence information. It introduces the relationship between the fuzzy relation and the local label of the region, which is the expression of the region's local constraints and can restrain the over-segmentation. λ ∈ [0,1] is the fuzzy degree, which is used to coordinate the relationship between two decision items. R = {r ik ; i = 1, 2, . . . , I S , k = 1,2, . . . , K} is a fuzzy membership set, where r ik describes the degree which region i belongs to the j-th cluster, and η = {η ik ; i = 1, 2, . . . , I S , k = 1,2, . . . , K} is a prior distribution of HMRF model. η ik is the pointwise prior probability of the RHMRF model states [41], which can be expressed as: where NF i is the neighborhood relationship of the region F i , i is the neighborhood region index of otherwise, t(L i , L i ) = 0. D ik is the information about the class k in the data of region i and can be defined as: where θ = {θ k = {µ k , k }; k = 1, 2, . . . , K} is the set of model parameters, µ k is the mean vector with B dimension, and k is the covariance matrix with B×B dimension. The global optimal segmentation result based on region partitioning is the minimum solution of the objective function (Equation (16)). The explicit class attributes of each region can be obtained by de-fuzzifying the fuzzy relation matrix R.
At present, defuzzification is usually solved by maximizing the fuzzy relation, such that: By substituting Equations (17) and (18) into Equation (16), the final form of the objective function was obtained: where µ k , k , and r ik are unknown parameters. µ k and k are explicit parameters, such that the partial derivatives of the objective function for these two parameters can be solved by the following equations: where The fuzzy membership r ik satisfies the constraint condition k r ik = 1, such that the sum of the fuzzy relations of a region belonging to all classes is equal to 1. Therefore, a Lagrange function was constructed to solve the problem based on the constraints, Let the partial derivative of f (r i1 , r i2 , . . . , r iK ) to r ik be 0, Equation (28) was obtained by taking Equation (25) divided by Equation (27), Equation (29) was obtained by substituting Equation (18) into Equation (28),

Parallel RHMRF-FCM Segmentation
When the number of pixels is huge, or the image information is complex, the number of regions generated will also be very large, which leads to protracted runtime of the serial RHMRF-FCM algorithm. In addition, for multicore parallel computing systems, the communication time between cores is extremely expensive. Effectively reducing the amount of data synchronization can significantly improve the performance of parallel algorithms [42]. Therefore, a parallel RHMRF-FCM algorithm was proposed based on a master-slave parallel mode. The algorithm is described in detail as follows.
For a parallel computing system consisting of P parallel computing units and a coordinated computing unit, in order to make full use of computing performance brought by parallel computing, the dataset Z was equally partitioned into P subsets according to the number of regions, such that Z = {Z l ; l = 1, 2, . . . , p}. Z l = {Z lq ; q = 1, 2, . . . , N l } is the l-th subset data, Z lq = {z lqa ; a = 1, 2, . . . , N q } is the data of region q in the subsetl, a is the pixel index of region q, and z lqa is the spectral measure of pixel a in the region q of subsetl. N l is the number of regions in subset l, which can be computed as: Similarly, other local variables in the calculation process were defined according to the same rules. By the master-slave parallel mode, the flow of the parallel RHMRF-FCM algorithm can be summarized as follows: (1) Initialization algorithm parameters include the iteration index d = 0, clustering number K, fuzzy degree λ, and iteration threshold T. Initiate the p parallel computing units and the coordinated computing unit. The p parallel computing units are responsible for the calculation of p subsets of corresponding data sets, while the coordination computing unit is responsible for the dataset partitioning and the allocation of the local parameters.
(2) Coordinate computing unit for data partition and data transmission include the partition data set Z l = {Z lq ; q = 1, 2, . . . , N l }. Randomly initialize and partition the region fuzzy membership set R (d) = {R l (d) ; l = 1, 2, . . . , p}. Compute and divide the neighborhood region index set U = {U l ; l = 1, 2, . . . , p}, then transmit the partitioned data sets Z l , R l (d , ) and U l to the corresponding threads.
(4) Reduce M l (d) and µ l (d) using the coordinate computing unit from p computation units, and calculate M k (d) and µ k (d) using Equations (35) and (36). Then, synchronize the parametersM k (d) and µ k (d) to p parallel computing units.
(5) Compute l (d) = { lk (d) ; k = 1, 2, . . . , K} using p parallel computing units by: (6) Reduce l (d) from p computation units, and calculate the k (d) by Equation (38). Then, synchronize the k (d) to p parallel computational units: (8) If the iterations threshold is reached, the partition result L (d) is reduced and the loop is withdrawn. Otherwise, d = d+1, and then return to step (3) in the iteration process.
In order to more clearly define the local parallel computation, the parallel RHMRF-FCM algorithm was reduced and synchronized. The flowchart of the algorithm is shown in Figure 13. The local variables in each parallel computing unit were Z l , U l , R l , L l , η l , M l , and θ l . The variables stored in the computation procedure of the coordinated computing unit were M k , µ k , and k . As shown in Figure 13, the parallel computing units were responsible for the calculation of parameters related to data and the main computing carriers of the algorithm. The coordination computing unit was responsible for the generation of local parameters, the calculation of global parameters, and the synchronization of global parameters to each parallel computing unit. In addition, the parallel RHMRF-FCM algorithm only needed to specify and synchronize the global model parameters, which mainly included M l

Results and Discussion
Our previous works (serial method), presented in Reference [10], discussed the influence of the parameters of serial MHR partitioning method to the final segmentation result, obtained available range of parameters, and compared the multiresolution and watershed segmentation method of the eCognition software, so we will not be repeatedly discussing them. In this study, we focused on the discussion of the cost of the parallel segmentation method and parallel performance analysis. The cost analysis was mainly concerned with the impact of the proposed parallel segmentation schemes on global segmentation results under different data tessellation scales and was evaluated qualitatively and quantitatively. The parallel performance analysis involved the assessment of runtime and speedup [43] of parallel MHR partitioning, parallel regionalization segmentation, and the proposed segmentation method with different cores.

Cost Analysis
In order to qualitatively and quantitatively evaluate the influence of block tessellation on the proposed parallel segmentation method, a WorldView-3 Pan-Sharpen multi-band image with 1024 × 1024 pixels and spatial resolution of 0.5 m was used as the experimental image (see Figure 14a). The image shows a factory in an urban area and contains a large number of peripheral objects, such as cars, landscape vegetation, and traffic signs. Figure 14b presents a standard segmented map,

Results and Discussion
Our previous works (serial method), presented in Reference [10], discussed the influence of the parameters of serial MHR partitioning method to the final segmentation result, obtained available range of parameters, and compared the multiresolution and watershed segmentation method of the eCognition software, so we will not be repeatedly discussing them. In this study, we focused on the discussion of the cost of the parallel segmentation method and parallel performance analysis. The cost analysis was mainly concerned with the impact of the proposed parallel segmentation schemes on global segmentation results under different data tessellation scales and was evaluated qualitatively and quantitatively. The parallel performance analysis involved the assessment of runtime and speedup [43] of parallel MHR partitioning, parallel regionalization segmentation, and the proposed segmentation method with different cores.

Cost Analysis
In order to qualitatively and quantitatively evaluate the influence of block tessellation on the proposed parallel segmentation method, a WorldView-3 Pan-Sharpen multi-band image with 1024 × 1024 pixels and spatial resolution of 0.5 m was used as the experimental image (see Figure 14a). The image shows a factory in an urban area and contains a large number of peripheral objects, such as cars, landscape vegetation, and traffic signs. Figure 14b presents a standard segmented map, which was interpreted artificially and contains five homogeneous regions marked in red (I), green (II), blue (III), yellow (IV), and pink (V). Figure 13c shows a superimposed image of the boundaries of the standard segmentation result and the original image. In order to quantitatively and qualitatively assess the influence of different block tessellation scales on the final segmentation results, four comparative experiments were implemented. For experiments 1-3, the image was tessellated into 64, 16, and 4 blocks with tessellation scales 128 × 128, 256 × 256, and 512 × 512, respectively. In experiment 4, the image was not tessellated. Parallel MHR partitioning technique for the first three experiments, while the fourth setup, which served as the control, was divided using the serial MHR partitioning method. According to Reference [10], the segmentation parameters used in the four experiments were as follows: Dividing scale S = 220, the weight of spectral measure similarity  = 0.9, the weight of region shape compactness γ = 0.1, the threshold of iteration number T = 40, fuzzy degree λ = 0.1, homogeneous regions K = 5, and the clique potentials of label field β = 0.1.
The results are presented in Figure 15.  In order to quantitatively and qualitatively assess the influence of different block tessellation scales on the final segmentation results, four comparative experiments were implemented. For experiments 1-3, the image was tessellated into 64, 16, and 4 blocks with tessellation scales 128 × 128, 256 × 256, and 512 × 512, respectively. In experiment 4, the image was not tessellated. Parallel MHR partitioning technique for the first three experiments, while the fourth setup, which served as the control, was divided using the serial MHR partitioning method. According to Reference [10], the segmentation parameters used in the four experiments were as follows: Dividing scale S = 220, the weight of spectral measure similarity α = 0.9, the weight of region shape compactness γ = 0.1, the threshold of iteration number T = 40, fuzzy degree λ = 0.1, homogeneous regions K = 5, and the clique potentials of label field β = 0.1.
The results are presented in Figure 15. Figure 15a-p show the results of experiments 1-4, respectively. Figure 15a,e,i,m show the results using MHR partitioning procedure, with different colors used to identify partitions. The resulting subdividing outlines overlaid over the original images are shown in Figure 15b,f,j,n. Figure 15c,g,k,o present the optimal segmentation results by parallel RHMRF-FCM, and Figure 15d,h,l,p are the optimal segmentation outlines overlaid on the original images.
As shown in Figure 15a,e,i,m, block merging can significantly reduce block tessellation. However, some regions were still affected by high heterogeneity in pixel spectral information, which results in difficulties to effectively process block merging. Areas (e.g., roads and parking lots) containing a substantial amount of peripheral objects (e.g., vehicles, traffic markers, and shadows) can exhibit high levels of spectral heterogeneity. When block tessellation boundaries encounter these types of objects, they may yield unpredictable effects on the merging results. For terrains with low spectral heterogeneity, such as the factory roof in the image, block merging can eliminate the effects of block tessellation. Comparing the MHR partitioning results in Figure 15a,e,i,m, the increase in block tessellation scale resulted in a significant reduction of the influence of block tessellation on segmentation. The result from the first experiment was visibly different from the other three. The effect of block tessellation could not be effectively eliminated in this setup, which was impacted by the very small tessellation scale (128 × 128). The results from the third experiment were found to be comparable with the control group. The superposition analysis shows that the outlines of Figure 15b,f,j,n matched the complex boundaries of objects well. Generally, the segmentation results were reasonably similar, with slight differences found in local areas. Also, excessive segmentation was found to decrease with the increase in the block tessellation scale.
Remote Sens. 2020, 12, 783 20 of 30 As shown in Figures 15a, e, i, m, block merging can significantly reduce block tessellation. However, some regions were still affected by high heterogeneity in pixel spectral information, which results in difficulties to effectively process block merging. Areas (e.g., roads and parking lots) containing a substantial amount of peripheral objects (e.g., vehicles, traffic markers, and shadows) can exhibit high levels of spectral heterogeneity. When block tessellation boundaries encounter these types of objects, they may yield unpredictable effects on the merging results. For terrains with low spectral heterogeneity, such as the factory roof in the image, block merging can eliminate the effects of block tessellation. Comparing the MHR partitioning results in Figures 15a, e, i, m, the increase in block tessellation scale resulted in a significant reduction of the influence of block tessellation on segmentation. The result from the first experiment was visibly different from the other three. The effect of block tessellation could not be effectively eliminated in this setup, which was impacted by the very small tessellation scale (128 × 128). The results from the third experiment were found to be comparable with the control group. The superposition analysis shows that the outlines of Figures 15b, f, j, n matched the complex boundaries of objects well. Generally, the segmentation results were reasonably similar, with slight differences found in local areas. Also, excessive segmentation was found to decrease with the increase in the block tessellation scale. Comparing the optimal segmentation results in Figure 15c,g,k,o, the four optimal segmentation results with different tessellation scales were generally similar but with slight variations locally. When the tessellation scale was increased, the difference between the local segmentation results and the control group showed a decreasing trend. These results show that the proposed RHMRF-FCM optimal segmentation method can accurately assess the class of homogeneous regions globally and can further reduce the influence of block tessellation. The optimal segmentation results had slight differences locally, mainly due to the influence of block tessellation. Based on the results of the superposition analysis (Figure 15d,h,l,p), the segmentation boundaries are able to delineate irregular boundaries more accurately, can overcome the influence of most peripheral objects, and are able to generate better segmentation results.
In order to quantitatively evaluate segmentation accuracies, the confusion matrix [44] was used in the quantitative segmentation accuracy evaluation model. The manual interpretation result of Figure 14b was used as the standard segmentation results, and the evaluation results are listed in Table 1. The results show that the segmentation accuracies of the first experiment were low, which suggests that the 128 × 128 block tessellation scale had a significant influence on the segmentation accuracy. The overall accuracies in the last three groups were at least 0.82, and the kappa values were at least 0.78. Also, the user's and producer's accuracies for each homogeneous region were mostly above 0.8. The results suggest the total segmentation accuracy was able to reach the standard for an excellent classifier. The overall accuracy and kappa value of the second group were lower than the control group by 0.03 and 0.04, respectively. The overall accuracy and kappa value of the third group were lower than the control group by 0.02 and 0.02, respectively. These suggest that the 256 × 256 and 512 × 512 tessellation scales had less effect on the final segmentation results. With the increasing block tessellation scale, the precision index generally increased in small amounts. In addition, the producer's accuracies in regions III and IV were only 0.64-0.67 and 0.44-0.48, respectively, which were different from the accuracies of other regions. For further comparison, the standard segmentation results and the control group segmentation results were compared, as shown in Figure 16. Figure 16a is the standard segmentation result. The result's outlines were extracted and overlaid on the original image, as shown in Figure 16c. Figure 16b shows the result of the control group segmentation. The result's outlines were extracted and overlaid on the original image, shown in Figure 16d. Homogeneous regions III and IV correspond to the blue and yellow areas in Figure 16a,b, respectively. Compared with Figure 16c, homogeneous area III corresponds to dark pavement and shadows and contains a large number of peripheral features, such as vehicles and traffic signs, while homogeneous area IV corresponds to independent trees, lawns, etc. Visually, the two homogeneous regions show strong heterogeneity in spectral characteristics, and the color difference between them is not distinctly apparent, which exacerbates the difficulty of distinguishing the two homogeneous regions. In order to more accurately reflect the distribution of spectral characteristics of the two homogeneous regions, their pixel indexes were extracted, and their spectral information was obtained from the original map. The spectral distributions of the two homogeneous regions are shown in Figure 17. The red points represent the spectral information of pixels in region III, while the blue points represent the spectral information for pixels in region IV. As seen from the graph, the entanglement of spectral features in the homogeneous regions III and IV was the fundamental reason for their low segmentation accuracy. homogeneous regions, their pixel indexes were extracted, and their spectral information was obtained from the original map. The spectral distributions of the two homogeneous regions are shown in Figure 17. The red points represent the spectral information of pixels in region III, while the blue points represent the spectral information for pixels in region IV. As seen from the graph, the entanglement of spectral features in the homogeneous regions III and IV was the fundamental reason for their low segmentation accuracy.

Performance Analysis
A Dell Power Edge R930 server with four sockets was used for the parallel performance analysis of the proposed parallel image segmentation method. The server had four Intel Xeon E7-8880 v3 @2.3GHz CPU, and each CPU had 18 cores (72 cores in total), 96 GB, 2133 MT/s memory, and a Samsung EVO Pro 512 GB SSD used as the operating system disk. The parallel computing toolbox of Matlab2018a software running on the Windows 10 operating system was used to implement the proposed method, which provided a parallel way through a multiworker mode.
In order to verify the parallel performance of the proposed parallel image segmentation method, two WorldView-3 Pan-Sharpened multiband images with 8192 × 8192 pixels/0.5-m spatial resolution/urban area of Saltlake, USA, and 7680 × 7680 pixels/0.3-m spatial resolution/Aeropueto de Madrid-Barajas Adolfo Suarez, Madrid, Spain, were used as experimental images (see Figure 18a pixel indexes were extracted, and their spectral information was obtained from the original map. The spectral distributions of the two homogeneous regions are shown in Figure 17. The red points represent the spectral information of pixels in region III, while the blue points represent the spectral information for pixels in region IV. As seen from the graph, the entanglement of spectral features in the homogeneous regions III and IV was the fundamental reason for their low segmentation accuracy.

Performance Analysis
A Dell Power Edge R930 server with four sockets was used for the parallel performance analysis of the proposed parallel image segmentation method. The server had four Intel Xeon E7-8880 v3 @2.3GHz CPU, and each CPU had 18 cores (72 cores in total), 96 GB, 2133 MT/s memory, and a Samsung EVO Pro 512 GB SSD used as the operating system disk. The parallel computing toolbox of Matlab2018a software running on the Windows 10 operating system was used to implement the proposed method, which provided a parallel way through a multiworker mode.
In order to verify the parallel performance of the proposed parallel image segmentation method, two WorldView-3 Pan-Sharpened multiband images with 8192 × 8192 pixels/0.5-m spatial resolution/urban area of Saltlake, USA, and 7680 × 7680 pixels/0.3-m spatial resolution/Aeropueto de Madrid-Barajas Adolfo Suarez, Madrid, Spain, were used as experimental images (see Figure 18a,

Performance Analysis
A Dell Power Edge R930 server with four sockets was used for the parallel performance analysis of the proposed parallel image segmentation method. The server had four Intel Xeon E7-8880 v3 @2.3GHz CPU, and each CPU had 18 cores (72 cores in total), 96 GB, 2133 MT/s memory, and a Samsung EVO Pro 512 GB SSD used as the operating system disk. The parallel computing toolbox of Matlab2018a software running on the Windows 10 operating system was used to implement the proposed method, which provided a parallel way through a multiworker mode.
In order to verify the parallel performance of the proposed parallel image segmentation method, two WorldView-3 Pan-Sharpened multiband images with 8192 × 8192 pixels/0.5-m spatial resolution/urban area of Saltlake, USA, and 7680 × 7680 pixels/0.3-m spatial resolution/Aeropueto de Madrid-Barajas Adolfo Suarez, Madrid, Spain, were used as experimental images (see Figure 18a,e), respectively. The segmentation parameters of the first experiments were set as same as those in Section 3.1, except the homogeneous regions K = 6 and the block tessellation scale 512 × 512. The segmentation parameters of second experiments were set as same as those in Section 3.1, except the homogeneous regions K = 6, the block tessellation scale was 512 × 512, the dividing scale S = 2000, the weight of spectral measure similarity α = 0.7, and the weight of region shape compactness γ = 0.4. In addition, considering the main land cover type was bare land, which had high heterogeneity in second experiment image (Figure 18e), the parameters (S, α, γ) of first experiment were different than the second experiment. S and γ were increased and α was decreased to improve the performance of the noice immunity of the parallel MHR method.  Figure 18. Segmentation results of the whole scene images.
In order to quantitatively evaluate the performance of the proposed parallel segmentation method, the runtime and speedup of the parallel MHR partitioning method, the parallel RHMRF-FCM segmentation method, and the proposed method were compared using different core numbers and segmentation parameters.
The runtime of the parallel MHR partitioning technique was tested using 8, 16, 24, 32, 40, 48, 56, and 64 cores. The results, as illustrated in Figure 19 by a blue curve, show the running time decreased rapidly with the increasing core number until it reached 48. At such point, the runtime became stable. This indicates that the best core number option for this experiment was 48. The speedup in the parallel MHR segmentation method is shown by the red curve in Figure 19. The results show that speedup increased rapidly as the core number increased until it reached 48, wherein the speedup tended to stabilize. The serial runtime of the minimum heterogeneous segmentation of an image with 8192 × 8192 pixels was then estimated by quadratic fitting. The runtime values for 128 × 128, 256 × 256, 512 × 512, 1024 × 1024, and 2048 × 2048 pixel images were used (see Figure 20a). The fitting function and fitting results are presented in Equation (40) and The experimental results are shown in Figure 18. The magnified areas highlighted in red and black boxes in Figure 18a,e are presented in Figure 18i,q (the images of Utah state office buildings and Salt Palace Convention Center) and Figure 18m,u (the images of airport terminal and parking lot), respectively. Figure 18b-d,f-h show the results of the segmentation outlines overlaid on the image, respectively. Figure 18j-l,n-p are the optimal segmentation results, respectively, while Figure 18r-t,v-x are the optimal segmentation outlines overlaid on the image, respectively. The findings suggest that visually, the results obtained by parallel MHR partitioning were relatively good, and the segmentation of complex object boundaries were complete and accurate.
In order to quantitatively evaluate the performance of the proposed parallel segmentation method, the runtime and speedup of the parallel MHR partitioning method, the parallel RHMRF-FCM segmentation method, and the proposed method were compared using different core numbers and segmentation parameters.
The runtime of the parallel MHR partitioning technique was tested using 8, 16, 24, 32, 40, 48, 56, and 64 cores. The results, as illustrated in Figure 19 by a blue curve, show the running time decreased rapidly with the increasing core number until it reached 48. At such point, the runtime became stable. This indicates that the best core number option for this experiment was 48. The speedup in the parallel MHR segmentation method is shown by the red curve in Figure 19. The results show that speedup increased rapidly as the core number increased until it reached 48, wherein the speedup tended to Remote Sens. 2020, 12, 783 25 of 29 stabilize. The serial runtime of the minimum heterogeneous segmentation of an image with 8192 × 8192 pixels was then estimated by quadratic fitting. The runtime values for 128 × 128, 256 × 256, 512 × 512, 1024 × 1024, and 2048 × 2048 pixel images were used (see Figure 20a). The fitting function and fitting results are presented in Equation (40) and Figure 20b, respectively. f(N) = 4.59 × 10 −10 N 2 + 0.02N − 1.04 × 10 4 (40) where N is the pixel number of the test image. The speedup values may have been overestimated in these experiments. On the one hand, using the theoretical analysis, the serial running time of the 8192 × 8192 pixel image was estimated using the fitting function, which may have resulted in significant deviations from the actual runtime. On the other hand, the time complexity of the serial algorithm in the proposed method was high given that the algorithm uses the MST model, and the calculation efficiency when implementing the model may not have been optimal. However, the speedup curve can also roughly reflect the proposed parallel MHR partitioning method, which had a significantly parallel performance.
Remote Sens. 2020, 12, 783 26 of 30 where N is the pixel number of the test image. The speedup values may have been overestimated in these experiments. On the one hand, using the theoretical analysis, the serial running time of the 8192 × 8192 pixel image was estimated using the fitting function, which may have resulted in significant deviations from the actual runtime. On the other hand, the time complexity of the serial algorithm in the proposed method was high given that the algorithm uses the MST model, and the calculation efficiency when implementing the model may not have been optimal. However, the speedup curve can also roughly reflect the proposed parallel MHR partitioning method, which had a significantly parallel performance. The runtime of the parallel RHMRF-FCM segmentation method was tested using 1,2,4,8,16,24,32,40,48,56, and 64 cores. The results, presented in a blue curve in Figure 20, show that runtime decreased rapidly with the increasing core number until the core number reached 16, and then the runtime tended to stabilize. When the core number was at 56, the runtime was minimum, which was only 24% of the serial algorithm, and the acceleration effect was significant. Therefore, for the parallel computing environment, using 56 cores for block parallel partitioning was the optimal choice, while 16 was the choice with the high-cost performance. The speedup of the parallel RHMRF-FCM segmentation method is shown by the red curve in Figure 21. As the number of cores increased, the speedup increased rapidly, until the core number reached 16. At such point, the speedup was at 3.5 and became stable, and the effect of parallel acceleration was also significant. In terms of the magnitude of speedup, the parallel optimal segmentation method was two orders where N is the pixel number of the test image. The speedup values may have been overestimated in these experiments. On the one hand, using the theoretical analysis, the serial running time of the 8192 × 8192 pixel image was estimated using the fitting function, which may have resulted in significant deviations from the actual runtime. On the other hand, the time complexity of the serial algorithm in the proposed method was high given that the algorithm uses the MST model, and the calculation efficiency when implementing the model may not have been optimal. However, the speedup curve can also roughly reflect the proposed parallel MHR partitioning method, which had a significantly parallel performance. The runtime of the parallel RHMRF-FCM segmentation method was tested using 1,2,4,8,16,24,32,40,48,56, and 64 cores. The results, presented in a blue curve in Figure 20, show that runtime decreased rapidly with the increasing core number until the core number reached 16, and then the runtime tended to stabilize. When the core number was at 56, the runtime was minimum, which was only 24% of the serial algorithm, and the acceleration effect was significant. Therefore, for the parallel computing environment, using 56 cores for block parallel partitioning was the optimal choice, while 16 was the choice with the high-cost performance. The speedup of the parallel RHMRF-FCM segmentation method is shown by the red curve in Figure 21. As the number of cores increased, the speedup increased rapidly, until the core number reached 16. At such point, the speedup was at 3.5 and became stable, and the effect of parallel acceleration was also significant. In terms of the magnitude of speedup, the parallel optimal segmentation method was two orders The runtime of the parallel RHMRF-FCM segmentation method was tested using 1,2,4,8,16,24,32,40,48,56, and 64 cores. The results, presented in a blue curve in Figure 20, show that runtime decreased rapidly with the increasing core number until the core number reached 16, and then the runtime tended to stabilize. When the core number was at 56, the runtime was minimum, which was only 24% of the serial algorithm, and the acceleration effect was significant. Therefore, for the parallel computing environment, using 56 cores for block parallel partitioning was the optimal choice, while 16 was the choice with the high-cost performance. The speedup of the parallel RHMRF-FCM segmentation method is shown by the red curve in Figure 21. As the number of cores increased, the speedup increased rapidly, until the core number reached 16. At such point, the speedup was at 3.5 and became stable, and the effect of parallel acceleration was also significant. In terms of the magnitude of speedup, the parallel optimal segmentation method was two orders lower than the parallel MHR segmentation technique. This is because the serial RHMRF-FCM algorithm was affected by its iterative solution method, and it was difficult to achieve large-scale performance improvements through parallel computing.
Remote Sens. 2020, 12,783 27 of 30 Figure 22 shows the runtime (in a blue curve) and speedup (in a red curve) of the proposed parallel segmentation method. The overall trend of the curves shows distinct similarities with the results of the MHR segmentation ( Figure 19). With increasing the core number, the runtime and the speedup rapidly decreased and increased, respectively. After 48 cores, the curves tended to stabilize. The best parallel performance was achieved at 56 cores, with a runtime of 3685 s (about 1 h).

Discussion
In the parallel cost analysis experiment, we qualitatively and quantitatively evaluated the influence of different block tessellation scale on the final segmentation results. The results suggest the following conclusions: (a) With the increasing block tessellation scale, the similarity increased between the parallel segmentation results and serial segmentation results in terms of visual examination and quantitative accuracy evaluation results. This demonstrates the effectiveness of the proposed parallel MHR partitioning method. (b) For a WorldView-3 remote sensing image with 0.5-m resolution, the most reasonable block tessellation scale was 512 × 512. Using such setting, the accuracy loss in segmentation was constrained within 2% compared with the control group, which is generally acceptable in most applications. In the parallel performance analysis experiment, relatively good results were obtained using the parallel MHR partitioning method. The segmentation results of complex object boundaries were found to be complete and accurate. For a WorldView-3 Pan-Sharpened multiband remote sensing image with 8192 × 8192 pixels and 0.5-m spatial resolution, the running time and speedup curve for the parallel MHR partitioning, parallel RHMRF-FCM, and the proposed method were examined using four sockets server with 72 Cores. The results demonstrate that the proposed method has a high parallel performance.  Figure 22 shows the runtime (in a blue curve) and speedup (in a red curve) of the proposed parallel segmentation method. The overall trend of the curves shows distinct similarities with the results of the MHR segmentation ( Figure 19). With increasing the core number, the runtime and the speedup rapidly decreased and increased, respectively. After 48 cores, the curves tended to stabilize. The best parallel performance was achieved at 56 cores, with a runtime of 3685 s (about 1 h).

Conclusions
Remote Sens. 2020, 12,783 27 of 30 Figure 22 shows the runtime (in a blue curve) and speedup (in a red curve) of the proposed parallel segmentation method. The overall trend of the curves shows distinct similarities with the results of the MHR segmentation ( Figure 19). With increasing the core number, the runtime and the speedup rapidly decreased and increased, respectively. After 48 cores, the curves tended to stabilize. The best parallel performance was achieved at 56 cores, with a runtime of 3685 s (about 1 h).

Discussion
In the parallel cost analysis experiment, we qualitatively and quantitatively evaluated the influence of different block tessellation scale on the final segmentation results. The results suggest the following conclusions: (a) With the increasing block tessellation scale, the similarity increased between the parallel segmentation results and serial segmentation results in terms of visual examination and quantitative accuracy evaluation results. This demonstrates the effectiveness of the proposed parallel MHR partitioning method. (b) For a WorldView-3 remote sensing image with 0.5-m resolution, the most reasonable block tessellation scale was 512 × 512. Using such setting, the accuracy loss in segmentation was constrained within 2% compared with the control group, which is generally acceptable in most applications. In the parallel performance analysis experiment, relatively good results were obtained using the parallel MHR partitioning method. The segmentation results of complex object boundaries were found to be complete and accurate. For a WorldView-3 Pan-Sharpened multiband remote sensing image with 8192 × 8192 pixels and 0.5-m

Discussion
In the parallel cost analysis experiment, we qualitatively and quantitatively evaluated the influence of different block tessellation scale on the final segmentation results. The results suggest the following conclusions: (a) With the increasing block tessellation scale, the similarity increased between the parallel segmentation results and serial segmentation results in terms of visual examination and quantitative accuracy evaluation results. This demonstrates the effectiveness of the proposed parallel MHR partitioning method. (b) For a WorldView-3 remote sensing image with 0.5-m resolution, the most reasonable block tessellation scale was 512 × 512. Using such setting, the accuracy loss in segmentation was constrained within 2% compared with the control group, which is generally acceptable in most applications. In the parallel performance analysis experiment, relatively good results were obtained using the parallel MHR partitioning method. The segmentation results of complex object boundaries were found to be complete and accurate. For a WorldView-3 Pan-Sharpened multiband remote sensing image with 8192 × 8192 pixels and 0.5-m spatial resolution, the running time and speedup curve for the parallel MHR partitioning, parallel RHMRF-FCM, and the proposed method were examined using four sockets server with 72 Cores. The results demonstrate that the proposed method has a high parallel performance.

Conclusions
High-resolution imagery offers rich geometric, texture, and spectral information, and is characterized by big data. These features bring new challenges to traditional image segmentation methods. To this end, this study proposed an efficient parallel high-resolution image segmentation method based on the MST and RHMRF-FCM model. There are two main contributions of this research: (a) On the basis of image domain decomposition by regular tessellation, block partitioning was executed in parallel by the MHR technique and integrated with the block IMSTF model, which effectively expressed the image's explicit and implicit features. The parallel block merging method was then proposed to decrease the impact of landscape caused by image tessellation. (b) A parallel RHMRF-FCM segmentation method based on master-slave parallel mode was proposed to improve the computing performance of the serial RHMRF-FCM method, which effectively decreased the amount of data synchronization in parallel algorithms.
The proposed method has the following advantages: (a) By introducing the shape proportionality expressed on the IMSTF model, it serves as a much improved version of the MHR technique that is able to solve the problems of geometric noise and increased heterogeneity in homogeneous regions of high-resolution images. (b) The new method provides an efficient segmentation framework for big-scale high-resolution images using a parallel MHR partitioning technique, a parallel merging method based on image tessellation, and a parallel RHMRF-FCM segmentation method with low data synchronization traffic parallel scheme.
Further research can focus on the following directions: (a) Other parallel partitioning methods for big-scale high-resolution images can be explored through the creation of a more effective coordination mechanism between blocks. (b) Also, the proposed method can be applied to examine different types of remote sensing data, such as hyperspectral image and SAR image, and other applications, such as feature extraction and target identification.
Author Contributions: This research mainly performed and prepared by W.L. and Y.L. W.L. contributed with ideas, conceived, and designed the study. Y.L. supervised the study and his comments were considered throughout the paper. All authors have read and agreed to the published version of the manuscript.