Integrating topographic features and patch matching into point cloud restoration for terrain modelling

ABSTRACT

Point clouds are widely used in Earth surface research but usually exhibit gaps of missing data. Previous point cloud restoration methods used in terrain modelling have not fully considered complex terrain characteristics, which can be summarised as the controlling role of topographic features in shaping terrain surfaces and the inherent similarities observed among these surfaces. This work introduces a novel method that integrates Topographic Features and Patch Matching (TFPM) into point cloud restoration processes for terrain modelling. The method mainly contains three steps. First, identifying gap boundary points. Second, topographic feature points are extracted and subsequently interpolated into the identified gaps. Third, searching other parts of the raw point cloud for patches resembling the gaps, and the identified patches are used as templates to restore the point cloud. The proposed method is benchmarked against three state-of-the-art point cloud restoration methods. The experimental results demonstrate that the TFPM method consistently exhibits superior accuracy in terrain modelling and analysis, as evidenced by low values of the root mean square error, average elevation difference, and average slope difference. This work endeavours to incorporate topographic features into point cloud restoration processes and can benefit future research related to terrain modelling and analysis.

KEYWORDS:

• Point clouds
• point cloud restoration
• topographic features
• patch matching
• terrain modelling

1. Introduction

Terrain models, such as digital elevation models (DEMs), have been extensively utilised to analyse and study surface processes on the Earth's terrain, which significantly influence various physical processes (Xiong et al. 2022; Zhao et al. 2023), particularly in relation to natural hazards such as landslides (Syzdykbayev, Karimi, and Karimi 2020), floods (Schumann 2014), and earthquakes (Okyay et al. 2019). With advancements in remote sensing acquisition techniques, acquiring high-density point cloud data, comprising three-dimensional coordinates and additional attributes such as colour or intensity, has become increasingly convenient (Xiong et al. 2022). These obtained point clouds have been broadly applied in infrastructure classification (Aljumaily et al. 2023; Mirzaei et al. 2022), object detection (Gharineiat, Kurdi, and Campbell 2022; Nahhas et al. 2018), 3D construction (Mirzaei et al. 2022), and terrain modelling (Callow, May, and Leopold 2018; Maguya, Junttila, and Kauranne 2013). Various operations, such as registration (Han, Zhang, and Zhang 2023), filtering (Dai et al. 2023a), denoising (Duan et al. 2021), restoration (Fu, Hu, and Guo 2018), and simplification (Chen et al. 2023b), are conducted with point clouds to enhance these applications. This study specifically concentrates on point cloud restoration for subsequent terrain modelling and analysis because the obtained point cloud data for terrain surfaces usually exhibit gaps of missing data, mainly due to the limitations of acquisition techniques, occlusion by vegetation or other obstructions, and the limited number of scans from different viewing directions (Dinesh, Bajic, and Cheung 2017; Rengers et al. 2016). For instance, in mountainous areas, complex elevation differences and terrain relief interfere with the signals emitted by sensors, resulting in gaps in ridge and valley areas (Boulton and Stokes 2018). As a result, errors and uncertainties tend to be amplified, especially in areas with complex geometry and significant variations in topography (Nourbakhshbeidokhti et al. 2019). Consequently, these factors introduce inaccuracies and uncertainties in terrain modelling, particularly in regions with complex topographic features. In addition, subsequent procedures such as filtering may also cause missing point cloud data (Podobnikar and Vrečko 2012), resulting in significant loss of topographic information (Li et al. 2022b), thus posing challenges for subsequent terrain analyses, such as drainage network analysis and terrain change detection (Lyu et al. 2021; Nourbakhshbeidokhti et al. 2019). Therefore, it is imperative to restore point clouds.

Many studies have focused on point cloud restoration, and the methods can be divided into mesh-based and point-based methods (Lin and Wang 2016). Mesh-based methods use 3D mesh models, such as triangular irregular networks (TINs), to reconstruct the target objects. Based on the mesh, it is easy to locate the missing region and obtain the local geometry and topological structure; then, this information is used to restore the missing area (Huang et al. 2022; Wang and Oliveira 2007). However, this approach highly depends on the quality of the generated mesh and can require considerable time and extensive computations due to the intricate process of constructing meshes from high-density point clouds (Huang et al. 2022). The second category involves directly performing point operations and can be divided into subcategories, including interpolation, patch matching, deep learning methods, etc. Specifically, interpolation methods, such as the inverse distance weighting (IDW), polynomial, kriging, and natural neighbour methods, approximate the surface elevation based on the elevation of the surrounding measured data, but they fail to restore terrain features and may not accurately reflect the actual surface conditions in regions of high topographic variation (Chen, Gao, and Devereux 2017). In contrast to the methods mentioned above, which only utilise the surrounding data for gap restoration and often result in restored surfaces that appear smoother than the actual surfaces, patch matching methods have the advantage of leveraging both local similarity and nonlocal self-similarity within point clouds (Hu, Fu, and Guo 2019). Gap restoration with the patch matching method is performed iteratively using templates near the gap boundaries to find the best matching regions elsewhere in the cloud, from which existing points are transferred to the gap (Dinesh, Bajic, and Cheung 2017; Fu, Hu, and Guo 2018; Hu, Fu, and Guo 2019). This method retrieves the most similar region for the gap area from elsewhere and is based on the similarity of different parts of the object to be restored (Fu, Hu, and Guo 2018). However, it is difficult to find a perfectly similar region for gaps in complex and irregular areas, such as terrain surfaces. Deep learning methods for point cloud restoration have been widely employed and extensively developed in computer graphics (Wen et al. 2020). However, their application to terrain surface point cloud restoration is relatively limited, with an initial focus on point cloud completion for 3D objects such as cars, rather than addressing gaps in surface reconstruction.

