A method for hierarchical weighted fitting of regular grid DSM with discrete points

A Digital Surface Model (DSM) is a crucial spatial geographic information data used to describe the shape of the earth’s surface in Geographic Information Systems (GIS). DSM is the core data used in terrain analysis in GIS. A regular grid DSM is generally generated by interpolating a large number of discrete point clouds. This paper proposes a method of using a hierarchical weighted strategy to fit a regular grid DSM with discrete points. This method uses a pyramid hierarchical strategy to refine the target regular grid from one grid with finer parameters of 3*3, until the nth level (the interval of the grid is equal to the expected interval), and then gradually places the discrete point cloud into the corresponding grid by weighted averaging, and uses the result of this level as the initial value of the next level. This algorithm can avoid the problem of low efficiency in retrieving a large number of discrete point clouds, and the indirect interpolation method not considering the contribution of distant neighboring point clouds. The operation of point cloud data is a stream operation, which does not require consideration of the topological information of point clouds, and has simple operation and no additional memory consumption. It is especially suitable for the production of regular grid DSM with massive point clouds. To verify the effectiveness of this method, the article selected six typical terrain data such as high mountains, mountains, hills, plains, urban areas, and lakes for experiments. The results show that compared with the construct-TIN method for producing DSM, this method has very good processing accuracy and processing efficiency.


Introduction
Digital Surface Model (DSM) is an important spatial geographic information data in Geographic Information Systems (GIS) and is the core data for terrain analysis in GIS.It has a wide range of applications in surveying and mapping, geomorphology and geology, engineering construction, military, and other fields (Tang, 2014).DSM mainly includes two forms: regular rectangular grid and Triangulated Irregular Network (TIN) (PETZOLD B, 1999).Among them, the regular rectangular grid DSM is more widely used in photogrammetry due to its simple data structure, easy algorithm implementation, and advantages in spatial operations and storage (Razak, 2011).
In photogrammetry, regular grid DEMs are mainly generated by interpolating massive and dense matching point clouds or laser point clouds.In actual processing, since large-scale terrain is complex and has large local differences, it is impossible to use a polynomial to fit the entire terrain.Because the fitting accuracy of low-degree polynomials is inevitably poor, and high-degree polynomials may produce unstable solutions.Therefore, global function interpolation is generally not used in DSM interpolation, but local function interpolation is used (Pan, 2023).Local function interpolation divides the entire area into several blocks, and fits each block with different functions, while considering the continuity between adjacent block functions.For the surface that is not smooth or even has a fracture line, further subdivision is required, and smooth or even continuous conditions cannot be used.In addition, point-bypoint interpolation methods are widely applied, such as inverse ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume X-1-2024 ISPRS TC I Mid-term Symposium "Intelligent Sensing and Remote Sensing Application", 13-17 May 2024, Changsha, China distance weighting (IDW) method (D Shepard, 1968), moving surface fitting method (SChut, 1976), finite element interpolation method (Ebner, 1978), etc.These algorithms define a local function to fit a certain number of discrete point clouds around the target regular grid point.These methods are flexible and effective, have high accuracy under general conditions, and have simple calculation methods and do not require a large amount of computer memory (Pan, 2023).
However, the disadvantages are also obvious.Since it is necessary to search for discrete points around each grid point, the calculation efficiency is extremely low.Especially in largescale DSM update tasks, it is difficult to meet the requirements.
Therefore, some algorithms have been developed to establish block indexes for discrete point clouds, which can improve search efficiency to a certain extent (Chen et al., 2008;Zhang, 2011;Yang, 2013;Chen et al., 2016).However, the calculation process is complex, the complex linked storage structure increases the difficulty of algorithm maintenance, and it is still limited by the size of the survey area and the number of point clouds.
In addition, there is another classic indirect algorithm.This type of algorithm first uses point clouds to construct Triangular Irregular Networks (TIN), and then interpolates the elevation of a regular grid inside the triangle (Dai, 2018).Although indirect algorithms reduce the number of point cloud retrievals and improve efficiency compared to direct interpolation algorithms, the generation of TIN and the determination of the triangle in which the point is located still require a long time.In particular, constructing TIN requires a large amount of computing memory, so indirect algorithms are limited by the number of point clouds.
Moreover, since the grid points are only interpolated from the three vertices of the triangle in which they are located, the interpolation accuracy of indirect algorithms is more likely to be affected by a triangle vertex with a large elevation error.
In summary, most of the current image matching direct interpolation algorithms start from ordered target regular grid points and search in massive unordered discrete point clouds.
Therefore, the efficiency of searching for discrete point clouds around regular grid points is greatly limited.Meanwhile, for some indirect algorithms, the elevation of regular grid points is only interpolated from the elevations of the three vertices of the triangle in which they are located, so their accuracy is limited.
To address these issues, we propose a method for fitting regular grid DSM with discrete points by a hierarchical weighted strategy.This method fully considers the characteristics of unordered discrete point clouds, only traverses the point clouds in order, places the point clouds into the corresponding grid according to the geographic coordinates, integrates the contributions of neighboring point clouds to the elevation of the grid point through a pyramid structure, and constrains the influence of distant discrete point clouds on the accuracy of the regular grid point through a weighted strategy.Finally, this method achieves fast production of regular grid DSM from discrete point clouds.

