SEMI-GLOBAL MERGING OF DIGITAL SURFACE MODELS FROM MULTIPLE STEREOPAIRS

The semi-global optimization algorithm, which approximates a global 2D smoothness constraint by combining several 1D constraints, has been widely used in the field of image dense matching for digital surface model (DSM) generation. However, due to occlusion, shadow and textureless area of the matching images, some inconsistency may exist in the overlapping areas of different DSMs. To address this problem, based on the DSMs generated by semi-global matching from multiple stereopairs, a novel semiglobal merging algorithm is proposed to generate a reliable and consistent DSM in this paper. Two datasets, each covering 1km, are used to validate the proposed method. Experimental results show that the optimal DSM after merging can effectively eliminate the inconsistency and reduce redundancy in the overlapping areas. * Corresponding author


INTRODUCTION
Digital surface model (DSM) plays an important role in many applications, for example, object extraction and change detection. Thanks to the outstanding performance of semiglobal matching (SGM) (Hirschmüller, 2005;Hirschmüller, 2008), many researchers (Yastikli, et al., 2014;Ghuffar, 2016) tended to achieve the generation of DSMs with the SGM algorithm. However, considering that DSMs (also referred to as point cloud data) generated by SGM are single stereopair-based, as shown in Figure 1, certain inconsistencies may exist in the overlapping areas of different stereopairs. Furthermore, DSMs generated by SGM are of high density and redundancy, especially for the aerial images with multiple overlap. Thus, merging of these DSMs to obtain a consistent and accurate DSM is very necessary. In the DSM merging, many scholars (Fratarcangeli, et al., 2016;Jaud, et al., 2016) tend to achieve it by using commercial software (e.g., INPHO, PhotoScan, MicMac and SURE). Meanwhile, some researchers cast the DSM merging as surface reconstruction problem, and screened poisson (Kazhdan and Hoppe, 2013), floating scale (Fuhrmann and Goesele, 2014), voronoi filtering (Amenta and Bern, 1999) and some probabilistic methods (Agrawal and Davis, 2001) were used to achieve a good surface reconstruction. In addition, Sadeq, et al. (2016) proposed a Bayesian approach to merge different DSMs from different sources. Furthermore, some researchers try to first divide the point cloud data into grids with regular intervals, then select an optimal point at each grid with the "winner takes all" strategy. For the images with good quality, the generated DSM is good. However, for the images with poor quality, considering that there is no consideration of spatial correlation, the generated DSM may still contain some outliers. And as Gong and Fritsch (2016) point out that merging of DSMs from multiple stereopairs can effectively remove outliers and further improve the quality of DSM generated by the SGM algorithm. Thus, based on the DSMs from multiple stereopairs generated by SGM, an automatic semi-global merging algorithm is proposed to obtain a more consistent and accurate DSM in this paper. This merging algorithm can not only remove the outliers in the point cloud data to obtain a reliable and consistent DSM, but also effectively reduce the redundancy in the overlapping areas of multiple DSMs.
The rest of this paper is organized as follows. Section 2 first describes the basic idea of semi-global optimization. Section 3 presents the proposed merging algorithm with semi-global optimization for point cloud data from multi-stereopairs. Then, Section 4 displays the experimental results. Finally, conclusions are drawn in Section 5.

SEMI -GLOBAL OPTIMIZATION
Semi-global optimization is inspired by the SGM algorithm (Hirschmüller, 2008). Its basic idea is to approximate a global 2D smoothness constraint by combining several 1D constraints. It mainly consists of three parts: establishment of 3D cost matrix, multi-directional aggregation of costs, and acquisition of optimal surface, as shown in Figure 2. In the 3D cost matrix, the first two dimensions represent X and Y on the two -dimensional plane, and the third dimension represents the label. The values in each grid (x, y, l) of the 3D cost matrix represent the cost of selecting the label l on the 2D plane (x, y). In general, the smaller the value, the greater the probability or possibility of selecting the label l. Then 8 or 16 directions of the dynamic programming algorithm are used for multi-directional aggregation of costs. During the dynamic programming process, it only needs to record the minimum cost of selecting each label in consideration of the cost of the label and the cost of the smooth constraint, but not need to record the optimal path. After the dynamic programming, the minimum cost of each grid will be accumulated to the accumulated cost matrix. And an optimal surface is calculated based on this accumulated cost matrix after all directions of the dynamic planning is completed. The optimal surface can be determined by the "winner takes all" strategy, that is, label of each grid on the 2D plane with the smallest value in the accumulated cost matrix is selected as the optimal label. Finally, a median filter is used to reduce the noise in the optimal surface.