Although some methods mentioned above have yielded ideal point cloud restoration results, they were initially designed primarily for 3D point clouds in computer graphics, and their effectiveness on 2.5D terrain point clouds is projected to be low due to the lack of consideration of terrain characteristics (Chen et al. 2023a). Currently, interpolation techniques are mostly used to address point cloud gaps during terrain modelling, which introduces uncertainties and inaccuracies in areas with intricate terrains and does not result in the restoration of topographic features (Xiong et al. 2022) because these interpolation techniques are elevation based rather than topographic feature based (Chen et al. 2022). Previous studies have proven that the terrain skeleton is crucial in determining the overall structure and characteristics of the terrain (Li et al. 2022a). The absence of skeleton information can lead to the loss of the structural details of terrain features during terrain data processing (Li et al. 2022b). Consequently, when performing terrain point cloud restoration, particular attention should be given to topographic features. Currently, sharp feature lines have been employed to restore comparatively regular 3D objects (Wang et al. 2017). A frequently employed approach for interpolating feature points within gaps involves utilising the endpoints of feature lines for curve fitting purposes. Subsequently, the fitted curves are sampled to obtain interpolated feature points (Tang et al. 2017). This idea can also be extended to restore terrain point clouds through the interpolation of topographic feature points (representative points of topographic features) within gap areas. Moreover, there are notable similarities in terrain surface and landforms. For instance, the Loess Plateau primarily consists of comparable landforms, such as Yuan, Liang, and Mao, and karst landforms consist of Fenglin and Fengcong (Li et al. 2020). Even on a small scale, different regions of hillslopes and terraces exhibit certain similarities.

Currently, the classic terrain point cloud restoration methods fail to adequately consider topographic features and inherent similarities among these surfaces, which deeply limits the accuracy of terrain modelling and analysis. Therefore, a restoration method for point clouds that integrates Topographic Features and Patch Matching (TFPM) is proposed for terrain modelling. The patch matching technique is adopted to leverage the similarity of terrain surfaces, and the identified actual terrain surfaces are used for the restoration of point clouds, which ensures that the restored surfaces bear a stronger resemblance to natural terrain surfaces. Moreover, topographic features are utilised to enhance the capability of the patch matching method in identifying similar regions for gap areas within complex terrain surfaces. Furthermore, enhancements have been made to the classic gap boundary and topographic feature detection method to optimise its performance in support of subsequent restoration processes. A set of metrics are employed to assess the accuracy of the TFPM method by comparing the discrepancies between the restored points and the corresponding ground truth surface.

2. Materials and methods

Figure 1 illustrates the process of the proposed method, which consists of five main steps:

1. Identifying the boundaries of gaps in the point clouds to determine the areas that need to be restored (Figure 1(b));

2. Extracting topographic feature points and interpolating them into the gap areas, where the extracted topographic features may possibly pass over (Figure 1(c, d, e));

3. Segmenting large gaps into multiple smaller gaps based on the interpolated topographic feature points (Figure 1(e));

4. Dividing the raw point cloud data into multiple patches and conducting a global search to find the patch that best matches the segmented smaller gaps (Figure 1(f)).

5. The identified patch is used as a template to fill gaps with a smoothing operation.

Figure 1. Main stages of the TFPM method. (a) Point clouds with a gap. (b) Detected gap boundary points. (c) Extracted topographic feature points (ridge). (d) A curve fitted using the extracted topographic feature points within the gap. (e) Topographic feature point interpolation and the segmented gap. (f) Patch matching for each segmented target patch.

Some of the steps mentioned above are influenced by specific processing scales. They are defined relative to the Average Point Distance (APD) of the raw point cloud data (formula 1), and some of them are set as 3 times the APD, which is a commonly utilised buffer distance during point cloud processing (Zhou et al. 2018). Notably, the ‘topographic feature points’ mentioned below refer to representative points of topographic features, such as ridges, valleys, and terrace risers. (1) $\begin{array}{c}APD={\displaystyle \frac{{\sum }_{i=1}^n\left(d_i\right)}{n}}\end{array}$(1) where $n$ is the quantity of raw point cloud data and $d_i$ is the distance between point $i$ and its nearest point.

2.1. Gap boundary point detection

Identifying gap areas is the initial stage in the restoration process. In this study, we made improvements to the boundary point detection method proposed by Kai et al. (2020), which is based on the differences in the distribution of surrounding points between gap boundary points and nongap boundary points. For example, the central point $O$ in Figure 2(e) is a boundary point because there are no neighbour points of point $O$ distributed to the southeast. The surrounding points of central point p refer to points that are within its search radius $R_{gap}$, and this parameter is set as 6 times the APD based on detailed trial (section 4.2). However, this approach fails to differentiate between external boundary points and internal gap boundary points, resulting in significant computational overhead when applied to the entire point cloud dataset. Thus, a preprocessing step is incorporated to ignore external boundary points and locate potential gap boundary points, thereby significantly enhancing the efficiency of the detection process. Specifically, first, the entire area is gridded into raster cells (elevation direction) with a size of 3 times the APD. Cells that do not contain any point cloud are assigned a value of 0, while others are set to 1. The gridded raster cells are then traversed row by row. For each row, the first encountered nonzero cell is regarded as the left boundary, and the last detected nonzero cell is considered the right boundary. The cells to the left of the left boundary and to the right of the right boundary are identified as empty areas, while cells with a value of 0 in between are classified as gap areas. The nearest nonzero value cells on the left and right sides of the gap areas are designated as potential gap boundaries (Figure 2(b)), which contain potential gap boundary points. After completing the row-wise traversal, the same strategy is applied for column-wise traversal. Ultimately, this strategy enables the identification of potential gap boundary points while excluding external boundaries. Then the traditional method is applied to the potential gap boundary points rather than the entire point cloud to obtain the final gap boundary points.

