An Efﬁcient Parallel Multi-Scale Segmentation Method for Remote Sensing Imagery

: Remote sensing (RS) image segmentation is an essential step in geographic object-based image analysis (GEOBIA) to ultimately derive “meaningful objects”. While many segmentation methods exist, most of them are not efﬁcient for large data sets. Thus, the goal of this research is to develop an efﬁcient parallel multi-scale segmentation method for RS imagery by combining graph theory and the fractal net evolution approach (FNEA). Speciﬁcally, a minimum spanning tree (MST) algorithm in graph theory is proposed to be combined with a minimum heterogeneity rule (MHR) algorithm that is used in FNEA. The MST algorithm is used for the initial segmentation while the MHR algorithm is used for object merging. An efﬁcient implementation of the segmentation strategy is presented using data partition and the “reverse searching-forward processing” chain based on message passing interface (MPI) parallel technology. Segmentation results of the proposed method using images from multiple sensors (airborne, SPECIM AISA EAGLE II, WorldView-2, RADARSAT-2) and different selected landscapes (residential/industrial, residential/agriculture) covering four test sites indicated its efﬁciency in accuracy and speed. We conclude that the proposed method is applicable and efﬁcient for the segmentation of a variety of RS imagery (airborne optical, satellite optical, SAR, high-spectral), while the accuracy is comparable with that of the FNEA method.


