3D POINT CLOUD COMPLETION USING TERRAIN-CONTINUOUS CONSTRAINTS AND DISTANCE-WEIGHTED INTERPOLATION FOR LUNAR TOPOGRAPHIC MAPPING

: Stereo vision has been proven to be an efficient tool for 3D reconstruction in lunar topographic mapping. However, point clouds reconstructed from pairs of stereo images always suffer from occlusions and illumination changes, especially in Lunar environments, resulting in incomplete geometric information. 3D point cloud completion is usually required for refining photogrammetric point clouds and enabling further applications. In this work, we address the problem of completing and refining 3D photogrammetric point clouds based on the assumption that 3D terrain should be continuous and with consistent slope change. We proposed a generalized strategy for 3D point cloud completion of lunar topographic mapping, including distance-weighted point cloud interpolation, terrian-continuous constrained outlier detection, and contour-based hole filling. We carried out experiments on two datasets of point clouds generated from 12 pairs and 6 pairs of stereo LROC NAC images covering the Apollo 17 and the Chang’E-4 landing sites, respectively. As a result, the holes in the initial DTM have been smoothly filled and the completeness of the whole DTM has been greatly improved. The incomplete area of the experimental areas has dropped by 100% and 93%, respectively. Finally, we constructed DTM with a resolution of 10 m covering a 33 km × 60 km area of the Apollo 17 landing site with RMSE of 4 m and a 12 km × 56 km area of Chang’E-4 landing site with RMSE of 4 m compared with LOLA laser points as a reference.


INTRODUCTION
As an essential task in lunar exploration, topographic mapping of lunar surfaces has drawn lots of attention (Wu et al., 2014, Di et al., 2012, Rosiek et al., 2001), especially in the task of geology analysis and landing site selection.Currently, different types of data acquisition techniques and data products have extensively promoted the resolution, precision, and coverage of lunar topographic mapping (Beyer et al., 2018).As the data source with the highest spatial resolution, Lunar Reconnaissance Orbiter Camera Narrow Angle Cameras (LROC NAC) images have been widely used for producing large-area digital terrain models (DTMs).For generating high-quality DTM enabling large-area coverage, pairs of stereo images should be utilized for making point clouds with high density and precision.However, due to occlusions and illumination changes when reconstructing point clouds, a loss of geometric information and incomplete point clouds always occur.To refine point clouds and meet the requirement of generating high-quality DTM for further applications, the process of postprocessing, mainly including point cloud completion is required (Huang et al., 2020).Thus, in this work, we address the problem of completing and refining point clouds generated from pairs of LROC NAC images.In this paper, we propose a generalized strategy for 3D point cloud completion, including point cloud interpolation, outlier detection, and hole filling.

RELATED WORK
To achieve large-scale topographic mapping, a common strategy involves two key steps, namely the point cloud fusion (Huang et al., 2021) and point cloud completion.* Corresponding author