Methodology
The hierarchical weighted fitting method shares fundamental principles with the inverse distance weighted fitting method, both being rapid distance-weighted fitting algorithms.The basic idea is shown in Figure 1.
The core of this algorithm is to hierarchically place point clouds into a pyramid-structured grid for create DSM from point cloud data.When the point cloud is placed in the first-level grid, the resulting elevation is the average elevation of all point clouds because the grid has only one cell.The second-level grid is based on the grid obtained in the first-level.The grid is first refined (the recommended refinement parameter is 3*3), and then the point cloud is re-placed into the second-level grid, and the initial value in the grid is weighted average (the weight of the initial value in the grid is recommended to be 0.5) to obtain the elevation value of the second-level grid.Then continue to refine the second-level grid.Iterate this process until the interval of the nth-level grid is equal to the expected interval.The specific processing steps are as follows: Step 1: Establish a pyramid-structure grid.The first-level grid only contains one cell.The first-level grid is divided according to the refinement parameter of 3 * 3 in the X and Y directions to obtain the second-level grid.Continue to refine until the interval of the Nth-level grid is equal to the expected interval (i.e., the resolution of the DSM).The central positions of cells in the Nth-level grid corresponds to the regular grid points of the target DSM.
Step 2: Starting from the first-level grid, place discrete points into the grid according to equal weight.Since there is only one Step 3: Divide the first-level grid according to the refinement parameter P*P to obtain the second-level grid.The second-level grid takes the elevation of the first-level grid as the initial value of the grid elevation.Then, the discrete point cloud is re-placed into the second-level grid and weighted averaged with the initial elevation in the grid (the weight of the initial elevation in the grid is recommended to be 0.5), to obtain the new second-level grid.
Step 4: Continue to refine the grid and re-place all discrete points in the grid to form a new grid.Iterate this process until the Nth-level grid is calculated.
Step 5: At this point, the elevation value of each cell in the Nth-level grid is the corresponding grid point elevation of the target DSM.
In the algorithm, the final interpolation value of the elevation of a target grid point is: where = the target grid point's elevation; ̃ = the contribution value of the j-th discrete point to the grid point; Z  = the elevation value of the j-th discrete point; represents that from the first-level until the K-level, the j-th discrete point remains in the grid; = the total number of discrete points in the grid Level 1 Level 2 Level n Level n-1 where the grid point is located at the i-th level grid.
So it can be seen that all point clouds contribute to the elevation of the grid point.The closer a discrete point is to the grid point, the higher its contribution value, and vice versa, the lower its contribution value, even approaching zero.

Experimental Data
To verify the effectiveness of the proposed method, the paper selected six typical terrain point clouds, including high mountains, mountains, hills, plains, urban areas, and lakes, for experiments.The information of each data is shown in Table 1. The

Experimental Scheme
To verify the efficiency and fidelity of the method proposed in this paper, a comparative experiment was designed.Two methods were used: the hierarchical weighted method in this paper (hereinafter referred to as Grd) and the indirect method of constructing TIN to generate DSM (hereinafter referred to as Tin).Six typical terrain point cloud data were used to produce DSMs with the same resolution.Then, the two methods were compared based on two criteria: processing speed and fidelity.
Figure 3 shows the experimental scheme.
The evaluation criteria for the two methods are as follows: Criterion 1: Processing speed Method: Generate DSMs with the same resolution using the two methods for each set of point cloud data, and record the processing time for each method.Obviously, the smaller the difference, the higher the fidelity of the method used.As shown in Figure 4 and Figure 5, the Grd method and the Tin method have comparable effects in displaying the details of mountainous terrain and serpentine roads.In Figure 6, the DSMs generated by the two methods show comparable effects for houses.For the top of the tower, the Tin method performs better than the Grd method.This is because Tin contains more terrain feature points, which is more conducive to recovering the top of the tower.As shown in Figure 7, for flat terrain, both methods can preserve the geomorphic features well.In Figure 8, for the urban data set, since there are no buildings with pointed roofs, the Grd method and the Tin method produce similar effects.Comparing Figure 9 The first evaluation indicator is processing speed.The processing times of the two methods are shown in Table 2.By comparison, it can be seen that the Grd method in this paper takes about half the time of the Tin method.This is mainly because constructing TIN and determining the triangle where the grid point is located take more time.Using the method proposed in this paper can improve the efficiency of generating DSMs from discrete point clouds.