Figure 2. Information on the TFPM method. (a) Topographic features. (b) Grid strategy used in gap boundary detection. (c) Target patch generation. (d) Source patch generation. (e) Central point and its neighbouring points. (f) Gently sloped areas in a valley. (g) The flowchart of the topographic feature extraction method.

2.2. Patch generation

After gap boundary point detection, a gap and its neighbouring points consist of a target patch (Figure 2(c)), which is used to define the processing scope for subsequent operations, such as topographic feature interpolation and gap segmentation. The source patch (Figure 2(d)) is matching resources for target patches during the patch matching process. The details of the generation process are as follows.

Regarding the generation of the target patch, the raw point cloud is gridded into raster cells (elevation direction) with a size of 3 times the APD. For each gap, the boundary cells are clustered outwards until a rectangular region encompassing the entire gap is generated, and the points including the gap boundary points and their neighbouring points within the rectangular region with a buffer distance of 3 times the APD (d in Figure 2(c)), constitute a target patch (points within outer rectangular regions in Figure 2(c)). Source patches (Figure 2(d)) are generated according to their corresponding target patch. To address the issue of missing similar regions during patch matching due to the fixed size of the generated source patch, a solution is proposed wherein multiple source patches of varying sizes are generated for each target patch. That is, the raw point cloud is gridded at sizes equivalent to 1, 1.5, and 2 times the size of the corresponding target patch’s longer side, i.e. the size value in Figure 2(d) is multiplied times by h in Figure 2(c), and point clouds within the raster cells are considered source patches (Figure 2(d)).

2.3. Topographic feature extraction and interpolation

2.3.1. Extraction of topographic features

Commonly employed topographic feature point extraction approaches involve identifying potential feature points, followed by thinning operations to obtain the final feature points (Zhou et al. 2018). Potential feature points are usually extracted according to factors that can well describe the morphology of topographic features. Zhou et al. (2018) utilised the Signed Surface Variation (SSV, formula 2 and 3) to identify the feature points of topographic features, such as feature points in valley and ridge areas (Figure 2(a)). In terrace areas, the steepest slope variations occur along the vertical direction of the terrace riser (Figure 2(a)), and the profile curvature refers to the curvature along the direction of the maximum slope gradient and is the rate of change in the terrain slope (Minár, Evans, and Jenčo 2020). Thus, the profile curvature is suitable for terrain riser extraction. In this study, the SSV is used to identify potential feature points in valley and ridge areas, and the profile curvature is utilised to identify potential points in terrace areas.