SEMI -GLOBAL OPTIMIZATION FOR MERGING OF MULTIPLE DSMS
Point cloud data obtained by SGM are single stereopair-based, and certain inconsistencies may exist in the overlap of different stereopairs. To obtain a consistent large range of DSM, a semiglobal optimization-based method is proposed to merge the point cloud data from multiple stereopairs in this paper, and details are as follows: Step 1: Divided the point cloud data into different blocks Considering that the amount of point cloud data may be too large to process at one time, the whole point cloud data is first divided into several DSMs with 1 km 2 area each. The division could also facilitate the subsequent processing (e.g., object extraction and change detection).
Step 2: Clustering of point cloud data The point cloud data in each kilometer is first assigned a grid index, and the grid size is set to a specific distance (e.g., 1 m). Then, the points in each grid cell are clustered into several clusters by a certain height distance (e.g., 0.9 m). The center and weight of each cluster are calculated. The center is represented by the mean of points in the cluster, and the weight is determined by the number of points in the cluster. The greater the number of points of the cluster, the greater the weight of each cluster. The process of clustering the points is shown in Figure 3. The benefits of the clustering mainly have two aspects. On the one hand, it can reduce the number of candidate labels in the subsequent semi-global optimization and increase the speed. On the other hand, it can guide the optimization algorithm to select those labels composed of dense points by increasing the weight of the labels.
Step 3 where l = possible cluster centers for all grid cells, l * = optimal cluster centers, V l = weight of the cluster obtained by step 2, dz = height difference between two neighbor centers, DZ 1 = threshold of height difference, DZ 2 = threshold of height difference.
In general, natural surface and artificial object surface are mostly continuous and smooth. The aim of defining DZ 1 is to make the generated DSM surface smooth. Meanwhile, for the discontinuities, e.g., breaklines, defining DZ 2 is to set a larger constant penalty. In this paper, DZ 1 and DZ 2 are set to 1.0 m and 2.0 m, respectively. After semi-global optimization, a median filter is used to reduce the noise in the generated DSM.
Step 4: Progressive TIN-based DSM densifying After achieving the grid DSM obtained by the above steps, the DSM can be densified by progressive TIN algorithm. In the densification process, a sparse TIN is first derived from the grid DSM, then the points are progressively added to the TIN if they are below the defined thresholds. More details about the progressive TIN algorithm can be seen in (Axelsson, 2000).

EXPERIMENTAL RESULTS
In this paper, two datasets with 1 km 2 area each are used for the experiments. Each dataset is composed of point cloud data from multiple stereopairs generated by the SGM algorithm (Hirschmüller, 2008). The stereopairs are from traditional aerial images obtained by DMC camera. The image size is 7,680 pixel × 13,824 pixel, pixel size is 12 um, the focal length is 120 mm, and forward lap and side lap are 65% and 35%, respectively. corresponding orthophotos of the two datasets are shown in Figure 4. From Figure 4, it can be seen that dataset 1 is a typical suburban area with sparse housing and dense farmland. Dataset 2 is a complex area; the right half of the dataset is a typical main urban area with very dense houses, and the left half is a mountain covered by dense trees.
To visually represent the merging results of point cloud data from multiple stereopairs in this study, comparisons before and after merging were selected to be enlarged for dataset 1 and dataset 2. For dataset 1, cross section and triangulation of the results are shown in Figures 5 and 6. For dataset 2, cross section and triangulation of the results are shown in Figures 7 and 8.    Figure 5 and 7 show that multilayered phenomenon in those building boundaries is significantly reduced after merging, and the accuracy of point cloud data is improved. Furthermore, from Figure 6 and 8, it can be seen that many spikes are effectively removed after merging, and a smoothed surface is generated.
To better show the results of the proposed method, median of each grid cell is used to compare with our semi-global merging results. Here, the size of grid cell is set to the same as the proposed method. Furthermore, considering that points composed of median of each grid cell are sparse, TIN algorithm is also used to add the points which are below the defined thresholds.
Comparisons are shown in Figure 9, where Area 1 and Area 2 are from dataset1, and Area 3 is from dataset 2. The first row is the triangulation of the raw data, the second row is the results from median and TIN algorithm, and the third row is the results of the proposed method. From Figure 9, it can be seen that the triangulation of the raw data includes many burr-like noises, and the results of the median and TIN algorithm are improved obviously, but there is still a small amount of burr-like noise.
The results from the proposed method are smooth and most of the burr-like noises have been further eliminated.

CONCLUSION
In this paper, a novel semi-global merging of DSMs from multiple stereopairs is proposed to remove the inconsistency in the point cloud data. And experimental results show that DSM after merging is more accurate and reliable, and its data redundancy is effectively reduced. However, considering the lack of reference datasets, e.g., corresponding LiDAR data, more quantitative evaluation is still unavailable, which is also our further research.