Introduction
Remote sensing (RS) image segmentation is an essential step to ultimately derive "meaningful objects" in geographic object-based image analysis (GEOBIA) [1,2].The optimal segmentation of an image should result in a balance of over-segmentation and under-segmentation of the scene.Many segmentation methods for RS imagery exist.They can be grouped into four categories: (a) pixel-based (thresholding, clustering, etc.) [3,4]; (b) edge detection (boundary tracking, curve fitting, etc.) [5]; (c) region-based (region growing, watershed, split and merge, level set, etc.) [6]; (d) others (wavelet, neural networks, mean-shift, fuzzy sets, etc.) [7,8].Region-based segmentation methods are the dominant methods that can be used to create compact regions and identify suitable scales [9].The fractal net evolution approach (FNEA) and graph method are the two dominant region-based methods.The FNEA, proposed by Baatz and Schäpe in 2000 [10] and commercialized by the software eCognition (http://www.ecognition.com), is a bottom-up region growing and merging technique based on local criteria.The neighboring image objects are merged to form bigger objects based on a minimum heterogeneity rule (MHR) algorithm.Thus, it builds a multi-scale hierarchy structure using a merging technique.In addition, spectral, shape and textural information can be extracted as well, which is helpful for advanced applications [10].However, this method-as with most local region growing techniques-has some limitations: it is time-consuming to create a large number of initial image objects starting with individual pixels [11]; and it is suboptimal to choose a sequence of starting points in a random sequence, due to the percolation patterns or clusters which increasingly appear, starting at a certain density of points [10].Graph methods explicitly organize the image elements into mathematically sound structure and make the formulation of the problem more flexible and the computation more efficient [12].There are four basic algorithms: best merge [13], minimum spanning tree (MST) [14], minimum mean cut [15,16], and normalized cut [17].The best merge and MST algorithms follow a bottom-up approach: they start from a segmented image, i.e., each pixel constitutes a separate segment; then, these algorithms iteratively merge two adjacent segments.The minimum mean cut and normalized cut algorithms follow a top-down approach: the whole image is one segment initially; then, in each step, one of the segments is split into parts [18].These graph-based algorithms have been applied to different RS imagery (such as hyperspectral, panchromatic, multispectral, LiDAR, SAR, etc.) [19][20][21][22].
Recently, two obvious trends within RS image segmentation can be observed: (a) new theories are constantly applied, which leads to more and more specific RS image segmentation methods and methodologies; (b) different segmentation algorithms have been combined to overcome the disadvantages of each individual segmentation algorithm [1].The key problem here is how to combine these algorithms to fully employ the advantages of various algorithms and complement their shortcomings.For instance, Li et al. proposed an efficient multi-scale segmentation method by combining a statistical region merging algorithm and a MHR algorithm, whereby the results showed that the combined algorithm is able to capture the main structural components of imagery by a simple but effective statistical analysis, and it is able to cope with significant noise corruption and handle occlusions with a sort function [23].Wang presented a multiresolution RS image segmentation method by combining a rainfall watershed algorithm and a fast region merging method, yielding a stable performance for high spatial resolution RS imagery, [24].
Our goal in this research is to develop an efficient segmentation method by combining MST and MHR, which should have the following properties: (a) be suitable for a variety of RS imagery, such as hyperspectral, multispectral, SAR, etc.; (b) be applicable to complex scenes (residential/industrial, residential/agriculture); (c) be efficient for large datasets; (d) the image objects should represent the sub-objects of geo-objects, whereby in an ideal situation, image objects are representations of geo-objects.The initial research was published in the 2015 International Workshop on Image and Data conference [25].In this paper, we extend the 2015 work significantly by presenting the method in detail and performing a comprehensive evaluation.The main parts of this study are as follows: (a) the MST algorithm for initial segmentation is proposed to be combined with the MHR algorithm for object merging which is adopted in FNEA; (b) an efficient implementation of the segmentation strategy is presented using data partition and a "reverse searching-forward processing" chain based on message passing interface (MPI) parallel technology; (c) high-resolution optical images, a hyperspectral image and a SAR image from four test sites are used to assess the proposed method.

Multi-Scale Segmentation Method Based on MST and MHR
The flowchart of our multi-scale segmentation method is shown in Figure 1.It consists of two main procedures: an initial segmentation using the MST algorithm based on the graph theory described in Section 2.1, and an object merging step using the MHR algorithm based on the FNEA method, detailed in Section 2.2.

Initial Segmentation Based on MST
Graph-based image segmentation techniques generally represent the problem in terms of a graph G = (V, E), where each node v i ∈ V corresponds to a pixel in the image, and each edge (v i , v j ) ∈ E connects a pair of neighboring pixels.Each edge is associated with a corresponding weight w((v i , v j )), which is a non-negative measure of the dissimilarity (e.g., the difference in intensity, color, motion, location or some other local attribute) between neighboring elements v i and v j [14].A segmentation S is a partition of V into components such that each component (or region) C ∈ V corresponds to a connected component in a graph G = (V, E ), where E ⊆ E. The pixel-based image is shown as a graph in Figure 2. The merging criterion D presented by Felzenszwalb and Huttenlocher [14] based on the MST algorithm is used to merge neighboring components, and is, The difference between two components C i , C j (C i , C j ∈ V), Di f (C i , C j ), is the minimum weight edge connecting the two components.That is, If there is no edge connecting C i and C j , Di f (C i , C j ) = 1.
The internal difference of a component C (C ∈ V), Int(C), is described as the largest weight in the MST of the component, MST(C, E), That is, The minimum internal difference, MInt(C i , C j ), is defined as, The threshold function is defined as, where |C| denotes the size of C and k is a constant parameter.The role of this parameter is to avoid small regions as excessive noise and to suppress the generation of a large connected region.
There are four steps in the initial segmentation procedure [14]: (1) A graph G = (V, E) with vertices V and edges E is built and the weight of each edge is the Euclidean distance d(v i , v j ) between neighboring elements v i and v j in terms of the intensities of all the bands.( 2) Sort E by non-decreasing edge weight.
(3) Start with an initial segmentation S 0 , where each vertex v i is in its own component.Compute a threshold function τ(C i ) for each component using Formula (5).Then Compute Di f (C i , C j ) and MInt(C i , C j ) for each component using Formulas (2) and (4).Subsequently we decide whether Di f (C i , C j ) <= MInt(C i , C j ).If the condition holds, the two components are merged; otherwise nothing would be done.Repeat the above steps until all the components are computed.(4) The output is a segmentation of V into components S = (C 1 , ..., C R ).

Object Merging Based on MHR
To implement the multi-scale segmentation, the MHR is introduced to merge two adjacent regions from the initial segmentation, which could reduce the disturbance from noise, reduce the fragmentized degree of object boundary, and obtain more regular objects [23].Not only the color heterogeneity (h color ) but also shape heterogeneity (h shape ) is considered.The heterogeneity (h) of MHR is defined as, h = w color h color + w shape h shape (6) where, h shape = w compt h compt + w smooth h smooth (8) The symbols w color , w shape , w compt , w smooth are weights of color, shape, compact and smooth, respectively, and The symbols n merge , n obj_A , n obj_B represent the number of objects after merging, the number of object A before merging, and the number of object B before merging.Other symbols l obj_A and l obj_B are denoted as the perimeters of object A and object B, respectively, where b obj_A and b obj_B are the bounding box perimeters of object A and object B. Lastly, the meaning of b merge is the bounding box perimeter of objects after merging.
There are three steps in the merging procedure: (1) Set the parameters of MHR, such as weights of color, shape, compact, and smooth (w color , w shape , w compt , w smooth ), and the scale parameter Q.Then compute the heterogeneity h of adjacent objects using Formula ( 6). ( 2) Decide whether h satisfies MHR.If h < (Q) 2 , the adjacent smaller objects are merged into a bigger one.Meanwhile, the average size, standard deviation, and mean of all the objects are calculated.Repeat this process until all the objects are merged.(3) Repeat steps 1-2 to accomplish the multi-scale segmentations.

Parallel Segmentation Based on MPI
One problem with applying this method is that it is very time-consuming for large datasets.Thus, in this section, we present a parallel implementation which relies on a data partition strategy and a "reverse searching-forward processing" chain strategy to boost the segmentation speed on a cluster platform using MPI parallel technology.The MPI has the advantages of portability, standardization, and efficiency, which could solve the communication problems of parallel processing [26].The computing resources of the cluster platform are divided into one master node and many slave nodes with the Master/Slave mode, which can handle dynamic load balancing [26].

Data Partition Strategy
Data partition is a common strategy for parallel processing [27]; the data is divided into rectangular blocks using the expand buffer strategy by the master node.The data partition strategy is shown in Figure 3.It has the following three steps: (1) The data is divided into rectangular blocks, and every block is labelled with "1, 2, 3, . . ., i, . . ." from left to right, up to down.(2) Each block is assigned to a different slave node; the block is labelled as i, that is, i = p + n × P where, (n, p) is the coordinate of the slave node, which is labeled with index i, and P is the number of slave nodes in the column direction.Here the column and row directions are equivalent to the coordinate in the x and y directions, respectively.
The slave nodes process their blocks according to the label of the blocks and send the result to the master node until all the data has been processed.

"Reverse Searching-Forward Processing" Chain
The "reverse searching-forward processing" chain is an interconnected system, which describes a series of steps.Each step can be performed simultaneously by different slave nodes [28].The chain is shown in Figure 4.At first, the slave nodes start to search for the relative data and parameters from the result until they find all the data, parameters, and algorithms.Then the slave nodes use the relative data, parameters, and algorithms to process from the source to the result until all the steps are completed.

Parallel Segmentation Method
The computing resources of the cluster platform are divided into one master node and many slave nodes.The master node is responsible for scheduling tasks and collecting results from the slave nodes.The slave nodes are responsible for processing and submitting the tasks.The source file is a series of images to be processed, and the target file is the segmentation result.The flowchart for the parallel segmentation method based on MPI is shown in Figure 5.The three steps in the parallel segmentation method are: (1) The master node is responsible for reading the image and dividing it into blocks using extended buffer strategy.A regular data division is used to divide the original image data into sub-rectangular data blocks, which will be used as the input data and assigned to the slave nodes for parallel computing.(2) The slave nodes are responsible for receiving blocks.First, each block is segmented separately, and the data block is segmented to obtain objects initially based on MST.Then, the objects are merged based on MHR to obtain the results.In the end, all the segmentation results are sent to the master node.
(3) The master node collects the results from each slave node, and outputs the segmentation image.
We apply two strategies to deal with the objects adjacent to borders.One is to use an extended buffer strategy to extend the size of the block under consideration to obtain an extended block.The extended block is processed, and the normal size of the block is output (Figure 6).The other strategy is to remove the boundary of the vectorized segmentation results by the mean feature of the neighboring objects (Figure 7); the mean is calculated using statistical methods.These two strategies do not deal with the boundary problem exhaustively.To solve this problem entirely would be a research question in itself and the subject of another study.Our two-fold solution can deal with the vast majority of specific problems and can be a trade-off between perfection and effort.

Test Sites
We tested the proposed method on a variety of images with different selected landscapes (Table 1).The test images were chosen with diverse scene complexities for exploring the general applicability of the proposed method.
Test data 1 (T1) covering the city of Potsdam in Germany was selected from the ISPRS benchmark project [29], which is an airborne true orthophoto with 0.05 m resolution and four bands (red, green, blue, and infrared).Test data 2 (T2) covering the city of DaFeng in China is a hyperspectral SPECIM AISA EAGLE II data acquired in November 2014 with a spatial resolution of 0.78 m, a spectral resolution of 9.6 nm, and 64 bands with a wavelength from 400 nm to 970 nm.The image was atmospherically corrected using the FLAASH model implemented in the ENVI software (http://www.harrisgeospatial.com/), and then 10 bands were selected using the minimum noise fraction (MNF) method from ENVI.Test data 3 (T3), covering the city of Lintong in the northwest of China, is a panchromatic (Pan) WorldView-2 image with a resolution of 0.5 m and a multispectral (MS) WorldView-2 image acquired in July 2011 with a resolution of 2.0 m (with eight bands, namely coastal, blue, green, yellow, red, red edge, near-infrared 1, and near-infrared 2).The WorldView-2 MS image and the Pan image were then fused using the PanSharp fusion method as implemented within the PCI Geomatica software (http://www.pcigeomatics.com/).Test data 4 (T4) covering the city of Genhe in China is a quad polarization RADARSAT-2 data acquired in September 2012 with a spatial resolution of 8.0 m.It was multi-looked using SARScape (http://www.harrisgeospatial.com/) and filtered using the enhanced Frost method provided by ENVI.

Results
Our segmentation method was evaluated and compared with the FNEA method, which has the same three parameters, namely, scale, weight of color, and weight of compactness.It is important to set these parameters properly.There are some methods for a fact-based determination of these parameters, such as estimation of scale parameters (ESP) [30], optimized image segmentation [31], segmentation parameter tuner (SPT) [32], plateau objective function (POF) [33] or the work of Stumpf and Kerle (2011), who optimized segmentation through the optimal use of the derived object features in a random forest framework [34].In FNEA, the color and the shape parameter work contrarily: The larger the weight of the color is, the better the spectral consistency of the resulting objects.The smaller the compactness weight is, the more complex the average shapes of the resulting objects are.In this study, the selection of image segmentation parameters is based on an iterative trial-and-error approach that is often utilized in object-based classification [35,36].
It should be noted that the scale parameter (Q) is the most important factor for segmentation as it controls the relative size of the resulting image objects and has a direct effect on the segmentation accuracy.The bigger the scale value is, the larger the object.The weight of color (w color ) is considered to be bigger than the weight of shape (w shape ), and the weight of smoothness (w smooth ) is considered to be similar to the weight of compactness (w compt ) [37].We set the same parameters (w color , w smooth ) for different experiments to compare different methods under the same conditions.For all four test data, w color is 0.9, w smooth is 0.5.The Q of the proposed method is the same as for the FNEA method.The segmentation parameters for the four experiments are shown in Table 2.The segmentation results of the proposed method and FNEA with Q of 240 for the first test data are shown in Figure 8. Visually, the image objects partially represent parts or sub-objects of the targeted geo-objects, such as buildings, trees, cars, and roads.The two segmentation methods were compared with manually interpreted buildings and individual trees.The segmentation results of the proposed method and the FNEA with Q of 360 for the second test data are shown in Figure 9. Buildings, waters, roads are well segmented.The two segmentation methods were compared to manually interpreted buildings as well.The third test was conducted in a residential/agriculture area.Figure 10 depicts the proposed method and the FNEA method with Q values of 260.As can be seen, buildings, roads and water bodies are well segmented.The two segmentation methods were also compared to manually interpreted objects.The segmentation results of the proposed method and FNEA with Q values of 2000 for the fourth test data are shown in Figure 11.The image objects partially represent sub-objects of geo-objects, such as water bodies, fields, and trees.The two results were compared to manually interpreted fields.Overall, the segmentation results of the two methods for the four experiments look similarly, but there are some differences such as those highlighted for areas A and B in Figures 8-11.For area A, highly detailed segments are found in Figure 9a

Accuracy Evaluation
The geometries of the image objects are compared through a supervised segmentation evaluation [38].This includes methods that assess the geometric differences between the generated image objects and the ground truth data [39].For this study, the SPT tool is used to evaluate the segmentation results with seven metrics, namely Hoover Index (H), Area-Fit-Index (AFI), Shape Index (SI), Rand Index (RI), F-measure (F), Segmentation Covering (C), Reference Bounded Segments Booster (RBSB) [32], as shown in Table 3.The ground truth data for the four experiments are shown in Table 4.The values of the respective metrics for the four experiments are shown in Table 5.
Table 3.The metrics for accuracy evaluation.

Metrics Formula Explanations
Hoover Index (H) Measures the number of correct detection based on the percentage of overlap between segmentation and reference ground truth (GT) [40].
CD is the number of correct detections and N GT represents the number of segments in the GT image.Range [0,1], "H = 0" stands for perfect segmentation.
Area-Fit-Index (AFI) Addresses over-/under-segmentation by analyzing the overlapping area between segmentation and reference GT [41].
where A k is the area, in pixels, of a reference segment C k in the GT image, and A l.i.k is the area, in pixels, of the segment in the segmentation outcome S, with the largest intersection with the reference segment C k .N GT is the number of segments in the GT image."AFI = 0" stands for perfect overlap.
Shape Index (SI) Addresses the shape conformity between segmentation and reference GT [42].
where N GT is the number of segments in the GT image, ρ i and ρ j are the perimeters of the segments C i and C j , and A i and A j are their respective areas.Range [0,1], "SI = 0" stands for perfect segmentation.
Rand Index (RI) Measures the ratio between pairs of pixels that were correctly classified and the total pairs of pixels [43].
Let I = {p 1 , . . .,p N } be the set of pixels of the original image and consider the set of all pairs of pixels P = {(p i ,p j ) I × I|i < j}.Moreover, C i , a segment in the segmentation S, and C j , a segment in the GT image, are considered as partitions of the image I.Then, P is divided into four different sets depending on where a pair (p i ,p j ) of pixels falls: P 11 : in the same segment both in C i and C j .P 10 : in the same segment in C i but different in C j .P 01 : in the same segment in C j but different in C i .P 00 : in different segments both C i and C j .Range [0,1], "RI = 0" stands for perfect segmentation.
Precision-Recall (F) Measures the trade-off between Precision and Recall considering segmentation as a classification process [44].Given a segment from the segmentation outcome S and a segment from its GT, four different regions can easily be differentiated: True positives (tp): pixels that belong to both S and GT.False positives (fp): pixels that belong to S but not to GT. False negatives (fn): pixels that belong to GT but not to S. True negatives (tn): pixels that do not belong to S or GT.Range [0,1], "F = 0" stands for perfect segmentation.

Metrics Formula Explanations
Segmentation Covering (C) Measures the number of pixels of the intersection of two segments [44].
where ΣN GT is the total number of pixels in the original image.
The overlap between two segments, C i in a segmentation S and C j in its GT, is defined as O(C i , C j ) = |Ci∩Cj| |Ci∪Cj| Range [0,1], "C = 0" stands for perfect segmentation.

Reference Bounded Segments Booster (RBSB)
Measures the ratio between the number of pixels outside the intersection of two segments with the area of the reference GT [45].where t represents a segment from GT and N GT is the number of segments in the GT image.fn, fp, tp is the same as Precision-Recall (F).Range [0,1], "RBSB = 0" stands for perfect segmentation.The values of metrics for Airborne, SPECIM AISA EAGLE II, WorldView-2, RADARSAT-2 image segmentation are shown in Figure 12.
All the values except AFI for the airborne image, SPECIM AISA EAGLE II image, and WorldView-2 image were relatively low.The lower the value, the better the result.The results shown in Table 4 demonstrate that the image objects do not over-estimate the ground truth polygons too much, which is desirable.However, the relatively high values suggest that the ground truth polygons are usually far larger in size compared to the evaluated image objects.For the RADARSAT-2 image, all the values are high except H and SI, but the ground truth polygons seem to match the image objects.We segmented the RADARSAT-2 image with different scale parameters and evaluated the segmentation results with the seven metrics.The object boundaries are not clear due to the speckle noise of SAR images, resulting in segments with rough boundaries compared to other images, so the values of RI, F, C, RBSB are close to 1, the AFI is very high, and the overlap is not good between segmentation and reference ground truth.Still, this radar specific phenomenon does not affect the comparative analysis as we use the same evaluation metrics and ground truth.The curve of the proposed method for the airborne image is the lowest, as shown in Figure 12a.The values of four metrics (AFI, RI, F, RBSB) of the proposed method are the lowest, which shows that the proposed method is better than the FNEA method for airborne images.There is a slight difference between the proposed method and the FNEA for the SPECIM AISA EAGLE II image, as shown in Figure 12b.The values of four metrics (SI, RI, F, C) of the proposed method are the same as the FNEA method, which indicates that the two methods perform similarly on the SPECIM AISA EAGLE II image.The curve of the proposed method is slightly higher than the FNEA method for the WorldView-2 image as shown in Figure 12c, while there is a little difference for the RADARSAT-2 image as shown in Figure 12d.The values of three metrics (AFI, RI, F) are slightly smaller using the proposed method compared to the FNEA, which indicates that the proposed method is better than the FNEA method for RADARSAT-2 image.In summary, the proposed method is generally applicable to different RS imagery and is at least comparable with the FNEA method in accuracy, or even slightly better.

Speed Evaluation
In this section, we analyze the segmentation speed using three metrics: speed-up, efficiency, and effectiveness.The speed evaluation metrics for image segmentation are shown in Table 6.
We use two WorldView-2 images with a volume of 1.1 GB and 4.3 GB to test the speed of our segmentation method.The segmentation efficiency is related to the segmentation parameters.The larger the block is, the faster the segmentation.The larger the scale is, the faster the segmentation.In the two tests, the block size is 1024 pixels, and the segmentation parameters are 200, 0.9, 0.5, respectively.The experimental processor is Intel(R) Xeon (R), CPU E5-2640, 24 CPU@2.50GHz with 16 GB of DDR3.The speed evaluation of parallel segmentation is shown in Figure 13.

Metrics Formula Remarks
Speed-up S P = T S T P T S is the processing time of single computing node.T P is the processing time of P computing node.
Efficiency S P is the speed-up of P computing node.
Effectiveness F P = S P P×T P = E P T P E P is the efficiency of P computing node.For the two test data sets, the running time (shown in Figure 13a) tends to decrease with the number of CPUs until it reaches 14, then it becomes stable and shows little variation because of the I/O (Input/Output) reaches a bottleneck.The speed-up (shown in Figure 13b) tends to increase with the number of CPUs until it reaches 14, at which point the maximum is 6 for 1.1 GB data, and 7.65 for 4.3 GB data, then it becomes stable.The efficiency (shown in Figure 13c) tends to decrease as the number of CPUs increases.The effectiveness (shown in Figure 13d) tends to increase with the number of CPUs until it reaches 11, after which it tends to decrease.In summary, when 14 CPUs are used for parallel computing, it is able to achieve the maximum speed-up for the two test data, when 11 CPUs are used for 1.1 GB data and 13 CPUs are used for 4.3 GB data, it is able to achieve the maximum effectiveness, and the parallel computing strategy is more efficient than stand-alone computing.The larger the dataset, the higher the speed-up is, and the higher the efficiency is.

Conclusions
This study proposes an efficient parallel multi-scale segmentation method for RS imagery based on graph theory and FNEA, which fully exploits the advantages of both the graph theory and the FNEA method.There are two main contributions: (a) the MST algorithm for initial segmentation is combined with the MHR algorithm which is used for object merging in FNEA; (b) an efficient implementation of the segmentation strategy is presented using data partition and the "reverse searching-forward processing" chain based on MPI parallel technology.Segmentation results using images from different sensors (Airborne, SPECIM AISA EAGLE II, WorldView-2, RADARSAT-2) and complex scenes (residential/industrial, residential/agriculture) covering four sites indicated a general applicability to different RS imagery (airborne optical, satellite optical, SAR, high-spectral) in accuracy using seven different metrics, namely Hoover Index, Area-Fit-Index, F-measure, Shape Index, Reference Bounded Segments Booster, Rand Index, Segmentation Covering.Segmentation results of the WorldView-2 large datasets demonstrate a high efficiency in speed in terms of three metrics: speed-up, efficiency, and effectiveness.
The proposed method has the following advantages: (a) it is an effective method using a fast graph segmentation algorithm to create the initial objects, yielding a linear complexity for the number of image pixels; (b) it is a multi-scale segmentation method for RS imagery using the MHR algorithm to merge the initial object, and it is suitable for a variety of landscapes (residential/industrial, residential/agriculture); (c) it is efficient for large datasets using data partition and the "reverse searching-forward processing" chain based on MPI parallel technology.We conclude that the proposed method has a general applicability and high efficiency for a variety of RS imagery (airborne optical, satellite optical, SAR, high-spectral) with a variety of landscapes (residential/industrial, residential/agriculture), and is at least comparable to the FNEA method in terms of accuracy, in some respects it even performs slightly better.
Nevertheless, as segmentation itself is a problem without a perfect solution-some scholars even claim it is an ill-posed problem-it is hard to find optimal segmentation parameters and there is no optimal solution across a variety of landscapes.Different landscapes may contain inherent features at different scales, thus requiring multi-scale consideration, but recent developments allow for the use of more objective comparisons between the various parameters and thresholds.We referred to the popular ESP tool [30] and to the work of Stumpf and Kerle who determined multiple optimal scales objectively and subsequently linked them to landscape fragments for image-based identification of landslides [34].Ma et al. [46] extensively analyzed feature selection methods for OBIA classifiers.These and other recent efforts pave the way for future workflows.The currently prevailing strategy of selecting segmentation parameters by an iterative trial-and-error approach still requires a lot of effort.Thus, automated methods need to be investigated for a variety of landscapes.In addition, the boundary problem of parallel segmentation could be tackled in this study satisfactorily but not perfectly.A full exploitation of this problem will be worth a study on its own.In fact, future work needs to integrate concepts from the OBIA community as laid out in [2] and concepts from computer vision, and the segmentation method will be implemented in an open source and made available to the whole RS community.

Figure 1 .
Figure 1.Flowchart of the multi-scale segmentation method based on MST and MHR.

Figure 2 .
Figure 2. Illustration of a graph-based image segmentation.

Figure 3 .
Figure 3. Data partition strategy based on MPI.

Figure 5 .
Figure 5. Flowchart for parallel segmentation based on MPI.
compared with Figure 9b.The opposite is true for area B. The same holds true for Figure 10a,b and for Figure 11a,b.

Figure 8 .
Figure 8. Segmentation results of the airborne true orthoimage at Q 240 (yellow outlines) and ground truth (104 buildings and individual trees, blue outlines), for area A, highly detailed segments are found in (a) compared with (b), the opposite is true for area B.

Figure 9 .
Figure 9. Segmentation results of SPECIM AISA EAGLE II image at Q 360 (yellow outlines) and ground truth (28 big buildings, blue outlines), for area A, highly detailed segments are found in (a) compared with (b), the opposite is true for area B.

Figure 10 .
Figure 10.Segmentation results of worldview-2 image at Q 260 (yellow outlines) and ground truth (28 fields, blue outlines), for area A and B, highly detailed segments are found in (b) compared with (a).

Figure 11 .
Figure 11.Segmentation results of RADARSAT-2 image at Q 2000 (yellow outlines) and ground truth (21 fields, blue outlines), for area A, highly detailed segments are found in (a) compared with (b), the opposite is true for area B.

Table 1 .
Summary of four test sites and imagery types.

Table 2 .
The segmentation parameters for four experiments.

Table 4 .
The ground truth data for the four experiments.

Table 5 .
The values of fitness functions for four experiments.

Table 6 .
Speed evaluation metrics for image segmentation.