The SSV of point $p$ is computed using formula 3, and the neighbouring points within a search radius $R$ of the centre point $p$ are used to calculate the eigenvalues of $C$ (Duan et al. 2021). (2) $\begin{array}{c}C={\displaystyle \frac{1}{n}\begin{array}{ccc}{\left[\begin{array}{ccc}{p}_1& -& \overline{p}\\ \cdots & \cdots & \cdots \\ p_k& -& \overline{p}\end{array}\right]}^T& \cdot & \left[\begin{array}{ccc}{p}_1& -& \overline{p}\\ \cdots & \cdots & \cdots \\ p_k& -& \overline{p}\end{array}\right]\end{array}}\end{array}$(2) where $n$ is the number of points within the $R$ of $p$ and $\overline{p}=\sum _{i=1}^n{\displaystyle \frac{p_i}{n}}$ is the centroid of $p$. The SSV of point $p$ is given by (Zhou et al. 2018). (3) $\begin{array}{c}SSV\left(p\right)=\{\begin{array}{cc}{\displaystyle \frac{{\lambda }_0}{{\lambda }_0+{\lambda }_1+{\lambda }_2},}& n_p\ast p\overline{p}\ge 0\\ -{\displaystyle \frac{{\lambda }_0}{{\lambda }_0+{\lambda }_1+{\lambda }_2},}& n_p\ast p\overline{p}<0\end{array}\end{array}$(3) where ${\lambda }_i$ $\left(i=0,1,2\right)$ are the eigenvalues of $C$, with ${\lambda }_0<{\lambda }_1<{\lambda }_2$, and $n_p$ is the normal vector of $p$.

The quadratic function (formula 4) is fitted using point $p$ and its neighbouring points within the search radius $R$. The coefficients $a$, $b$, $c$, $d$, $e$, and $f$ are calculated using the least squares method; details are provided in (Chen et al. 2023b). (4) $\begin{array}{c}f\left(x,y\right)=a{x}_j^2+b{y}_j^2+c{x}_j{y}_j+d{x}_j+e{y}_j+f\end{array}$(4) Then, the first and second derivatives of equation 4 are: (5) $\begin{array}{c}p={\displaystyle \frac{\delta f\left(x,y\right)}{\delta x},q={\displaystyle \frac{\delta f\left(x,y\right)}{\delta y},r={\displaystyle \frac{{\mathrm{\partial }}^{2}f\left(x,y\right)}{\mathrm{\partial }{x}^2},s={\displaystyle \frac{{\mathrm{\partial }}^{2}f\left(x,y\right)}{\delta x\delta y},t={\displaystyle \frac{{\mathrm{\partial }}^{2}f\left(x,y\right)}{\mathrm{\partial }{y}^2}}}}}}\end{array}$(5) The profile curvature and slope of point $p$ are calculated as follows. They are calculated based on formula 5. (6) $\begin{array}{c}{C}_p=-{\displaystyle \frac{p^{2}r+2pqs+q^{2}r}{\left(p^2+q^2\right){\left(1+p^2+q^2\right)}^{{\displaystyle \frac{3}{2}}}}}\end{array}$(6) (7) $\begin{array}{c}Slop{e}_p={\displaystyle \frac{180}{\pi }arctan\sqrt{p^2+q^2}}\end{array}$(7) The potential topographic feature points $P_f$ are extracted using formula 8 (Zhou et al. 2018), which extracts potential feature points on a local scale with $Factor$ values larger than the Local Mean $Factor$ (LMF), and the Global Mean $Factor$ (GMF) is applied to filter out the extracted points that have smaller $Factor$ values in the global range. Parameters α and β are used to fine-tune the quantitative relationships. Specifically, for a given point p, its LMF is the average $Factor$ value for p and its nearest k neighbouring points (k is set to 50 in accordance with Zhou et al. (2018) based on detailed trail in section 4.2), and the GMF is the average $\mathrm{F}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}$ value for all point clouds. Point p is considered a potential feature point if its $\mathrm{F}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{o}\mathrm{r}$ value exceeds both α times the LMF and β times the GMF. In Zhou et al. (2018), the SSV calculated using points within a search radius $\mathrm{R}$ of a central point is used to extract potential topographic points, and α and β are used to refine the extraction process. However, these key parameters are fixed without detailed explanations or analysis, which limits the generality of the method and may yield poor results when applied in other areas. In this work, a multiscale search radius method that computes robust geometric features on point clouds to retrieve the optimal neighbourhood size for each point (Demantké et al. 2012) is used for search radius $\mathrm{R}$ determination. Furthermore, α and β are set to 1, which will result in missing potential topographic feature points that are mainly distributed in areas with relatively flat terrains, such as river channels in valleys (Figure 2(f)). These regions normally exhibit low slope values. Thus, two steps are adapted to solve this problem. The details are shown in Figure 2(g). First, the extracted feature points (formula 8) with slope values (formula 7) lower than the average slope value of the extracted feature points (slope 1 in Figure 2(g)) are defined as seed points to expand neighbouring points of each seed point within its search radius R (obtained above). Second, the expanded points that have larger slope values than the average slope of the seed points (slope 2 in Figure 2(g)) are filtered. Finally, the retained points are combined with the extracted feature points (formula 8) as potential topographic feature points. In addition, a noise removal operation is added to reduce noise points and enhance the quality of the extracted potential points, and this approach is based on the density-based spatial clustering method, which allows the grouping of closely located points (Ester et al. 1996). Noise removal is conducted by eliminating groups with fewer points than the average number, which are considered noise points, to preserve the main topographic features.

The second step is thinning the extracted potential points. An improved HC-Laplacian method (Taubin 1995) (formula 10) is adopted, which uses $Factor$ values as weights (formula 9) to evenly distribute the thinned points among the topographic features. (8) $\begin{array}{c}{P}_f=\{p|p\in P,|Factor\left(p\right)|>\alpha \cdot \sum _{i=1}^k{\displaystyle \frac{\left|Factor\left(p_i\right)\right|}{k}\cap \left|Factor\left(p\right)\right|>\text{β}\cdot \overline{\left|Factor\right|}\}}\end{array}$(8) (9) $\begin{array}{c}w\left(p_i\right)={\displaystyle \frac{Factor\left(p_i\right)}{{\sum }_{j=1}^{n}Factor\left(p_j\right)}}\end{array}$(9) (10) $\begin{array}{c}q={\displaystyle \frac{1}{{\sum }_{i=1}^{n}w\left(p_i\right)}\sum _{i=1}^{k}w\left(p_i\right)\left(p_i\right)}\end{array}$(10) where $Factor$ is determined based on the terrain characteristics; i.e. the SSV is for hillslope areas with valleys or ridges, and the profile curvature is for terrace areas. $p_i$ denotes the neighbouring points of point $p$, as determined by the adaptive search radius method mentioned above, and $Factor\left(p_i\right)$ denotes their factor values. $\overline{\left|Factor\right|}$ is the averaged $Factor$ of the raw point clouds. $\alpha$ and $\text{β}$ are 1 in this study. $n$ is the number of neighbouring points for point $p$ within the search radius of 3 times the APD, $w\left(p_i\right)$ is the calculated weight for the neighbouring points, and $q$ is the final position for point $p$.

2.3.2. Interpolation of topographic features within gap areas

Prior to point cloud restoration, topographic feature lines, i.e. the feature lines that are supposed to pass over gap areas, are recovered to facilitate the patch matching process, which also improves the capability of the TFPM method in restoring large gaps. Topographic feature points within gaps are interpolated as follows.