Point cloud fusion
3D point cloud fusion can be categorized into the following three main types of methods, voxel-based, TIN-based, and point-based methods.Voxel-based fusion partitions the space into a voxel grid and uses different algorithms to fuse data points that are in the same voxel.The confidence of the fused points can be assessed by calculating the closest distance from the voxel to the target point as well as information such as the normal direction of the point and the angle of the inter-point vector (Curless andLevoy, 1996, Yemez andWetherilt, 2007).Although this method is sensitive to splicing errors and point cloud noise and is relatively computationally and storage intensive, it is still an effective data fusion method in some applications.Triangular mesh-based fusion methods first require the construction of a triangular mesh of the point cloud, and then the overlapping regions are determined by various methods.(Sun et al., 2003) retained the exact triangles of the overlapping regions and reconnected the remaining triangles for data fusion.However, this method does not handle the overlap of the fusion region well and the process of constructing the triangular mesh is time-consuming.Point-based fusion methods deal with overlapping points directly and perform data fusion by removing or covering the overlapping points involved in the fusion.(Zhou and Liu, 2008) proposed a method for point cloud fusion based on the K-means clustering method, which leads to a redundancy-free single-layer point cloud model.DTM fusion is an important application scenario for data fusion in remote sensing (Schmitt andZhu, 2016, Bagheri et al., 2017).First, the simple averaging method consists of calculating the mean elevation of the input DTM on a per raster basis (Leitão and De Sousa, 2018).The fusion of homologous data works well (Banu, 2011, Kaur et al., 2021).Based on this, the weighted average method was developed, which is currently the most commonly used DTM fusion method.The optimal solution does not always obtain new elevation values by averaging all the corresponding elevation points (Hoja andd'Angelo, 2009, Tian et al., 2018).In weighted averaging (WA), weights are used to quantify the influence of the input DTM at each grid or surface location and are a function of relative elevation accuracy (Schindler et al., 2011).In practice, weighted averaging has been shown to give desirable results (Tran et al., 2014, Deng et al., 2019).Variations in weights are influenced by several factors such as scene characteristics, sensor technology, and the method used to generate the original DTM (Schindler et al., 2011).However, even slight elevation differences (on the order of a few meters) between datasets can lead to rough edges at the fusion interface (Deng et al., 2019).In weighted averaging fusion, weights are usually estimated based on the error between the DTM and the reference data.However, since reference data sources (and height error maps) are not always readily available (Papasaika et al., 2011, Pham et al., 2018), this work is still challenging.

Point cloud completion
Although the problem of missing data can be nearly solved by fusing data from different viewpoints and sources, incomplete geometric information is still unavoidable due to occlusions and environmental changes.In this case, the problem of completing 3D geometry should be addressed.Recently, point cloud completion has been mainly innovated from two directions, namely geometry-based and learning-based methods.
The geometry-based method utilizes interpolation for hilling holes of 3D point clouds under the assumption that surfaces should be continuous (Berger et al., 2014, Davis et al., 2002).The geometry-based method is effective in most cases of reconstructing surfaces without complex structures.The other obvious drawback of the geometry-based method is that it cannot deal with large missing areas.The learning-based methods are widely investigated in recent years and achieved excellent performances (Yang et al., 2018, Huang et al., 2020, Wen et al., 2020, Yu et al., 2021, Wang et al., 2020).By utilizing the deep network, the complex local geometry can be learned and the latent denoising optimization method can be involved.
Although many learning-based point cloud completion methods have shown great potential, for the lunar surfaces, there lack of large-scale training data.In this paper, we focus on proposing a generalized procedure of 3D point cloud completion for precise topographic mapping.The major contributions of this paper are abstracted as follows: • We proposed a general workflow that involves point cloud interpolation, outlier detection, and hole filling for completing and refining the topographic point cloud.In the proposed workflow, terrain continuity is addressed.
• We tested the proposed method using datasets covering two test areas, including the Apollo 17 and the Chang'E-4 landing sites, and proved the effectiveness of the proposed method.

METHODOLOGY
The workflow of 3D point cloud completion consisted of three steps, including distance-weighted point cloud interpolation, terrian-continuous outlier detection, and hole filling (as shown in Fig. 1).
Figure 1.The processing workflow.

Point cloud interpolation
Point cloud interpolation is applied to transform discrete point cloud data into a continuous terrain surface model, namely DTM.In this paper, we used the distance-weighted interpolation method.The basic idea is that in the point cloud data, the closer the point to the target point, the greater the weight of the point, and the farther the distance, the smaller the weight of the point.The distance-weighted interpolation method determines the weight of each point by calculating the distance between the target point and its surrounding points, and then the attribute values of these points are weighted and averaged to obtain the interpolation result of the target point.The inverse distance interpolation method requires a series of preprocessing and parameter settings: where zp is the interpolated point elevation value, di is the distance between the point to be interpolated and the sampling point, and u is the weight index.The number of nearest neighbor points used for interpolation should be determined first, as well as the weighting formula.
Figure 2. The neighboring grid points used for calculating slopes.

Outlier detection
The second step is the detection of outliers of DTM.Although the constraints based on elevation deviation have been applied in the filtering of individual point clouds (Balta et al., 2018) and the global registration of multiple point clouds (Xu et al., 2023), there are still fewer coarse outliers in the fused DTM, which are mainly distributed near the hole area (dome, impact crater, etc.).Thus, before filling the missing area in the fused DTM, a post-processing step should be conducted to detect and reject outliers, facilitating the subsequent hole-filling process.In this paper, the outlier detection algorithm based on slope information (Hannah, 1981) is used to detect outliers.
As a basic attribute to characterize the terrain, the slope is effective to detect the coarse pixels in the DTM data through the slope information.As shown in Fig. 2, take grid point A as an example, its row and column number is (i, j), and it has 8 neighboring grid points 1 point (i − 1, j − 1), 2 points (i, j − 1), ..., 8 points (i + 1, j + 1), and 6 slope values are computed for point A and the 8 neighboring grid points in the row and column directions.Slopes of the row and column directions can be calculated as: where H(i, j) is the elevation value of the raster pixel, Dist(i− 1, i) is the spatial grid size of the DTM, and 6 slope values for each of the row and column directions can be obtained according to the equation as above.Using the slopes of the row and column directions, the differences in slope change (DSC) can be calculated: Dslopey(i, j, 1) = slopey(i, j) − slopey(i, j − 1) Dslopey(i, j, 2) = slopey(i, j) − slopey(i, j + 1) (3) By using the DSC values of all data points, if grid point A is not a coarse point, its DSC in the same direction (e.g., row direction) should be consistent, and this property is utilized to determine the coarse point by using the information of the relevant statistics of DSC.By calculating the Root Mean Square Error (RMSE) in the row and column directions for each grid point, and the sum of the DSC values for each data point in the same direction is used to calculate the RMSE value for this summation since the sum of the DSC coefficients for the same point in the same direction in the presence of a consistent Slope Change should be close to zero, and inconsistent slope changes will present a larger value for this value.
The threshold is set as K times of RMSE.The threshold value of roughness rejection is obtained by the above equation, and the grid cells whose DSC values in both row and column directions are greater than the threshold value are determined to be roughness, which is rejected, and the raster value is set to null.
Taking a text area of the Apollo 17 landing site as an example, k is set as 3, and T hresholdx=0.196135and T hresholdy=0.197602.865 pixels are retrieved as the coarse difference values in a total of 18,605,325 rasters of the fused DTM.As the proportion of coarse difference is only 0.004%, the slope grading map can be found that the coarse difference is distributed at the edge of the impact crater in this localized area, and there are more pixels with inconsistency in its slope change in the area with larger slope, and the coarse difference elevation values also show discontinuity, the above pixels are effectively eliminated and nulled.

Hole filling
After removing outliers in the fused DTM, the last step is hole filling.In this paper, the idea of hole filling is to retrieve the hole area by detecting a set of hole contour points and fill the holes by linear interpolation.The algorithm for hole detection in this paper is the findContours function (Suzuki and Be, 1985), which is used to determine the location and size of holes and repair these holes by filling the contours.In the detection of contours, all contours, including both inner and outer contours, are detected.However, instead of establishing a hierarchical relationship between contours, all consecutive contour points on the object boundary are saved into the contours vector.With the above contours retrieval rules, holes can be accurately detected.
In hole filling, contours can be used as masks to fill the area surrounded by contours using linear interpolation in both directions of ranks and columns to achieve the purpose of filling holes.Fig. 3 shows the binarization of DTM, the result of contour detection, and the result of hole hilling of a test area.It can be seen that the holes have been effectively retrieved and smoothly filled.

Experimental data
In this paper, we tested the proposed point cloud completion strategy using data obtained in two areas.The first one is the Apollo 17 landing site (19.4 • -21.3 • N, 29.9 • -31.6 • E) with an area of 50 km × 57 km.A total of 12 photogrammetric point

Experimental results
Fig. 4 shows the results of point cloud completion of the Apollo 17 landing site.From the figure, we can see that the elevation of the area is highly variable, with a difference of nearly 2800 m between the highest and lowest elevations.Both the elevation map and the rendering of the mountain shadows show that the Apollo 17 area is characterized by steep and complex topography, with impact craters of various sizes and depths, as well as numerous domed landforms.The incomplete area mainly lies on the border of the mountain.After completing and refining point clouds, all the areas with holes are filled.Evaluating the final DTM with LOLA laser points, the RMSE between the refined DTM and the checkpoint is about 4 m.
Fig. 5 shows the original DTM, and DTM after postprocessing of point cloud completion of Chang'E-4 area.From the figure, we can see that the terrain of the test area is relatively flat, and the difference between the highest and lowest elevations is less than 500 m.From the elevation map and the mountain shadow rendering map, we can see that there are some smallscale impact craters in the Chang'E-4 area, and there are also a few domed landforms above the Lunar surface.The incomplete areas mainly lie around or inside the craters.After the process of point cloud completion, the holes sparsely distributed in the experimental area are effectively retrieved and smoothly filled.
The area of holes dropped by 93%.Evaluating the final DTM with LOLA laser point, the RMSE between the refined DTM and the checkpoints is also about 4 m, which shows that the final DTM is of high quality.
To further evaluate the quality of the produced DTM, we compare the generated DTM with a resolution of 10 m with SLDEM with a resolution of 60 m.For a better illustration, we obtained cross-sections of the two DTMs.As shown in Fig. 6, in the terrain profile with a length of about 3000m, the overall degree of agreement between the two kinds of terrain data reflects the terrain undulation, which fully demonstrates the reliability and effectiveness of the terrain construction method of this paper; secondly, the DTM of this paper reflects more fully in the degree of detail, and the more minor terrain undulation can be reflected in the DTM constructed in this paper, especially in the region of Chang'E-4, due to the fact that the elevation is relatively gentle, the DTM constructed in this paper is more accurate, and it is more accurate than the DTM constructed in this paper, especially in the Chang'E-4 region.Especially in the Chang'E-4 area, due to the gentle elevation, the DTM constructed in this paper is obviously more continuous, while the SLDEM shows poorer continuity.

CONCLUSION
In this work, we proposed a 3D point cloud completion strategy based on the assumption that terrain should be continuous and with consistent slope change.The strategy involved distanceweighted point cloud interpolation, terrian-continuous constrained outlier detection, and contour-based hole filling.By completing 3D point clouds, the continuity of the Apollo 17 and Chang'E-4 landing sites was improved and the incomplete areas dropped by 100% and 93%, respectively.Finally, we achieved a DTM mapping of the experimental area with a resolution of 10 m and the RMSE between the DTM with LOLA laser point was about 4 m.
Although the proposed method provided satisfying results for the two test areas, there are still some limitations of the proposed method.First, the proposed method is limited to dealing with small missing areas.Second, there is a lack of uncertainty evaluation of the generated 3D topography, especially for the area of hole filling.In future work, learning-based methods can be involved in solving the filling of large areas and a process of uncertainty measurement should also be considered.

Figure 3 .
Figure 3.The process of hole filling.(a) is the binarization of DTM.(b) is the contour detection.(c) is the original DTM.(d) is the completed DTM.

Figure 4 .
Figure 4. Results of terrain completion of the Apollo 17 landing site.(a) and (b) are the initial DTM and the DTM after point cloud completion.(c) and (d) are corresponding hill-shaded results.clouds were generated after the stereo image pairs were selected by image filtering to obtain 79 images.The second one is the Chang'E-4 landing area (45.457 • S, 177.588 • E) with a range of 12 km × 56 km.Here, 6 stereo pairs with 72 raw images were used for generating photogrammetric point clouds covering the experimental area by 3D reconstruction of LROC NAC images and global fusion of multiple point clouds.We treated the fused photogrammetric point cloud as input for the further processing of 3D point cloud completion.

Figure 5 .
Figure 5. Results of terrain completion of Chang'E-4 landing site.(a) and (b) are the initial DTM and the DTM after point cloud completion.(c) and (d) are corresponding hill-shaded results.