Data
The The second evaluation indicator is fidelity.Calculate dZs for all points in the six sets of point cloud data.And then the distribution of dZ values is statistically analyzed.Figure 10 shows the violin plot of the dZ distribution for six sets of point cloud data.The body of the violin plot shows the frequency distribution, and the wider it is, the higher the frequency near a certain value.From Figure 10, it can be seen that the dZ distributions obtained by the Grd method and Tin method are approximately the same.The fidelity of the six sets of point cloud data can be ranked from high to low as follows: lakes > plains > hills > mountains > high mountains > urban areas.
Table 3 lists the mean, variance, and percentages of dZ absolute values within 1m, between 1m and 5m, and above 5m for each set of point cloud data.According to Table 3, the means of dZ obtained using both methods are close to 0, the variances are similar, and most of the dZs are within 1m.Based on Figure 10 and Table 3, it can be concluded that the fidelity of the Grd method and Tin method is comparable.

Figure 1 .
Figure 1.The basic idea of the hierarchical weighted fitting method

Figure 2 .
Figure 2. Point cloud data for 6 typical terrains: (a) high mountains, (b) mountains, (c) hills, (d) plains, (e) urban areas, (f) lakes the XY coordinates of each point in the point cloud to interpolate  ′ in the resulting DSM.Calculate the difference dZ between  ′ and the point cloud coordinate Z.

Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .
(b) and Figure 9(d), it can be seen that the Grd method is superior to the Tin method.In the DSM generated by the Grd method, the surface of the lake in the upper left corner is smoother.At the bottom of Figure 9(d), there are many obvious triangles, and the Tin method does not show the original topography.Overall, the DSMs generated by two methods have comparable effects.The Tin method is more conducive to preserving sharply changing geomorphic features.The Grd method is more suitable for the production of DSMs in flat areas.Resulting DSMs of high mountains: (a) DSM (Grd), (b)Partial image of DSM(Grd), (c) DSM (Tin), (d) Partial image of DSM(Tin) Resulting DSMs of mountains: (a) DSM (Grd), (b)Partial image of DSM(Grd), (c) DSM (Tin), (d) Partial image of DSM(Tin) ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume X-1-2024 ISPRS TC I Mid-term Symposium "Intelligent Sensing and Remote Sensing Application", 13-17 May 2024, Changsha, China This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.https://doi.org/10.5194/isprs-annals-X-1-2024-91-2024| © Author(s) 2024.CC BY 4.0 License.Resulting DSMs of hills: (a) DSM (Grd), (b)Partial image of DSM(Grd), (c) DSM (Tin), (d) Partial image of DSM(Tin) Resulting DSMs of plains: (a) DSM (Grd), (b)Partial image of DSM(Grd), (c) DSM (Tin), (d) Partial image of DSM(Tin) Resulting DSMs of urban areas: (a) DSM (Grd), (b)Partial image of DSM(Grd), (c) DSM (Tin), (d) Partial image of DSM(Tin) Resulting DSMs of lakes: (a) DSM (Grd), (b)Partial image of DSM(Grd), (c) DSM (Tin), (d) Partial image of DSM(Tin)

Figure 10 .
Figure 10.The violin plot of the dZ distribution

Table 1 .
Experimental data information

Table 2 .
Process time of two methods

Table 3 .
Statistical table for the distribution of dZs4.ConclusionThis paper proposes a method for fitting regular grid DSM with discrete points by a hierarchical weighted strategy to address the low efficiency of existing technologies in retrieving massive discrete point clouds and the issue of indirect interpolation methods not considering the contributions of distant neighboring point clouds.The paper conducts experiments using six types of typical terrain data, including high mountains, mountains, hills, plains, urban areas, and lakes.The results demonstrate that the proposed method in this paper to fit regular grid DSM improves processing efficiency by more than twice compared to the indirect method of constructing TIN to generate DSM through interpolation, while maintaining comparable data fidelity.Therefore, we believe that the hierarchical weighted method in this paper has better prospects for applications.Additionally, if computing resources are limited, such as mobile devices, edge computing devices, etc., this method has better performance.