Topographic feature points within a target patch are first grouped into multiple point sets. Second, for each target patch, the algorithm checks the lines ($y=kx+b$) fitted using the X and Y coordinates of the grouped point sets, such as line 1 and line 2 in Figure 3(b), and if two of them intersect within the gap (Figure 3(b)) or are approximately parallel with a distance (d in Figure 3(c)) less than 3 times the APD, then they are matched and clustered to form a pair of candidate points. Moreover, if gaps are located at terrain surfaces with complex terrain features (Figure 3(d)) or noise points that are commonly caused by slightly undulating topography (Figure 3(e)), a feature point set will be matched twice or more with other feature point sets. For example, in Figure 3(d, e), line 2 is matched with line 1 and line 3. In these situations, two point sets of these matched point sets with similar elevation values (line 1 and line 2 in Figure 3(d, e)) are selected as a pair of candidate points because feature points belonging to the same topographic feature normally have similar elevation values compared to others. Third, 2D curves ($y=a{x}^2+b\mathrm{x}+c$) are fitted using the X and Y coordinates of each pair of clustered candidate points (Figure 3(f(2))), and the X and Y coordinates of interpolated feature points are obtained by sampling the fitted 2D curves at intervals of the APD within gaps (Figure 3(f(3))). Moreover, the X and Z coordinates of each pair of candidate points are fitted to 2D curves ($z=a{x}^2+b\mathrm{x}+c$, sharing the same X coordinates), yielding the corresponding Z coordinates, which are used to obtain the final interpolated topographic feature points (Figure 3(f)).

Figure 3. Topographic feature point interpolation and patch matching processes. (a) A similar region for a gap. (b), (c), (d) and (e) The rules to determine if two point sets match. (f) Flowchart for the topographic feature recovery process.

2.3.3. Gap and target patch segmentation based on topographic features

After topographic feature point interpolation, the target patch is divided into smaller patches using the fitted curves discussed in section 2.3.2. Simultaneously, the large gap is segmented into smaller gaps. As Figure 3(f(4)) shows, the target patch and its gap inside are segmented into two parts by the fitted curve. Notably, for each segmented target patch, corresponding topographic feature points are also employed as boundaries during matching processes.

2.4. Gap restoration using the patch matching method

After topographic features are incorporated into the generated target patch, the patch matching method can be utilised to explore nonlocal similarity, which considers terrain surface similarities. As illustrated in Figure 3(a), a nongap ridge area is utilised to restore the gap within the ridge area. This approach encompasses several key operations, including patch generation, patch matching, similarity computation, and gap restoration (Fu, Hu, and Guo 2018).

2.4.1. Patch matching

When performing patch generation and segmentation, the segmented target patches with small gaps are matched with source patches that are generated with multiple sizes (section 2.2). As Figure 1(f) illustrates, a target patch is divided into two segments, and each segment is used to search for the best-matched source patch that exhibits the highest similarity (quantified in section 2.3.2). Furthermore, during the matching process, the geometric structure in a target patch and its source patches is registered to strengthen the similarity between them for better restoration, which includes the translation obtained by the difference in location between the target and source patches and the rotation obtained by the simplified iterative closest points (ICP) method (Chetverikov et al. 2002; Hu, Fu, and Guo 2019). During the rotation implemented by the ICP method, gap boundary points and their surrounding points within 3 times the APD are used as templates for the rotation of source patches.

2.4.2. Gap restoration

After the patch matching process mentioned above is applied, the source patch that exhibits the highest similarity with the corresponding target patch is utilised for gap restoration. The similarity is quantified based on the average elevation difference (formula 11). For each segmented target patch, the source patch exhibiting the lowest average elevation difference is chosen as the template for gap restoration. Specifically, each point in the found source patch and is identified and retained it if it falls within the gap area. (11) $\begin{array}{c}AED={\displaystyle \frac{{\sum }_{i=1}^n\left(d_i\right)}{n}}\end{array}$(11) where $n$ is the number of gap boundary points of a target patch and $d_i$ is the elevation difference from a gap boundary point $p_i$ to its closest point in the source patches.

2.4.3. Smoothing of the restored results

Additionally, the identified source patch is usually not completely matched with the corresponding target patch due to the inherent difference in terrain surfaces, thus, a smoothing technique is applied to enhance the overall continuity of the restored surfaces. Specifically, the restored points are adjusted according to quadratic surfaces fitted using the gap boundary points. The underlying principle is that more smoothing is applied as the restored points approach gap boundaries (formula 12). (12) $\begin{array}{c}{Z}_{adjusted}=Z_{original}+\left(Z_{surface}-Z_{original}\right)\ast e^{-d}\end{array}$(12) where $d$ is the distance from a restored point $p$ to its nearest gap boundary point, $Z_{original}$ is the original elevation, and $Z_{surface}$ is its corresponding elevation on the fitted quadratic surface.

2.5. Evaluation methods

The TFPM method was compared with the IDW (Shepard 1968), second-order polynomial method, which is based on quadratic surface fitting (QSF), and a patch matching (PM) method with the same implementation as the TFPM method but without the aid of topographic features. Gap boundary points and their neighbouring points within a search radius of 3 times the APD are utilised to implement these methods, as in the TFPM method. A comparison is performed from three perspectives: (1) visual inspection, (2) a performance assessment of the restored results in terrain modelling, specifically for the generation of grid-DEMs and TINs, and (3) an assessment of the capability to effectively use the restored results in terrain analysis, including terrain derivative extraction and re-extraction of topographic features. Specifically, first, grid-DEMs are generated by the linear interpolation method integrated by the open source point cloud processing software CloudCompare (https://www.danielgm.net/cc/) using the restored results of different restoration methods within the target patch extent; then, the Root Mean Square Error (RMSE) (Reuter, Nelson, and Jarvis 2007) is utilised to determine the accuracy of the generated grid-DEMs (formula 13). To comprehensively analyse the modelling error distributions, the frequencies of errors in the generated grid-DEMs are considered. Second, the average elevation difference (formula 14) is used to quantify point elevation distances. This metric reflects the accuracy of the generated TINs because they are constructed directly from points. Third, the average slope (formula 7) difference is calculated to measure the slope distance (formula 14). Finally, topographic features are re-extracted to evaluate the quality of the restored results of different methods in topographic feature restoration. (13) $\begin{array}{c}RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{\left({Z^{\mathrm{\prime }}}_i-Z_i\right)}^2}}\end{array}$(13) where $Z_i^{\mathrm{\prime }}$ and $Z_i$ represent the elevation value of the gird-DEM constructed using the restored points and the reference points within the target patch at the position of pixel $i$, respectively, and n is the total number of pixels. (14) $\begin{array}{c}H={\displaystyle \frac{{\sum }_{i=1}^n\left(d_i\right)}{n}}\end{array}$(14) where $n$ is the number of restored points, and $d_i$ represents the elevation or slope difference between the restored points and the nearest reference points.

2.6. Test areas and data

Four test areas, namely, Area 1, Area 2, Area 3 and Area 4 (Figure 4(a, d, g, j)), were sampled from the Loess Plateau, where suffers from severe soil erosion (Dai et al. 2023b; Xiong et al. 2023). Point clouds in the above four test areas were generated using photogrammetric methods. Original photogrammetric imagery was generated using unmanned aerial vehicles. Area 1 and Area 3 are terrace areas with APDs of 0.31 and 0.30 m respectively, and Area 2 and Area 4 are hillslope areas with APDs of 0.44 and 0.33 m, respectively. To quantitatively validate the restored results obtained with multiple restoration methods, 5 gaps were randomly and manually dug in both Area 1 and Area 2. Specifically, in Area 1, gap 1 was in a flat region, and gaps 2, 3, 4, and 5 were dug in terrace risers. In Area 2, gaps 1 and 2 were dug in valley areas, gaps 3 and 4 were in ridge areas, and gap 5 was in a gentle slope area. The dug points were used as references for accuracy validation. In addition, two test areas, Area 3 and Area 4, with real terrain gaps (Figure 4(i, l)) caused by filtering out vegetation, were used to assess the performance of the TFPM method in gap restoration in real scenarios.

Figure 4. Information for the test areas. (a), (d), (g) and (j) are point clouds. (b) and (c) are point clouds with gaps created by digging. (c) and (f) are rendered maps. (h) and (k) are partial 3D visualisations of Area 3 and Area 4, respectively. (i) and (l) are real terrain gaps caused by filtering out vegetation.

3. Results

3.1. Restored results in the test areas

During the gap restoration processes, in addition to the final restored point clouds, various intermediate results are produced, mainly including gap boundary points, as well as extracted and interpolated topographic feature points.

3.1.1. Detected gap boundary points and interpolated points within the gaps

The partially enlarged views in Figure 5(a, b) clearly demonstrate that gap areas are accurately detected, which ensures effective target patch generation and gap restoration. Figure 5(c, d) shows the extracted topographic feature points and their interpolated topographic points within the gap areas. In the terrace area, the extracted topographic feature points are well distributed across the terrace risers. Topographic feature points within gaps appear to be effectively interpolated, displaying a high degree of continuity with the surrounding terrace riser points. The partially enlarged views of the hillslope indicate that the valleys and ridges are effectively reconstructed within gap areas. These results ensure the quality of the segmented gaps.

Figure 5. Gap boundary points and topographic feature points in Area 1 and Area 2. (a) and (b) Detected gap boundary points. (c) and (d) are the extracted and interpolated topographic feature points.

3.1.2. Restored results of the different methods

To better visualise and compare the quality of the gap restoration results, the restored points are converted into TINs. From Figures 6 and 7, it is evident that, first, the restored outcomes using the TFPM and PM methods exhibit a close resemblance to actual surfaces with topographic relief. Conversely, the restored results obtained with the QSF and IDW methods tend to appear smoother compared to real surfaces, such as in Area 1-gap 3 and Area 2-gap 5. Second, the PM method has the potential to identify similar regions to gaps, but with a certain degree of uncertainty. This is the reason why topographic features are considered and incorporated into the patch matching process. For instance, the PM method achieves better results in Area 1-gap 4 and Area 2-gap 3 but yield poor results in Area 1-gap 2 and Area 2-gap 4. Third, the restored results of the TFPM method are closer to the reference data due to the restoration of topographic features. Specifically, the TFPM method excels in effectively restoring essential skeletal information, i.e. terrace risers, ridges and valleys, while the other methods fail to achieve this result.

Figure 6. Gap restoration results in Area 1.

Figure 7. Gap restoration results in Area 2.

3.2. Terrain modelling using the restored point clouds

The restored results within the target patch extents of various methods were employed to generate grid-DEMs with a resolution of 1 m and TINs to showcase their effectiveness in terrain modelling. The accuracy validation outcomes are as follows.

3.2.1. Accuracy of the generated grid-DEMs

Figure 8(a) shows the calculated RMSE values in the terrace and hillslope areas. It is evident that, first, all methods exhibit similar accuracy when dealing with gaps in relatively gently sloped regions (Area1-gap 1 and Area 2-gap 5). Second, the TFPM method consistently demonstrates superior performance in cases with multiple gaps, as evidenced by consistently lower RMSE values. Although the PM method exhibits slightly better performance for Area 2-gap 3, it performs the worst for Area 2-gap 4, possibly because the method fails to identify a similar region for the gap. The conclusions demonstrate the superior capability of the TFPM method in grid-DEM generation, with higher accuracy than the other methods. Furthermore, Figure 8(b, c) shows that the modelling errors of the TFPM method are predominantly and consistently distributed within low error intervals in the terrace and hillslope areas, which highlights the superior performance of the TFPM method in terrain modelling.

Figure 8. Accuracy verifications of the generated grid-DEMs in Area 1 and Area 2. (a) RMSE values of the generated grid-DEMs. (b) and (c) show the modelling error frequencies of the grid-DEMs.

3.2.2. Accuracy of the generated TINs

Figure 9 shows the average point elevation difference from the restored points to the reference data. Obviously, the TFPM method consistently achieves high accuracy and maintains stable performance, which indicates that the restored points generated by the TFPM method closely resemble their corresponding real surface counterparts, highlighting the superior performance of the TFPM method in TIN generation.

Figure 9. Average point elevation differences in Area 1 and Area 2.

3.3. Terrain analysis using the restored point clouds

The restored results are utilised in the extraction of slope and topographic features to assess the ability of the TFPM method in terrain analysis.

3.3.1. Terrain derivative extraction using the restored results

Figure 10(a, b) shows that the TFPM method exhibits the lowest average slope difference with the reference data, which indicates that the slope values of the points restored with the TFPM are closer to the true values than are those obtained with other methods and suggests that the spatial positions of and relationships among restored points and the surrounding points are also similar to those for the original surface points.

Figure 10. Average slope differences in Area 1 and Area 2.

3.3.2. Topographic feature extraction using the restored results

Figure 11 shows the topographic features re-extracted using the restored results obtained with different methods. Figure 11(a, b) shows that the TFPM method effectively re-extracts topographic features of terrace risers, valleys, and ridges, and the results resemble the reference data extracted using point clouds without any gaps. Conversely, the other methods nearly fail to extract points associated with terrace risers and ridges due to their poor restoration of these crucial topographic features. Partial visualisations in Figure 11(a) clearly demonstrate that the TFPM method consistently yields terrace riser points that closely match the ground truth data, whereas other methods fail to accurately recover terrace risers. In Figure 11(c), the ridge lines generated from the restored points of these methods show that the TFPM method outperforms the other methods in ridge restoration. These findings demonstrate the superior performance of the TFPM method in restoring topographic features.

Figure 11. Re-extracted topographic features in Area 1 and Area 2 using the restored results obtained with different methods.

3.4. Gap restoration in real scenarios

Figure 12 illustrates gradient terrain surfaces in real-world scenarios for Area 3 and Area 4. It is apparent that the restoration results obtained with the TFPM are more consistent with the actual terrain surfaces, and topographic features, namely, ridges, valleys, and terrace risers, are effectively and stably restored. In contrast, the other methods struggle to perform well in complex and realistic terrain cases, particularly in terms of preserving terrain features. Figure 12(b) shows that the restored valley line of the TFPM method is consistent with the nongap areas (reference data), while the other methods yield poor results. These results demonstrate the superior performance of the TFPM method in the restoration in real terrain gaps, which underscores the importance of considering topographic features and terrain similarities when performing terrain restoration in practical applications.

Figure 12. (a), (c) and (d) Gradient maps of restored terrain surfaces and their partial visualisations in Area 3 and Area 4. (b) Valley lines in Area 4 generated by the restored points of these methods.

4. Discussion

4.1. Functionality of topographic features in point cloud restoration

Topographic features have been proven to play key roles in terrain modelling (Jiang et al. 2023; Zhang, Yu, and Zhu 2022). In this study, they are incorporated into point cloud restoration. From the comparison of results, the TFPM method achieves superior accuracy in terrain modelling and analysis. To further investigate the functionality of topographic features in point cloud restoration, the QSF and IDW methods are also refined by incorporating topographic features. The interpolated topographic feature points are initially inserted into gaps, after which conventional methods are employed for gap restoration. The TFPM, TFQSF, and TFIDW methods are improved versions of the PM, QSF and IDW methods, respectively. Figure 13(a, b) shows that the improved methods can better repair gaps and recover ridge information. Furthermore, the improved methods generally yield significant advancements in the accuracy of terrain modelling and analysis (Figure 13(c, d)). These findings illustrate the crucial and positive impact of topographic features on point cloud restoration.

Figure 13. Comparisons of the restoration methods with and without topographic features based on visualisation of results, elevation distance and slope distance.

4.2. Key parameter setting in the TFPM method

The TFPM method includes three main steps, namely, gap boundary point detection, topographic feature extraction and patch matching, inevitably introducing some parameters, which increases the uncertainty of the TFPM method. Some improvements were made to enhance each step, and some parameters were optimised or removed, such as the search radius, $\alpha$ and $\text{β}$ during topographic feature extraction, and patch size during patch generation. Moreover, the details for the determination of other key parameters are as follows.

For gap boundary point detection, a search radius $R_{gap}$ from 4 times the APD to 6 times the APD was tested in Area 1 and Area 2 (Figure 14(a, d)). When $R_{gap}$ was set to 5 times the APD, some noise points appeared around gaps in both Area 1 and Area 2 (view 1 and 4). When $R_{gap}$ was set to 5 times the APD, ideal gap boundary points were obtained in Area1 and few noise points appeared in Area 2 (view 2 and 5). When $R_{gap}$ was set to 6 times the APD, ideal results were obtained in both Area 1 and Area 2 without noise points, obtaining similar results compared to the results of those in the case for 5 times APD (view 3 and 6). Thus, 6 times the APD is best threshold for $R_{gap}$ in both terrace and hillslope areas. In addition, Figure 14(b, e) shows that when $R_{gap}$ is set to be 6 times the APD, ideal results were also obtained in Area 1 and Area 2 with multiple APDs, as well as real terrain gaps in Area 3 and Area 4 (Figure 14 (c, f)). In conclusion, 6 times the APD is an ideal search radius for $R_{gap}$. Moreover, during topographic feature extraction, the parameter k used to define the local nearest k neighbouring points directly affects the extracted potential feature points. Figure 15 shows that there is no significant difference in the extracted feature points when k is set to 40, 50 or 60. Thus, k is set as 50 in accordance with Zhou et al. (2018).

Figure 14. Detected gap boundary points obtained by multiple search radii, APDs, and test areas.

Figure 15. Extracted potential feature points under multiple k values.

4.3. Topographic knowledge applications in point cloud processing

Point cloud data are widely used in various domains, such as 3D modelling (Wang, Peethambaran, and Chen 2018), feature extraction (Li et al. 2018), and Earth surface research, including terrain modelling (Chen et al. 2023b), hydrological analysis (Lyu et al. 2021), and change detection (Nourbakhshbeidokhti et al. 2019), etc. During terrain point cloud processing, topographic knowledge, such as terrain derivatives and topographic features, should be fully considered to avoid the loss of terrain information. Previous studies have attempted to incorporate topographic knowledge during terrain data processing, such as grid-DEM super-resolution (Jiang et al. 2023) and point cloud simplification (Chen et al. 2023b). In this study, topographic features are considered for point cloud restoration, which improves the accuracy of terrain modelling and analysis, as well as the restoration of topographic features. Consequently, corresponding knowledge systems for various domains should be recognised and mined for specific knowledge to control point cloud processing. For instance, road boundaries and house outlines are types of knowledge associated with point clouds for urban infrastructure, and thus deserve significant attention in 3D building modelling.

5. Conclusions

In this paper, a gap restoration method that integrates topographic features and patch matching is proposed. The TFPM method is mainly divided into three parts: (1) an improved gap boundary detection method, (2) topographic feature point extraction and interpolation, and (3) patch matching processes that find and use similar regions as templates to repair gaps.

The comparison of the TFPM method with the other methods demonstrates the TFPM method’s superior performance in terrain modelling and analysis. The improved gap boundary and topographic feature detection methods help the TFPM method accurately extract gap boundary points and topographic features. From the final restored results, the visualisation comparison shows that the TFPM method’s restored results are closer to the reference data (real surfaces); the method also commendably restores topographic features. In terms of terrain modelling, the TFPM method achieves superior performance, with lower RMSE values, good error frequency distributions, lower elevation differences relative to real surface points, and stable performance. In addition, the TFPM method’s restored results are verified as accurate in terrain analysis. Specifically, the results display smaller slope differences as the real terrain and accurately recreate topographic features.

The results demonstrate that the consideration of topographic features can help algorithms repair gaps, and stable performance is achieved, especially in areas with topographic features. Additionally, the patch matching method is utilised to consider terrain surface similarities, making the restored points closely resemble natural terrain surfaces. These two concepts provide a broad application range for this restoration strategy. Overall, this integration strategy can be easily extended to other fields for data restoration. Although the proposed method displays better results than other methods in the context of terrain modelling and analysis, it also suffers from some limitations. First, the final restoration results rely on the outcomes of preceding steps, such as the extraction of topographic features. More attention should be given to establishing effective and high-performance methods for topographic feature extraction. Second, the proposed method is not as efficient as interpolation methods because of the topographic extraction step and the consideration of global similarities rather than local information. In further studies, parallel computing, which is commonly used for high-density point cloud data processing (Lyu et al. 2021), can be applied to speed up these processes. Additionally, previous research has demonstrated the effectiveness of deep learning techniques in restoring terrain data, specifically in grid-DEM restoration, by integrating topographic knowledge (Li et al. 2022b). Consequently, considerable emphasis should be placed on designing a network that effectively incorporates topographic knowledge for the restoration of terrain point clouds.

Acknowledgments

The authors would take to thank editors and two anonymous reviewers for the useful comments on the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The test data that supports this work is available in ‘figshare’ repository with the private link ‘https://figshare.com/s/d1e7ccb5d56793ed5e89’.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China under Grant [41971333, 41930102, 42371407]; Priority Academic Programme Development of Jiangsu Higher Education Institutions under Grant [164320H116]; The Priority Academic Program Development of Jiangsu Higher Education Institutions and the Deep-time Digital Earth (DDE) Big Science Program.