LUNAR TERRAIN RECONSTRUCTION FROM MULTI-VIEW LROC NAC IMAGES BASED ON SEMI-GLOBAL MATCHING IN OBJECT SPACE

Most of the lunar surface area has been observed from different viewing conditions thanks to the on-orbit work of lunar orbiters, a large amount of images are available for photogrammetric three-dimensional mapping, which is an important issue for lunar exploration. Theoretically, multi-view images contain more information than a single stereo pair and can get better 3D mapping results. In this paper, the semi-global matching method is applied to the object space, and the steps of cost calculation, cost aggregation, and elevation calculation are performed to obtain the three-dimensional coordinates directly. Compared with the traditional image-based semi-global matching method, the object-based semi-global method is more easily extended to multi-view images, which is beneficial for applying multi-view image information. In addition, it does not require steps such as stereo rectification and forward intersection, that is, the overall pipeline is more elegant. Using the LRO NAC images covering Apollo 11 landing area as the experimental data, the result shows that the object-based semi-global matching is competent for the multi-view image matching and the multi-view image result achieves higher accuracy and more details than the single stereo pair. Furthermore, the experimental results of Zhinyu crater data show that this method can also alleviate the uncertainty of the lunar orbiter's positioning to some extent.


INTRODUCTION
Topographic reconstruction of the lunar surface is critical to engineering applications and scientific research (Karachevtseva et al., 2013;Wu et al., 2014). Generally, there are two types of data sources including laser altimetry data and stereo image data for lunar topographic reconstruction. Photogrammetrically derived DEM from stereo images can achieve better spatial resolution relative to laser altimetry data.
In the past two decades, the United States, China, Japan, India, and other countries have carried out many lunar exploration missions and obtained a large amount of image data that can be used for lunar terrain reconstruction. With continuous data acquisition from lunar orbiters, many areas of the lunar surface have been observed repeatedly under different viewing conditions. Among the recent orbital data, the highest-resolution lunar orbital imagery is achieved by the Lunar Reconnaissance Orbiter Camera (LROC) Narrow-Angle Camera (NAC) (Robinson et al., 2010). Notably, the LROC NAC images have covered nearly the entire lunar surface with a resolution of 0.5-2.0 m.
Compared with the earth, the lunar surface changes little at long durations. Thus, LROC NAC images obtained at different times and from different views can be used for multi-view reconstruction; it also poses challenges to use those data to generate high-quality DEMs.
The photogrammetric processing procedures of lunar orbital images mainly include camera geometric model construction, block adjustment, dense matching, 3D point generation through the forward intersection, DEM interpolation, etc. However, there * Corresponding author exist many limitations on the lunar surface, such as the illumination changes, geopositioning uncertainty and so on.
Semi-global matching (SGM) is a powerful stereo matching method which can achieve a good compromise between accuracy and efficiency (Hirschmuller, 2007), the matching problem caused by geopositioning uncertainty conditions can be solved by transferring semi-global matching from image space to object space. Moreover, semi-global matching in object space can be extended to multi-view images. In this paper, we use a semiglobal optimization strategy in object space for multi-view LROC NAC image matching to generate high-resolution DEM.

Matching in Image Space and Matching in Object Space
A large amount of in-house photogrammetric software was developed by various research institutions at home and abroad. The Intelligent Robotics Group (IRG) at the NASA Ames Research Center developed an open-source image matching program Ames Stereo Pipeline (ASP). In ASP, Normalized Cross Correlation (NCC) was regarded as the metric to calculate the matching cost, and WTA (winner takes all) strategy is used to find the disparity corresponding to the minimum matching cost as the final disparity. Finally, the 3D point cloud is obtained by space intersection. A DEM with 40m / pixel resolution was produced with 71 pairs of stereoscopic images obtained by APOLLO 15 selected (Broxton et al., 2009). The research team from the University of Parma in Italy developed a software called Dense Matcher to generate DEM. This software implements the NCC matching method, least-squares matching method, and other matching methods involving geometric constraints such as Multiphoto Geometrically Constrained Matching (MGCM). Using HiRISE and LROC NAC images as the experimental data, the results show that the accuracies of the DEM generated by the least-squares matching method developed in DM and the DEM generated by ASP are similar (Re et al., 2012a(Re et al., , 2012b. In 2019, the research team of the Hong Kong Polytechnic University developed a photogrammetric software named PLANETARY3D software for processing planetary images, which uses semiglobal matching to find image correspondences. Experimental results using three planetary data sets, i.e., MRO-CTX, MRO-HiRISE, and Chang'e-2, show that the software can generate high-quality planetary DEM (Hu et al., 2019).
Because there is no atmosphere on the lunar surface, the image matching results are greatly affected by the change of illumination conditions. Besides, the texture information of most regions is weak. Therefore, matching strategies such as epipolar constraint, coarse to fine strategy, and triangular network constraint are often adopted in stereo matching to improve the success rate and accuracy of matching. Wu et al. (2014) integrated 7m-resolution Chang'E-2 imagery and Lunar Orbiter Laser Altimeter data to generate a digital elevation model (DEM) with a resolution of 20 m for the entire Sinus Iridum landing area and thus provided data support for the selection of Chang 'E-3 landing site. In the process of terrain reconstruction, the search range is constrained by using the triangle network. Li et al. (2018) used the area-based image matching method to find Chang'e-2 image correspondences and finally produced a full-month DEM.
As described abov，the methods used for lunar image matching are generally in image space. The applications of matching methods in object space are mainly for earth remote sensing images, but they should be applicable in planetary 3D reconstruction. Ji et al. (2009Ji et al. ( , 2011 extended the traditional stereo vertical line locus method (VLL) to multiple images, proposed an algorithm called MVLL (multi-image VLL) which is suitable for multi-view image matching. Zhang et al. (2007) used the average value of the multi-image pair window correlation coefficient as a metric, and indirectly used epipolar constraints to perform multi-view matching in object space, effectively integrating the information of different stereo image pairs. Yuan et al. (2009) proposed a method for multi-image matching by integrating image and space information. Matching methods in object space usually uses a matching window of a certain size to find image correspondences. There is no global constraint in the matching process, which will influence the reliability of matching to some extent.

Global Methods and Local Methods
Stereo matching methods can be subdivided into two categories, global optimization algorithms, and local window algorithms. Some researchers use the local matching method to find lunar orbital image correspondences, e.g. the NCC and the LSM (Least Squares Matching) (Re et al., 2012b). However, local methods lack global constraints, which will result in poor matching quality in occlusion and insufficient texture area.
Global methods often can achieve a better result than local methods due to global constraints. Typical global methods include matching methods using Markov model, semi-global matching, etc. Peng et al (2014) proposed an adaptive Markov random Field (aMRF) model for the dense matching, real rover images from the Mars Exploration Rover mission and Chang'e lunar orbiter images experiments demonstrate the effectiveness of the method, but the method is computationally expensive. Semi-global matching applies two parameters P1 and P2 to impose smoothing constraints on different disparity changes and uses multi-path one-dimensional constraints to approximate the two-dimensional constraints, which greatly improves efficiency while ensuring accuracy. Semi-global matching has been widely used in many software and algorithms to generate Lunar DEM. For example, the algorithm used in the dense matching module of PLANETARY3D is based on the semi-global matching method. Besides, ASP also introduced semi-global matching in its recent versions (Hirschmuller et al., 2007;Beyer et al., 2018). Bethmann et al. (2015) applied semi-global matching in object space, which essentially simplified the multi-view matching process and used indoor and aerial images for experiments. However, they didn't perform quantitative analysis on aerial image reconstruction results.
As far as lunar images are concerned, there is little research on integrating semi-global matching with stereo matching in object space. when there exist geopositioning uncertainty conditions such as a small overlapping condition, the in-house software may not be available. Additionally, it's hard to find lunar image correspondences due to insufficient texture and other difficulties. However, those problems can be alleviated to some extent by applying semi-global matching in object space with constraints such as height search range constraint imposed by initial DEM and a global constraint imposed by semi-global matching optimization strategy.

Overall Workflow
The overall workflow of lunar terrain reconstruction from multiview images is shown in Figure 1. The inputs for the method are multi-view lunar images and their corresponding rational polynomial coefficients (RPCs) and an initial low-resolution DEM. RPCs are used to establish the relationship between imagespace coordinates and object-space coordinates with ratios of polynomials . Having subdivided the object space into dense voxels, the coordinates of each voxel are backprojected to images to get corresponding positions in image space. By using pre-defined correlation coefficient window size and its window-side length, the correlation coefficient between these image points can be calculated, the maximum value of the results is denoted as ρ max , and (1 − ρ max ) can be assigned to the corresponding voxel, which is the so-called voxel cost in this paper. After that, semi-global optimization is used to aggregate cost in object space, and the spline function of the cost values of specific planar coordinates is fit to find the final height corresponds to the minimum cost. Finally, the final DEM is obtained after outlier rejection.

Traditional Semi-global Matching
The main idea of the traditional semi-global matching method is the approximate minimization of the global 2D energy function: where ∑ ( , D p ) p is the data term, which represents the sum of all pixelwise matching cost when the disparity map is D. The second term and the third term of the right side in Eq.1 are the 2D smoothness terms, which are used to add penalties for all pixels q in the neighborhood N p of p. T[] is a logical operator, which returns 1 when the condition is satisfied and 0 otherwise. P 1 and P 2 represent small penalty and large penalty respectively, which are used to penalize small disparity changes and large disparity changes.
However, the minimization of this energy function is NP-hard (Boykov et al., 2001). Thus, SGM aggregates costs along several different 1D paths to approximate all directions. The cost L r (p, d) along the path r is defined recursively as: The first term is the data term, which means the matching cost. The second term is the smoothness term, which equals to the lowest cost of the previous p-r of the path, including the appropriate penalty for discontinuities. The remainder of the equation subtracts the minimum path cost of the previous pixel to make L r (p, d) not too large. The so-called path can be understood as a pixel at a specific position in the neighborhood of the pixel p (x, y). Figure 2 shows the way that SGM aggregates costs from several 1D paths. In Figure 2, the number of paths is 16. Generally, the more the number of paths, the greater the time consumes, but the result is relatively better. It's necessary to choose an appropriate number of paths.

Figure 2. Aggregation of costs
Finally, the results of the cost aggregation from several paths are fused as: The final disparity will be derived from ( , ) by searching the minima of S(p, d) for each pixel p.

Cost Calculation in Object Space
When transferring image-based SGM to object-based SGM, the first step is to subdivide the object space into dense 3D voxels. And the energy function was modified correspondingly: The first term is the sum of the matching cost of each voxel. The second and the third term are the smoothness terms, used to add penalties for all voxel q in the neighborhood p of p, similar to Eq.1. T[] is a logical operator ， which returns 1 when the condition is satisfied and 0 otherwise. P 1 and P 2 represent small penalty and large penalty respectively, used to penalize small height changes and large height changes instead of disparity changes.
The purpose of calculating cost for each voxel is to get data term ( , , ).

Determining Adaptive Height Search Range:
Before calculating cost, object-based SGM should subdivide the object space into dense voxels. The size of voxel (∆ , ∆ , ∆ ) should be predefined considering mean ground sampling distance and base-to-height ratios. It's necessary to constraint the height search range because calculating cost for voxels will consume large memory and time. Lunar global DEM such as SLDEM2015 (a combined product of Lunar Reconnaissance Orbiter Laser Altimeter (LOLA) and a DEM generated by the Japanese Selenological and Engineering Explorer (SELENE) terrain camera) (Barker et al., 2016) can be used as initial DEM to constraint the height search range.
Erosion and dilatation operators can be used to find local maximum and minimum values. Accordingly, the height search range can be determined by applying erosion and dilatation operators. As shown in Figure 3, after using the initial DEM to constraint the height search range, the number of voxels needed to be calculated is greatly reduced.

Figure 3. Height search range for dense voxels
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition)

Cost Function:
There are many cost functions often used in photogrammetric processing, such as Census or NCC. Correlation coefficient is widely used as a similarity measure for lunar image matching due to its robustness. Therefore, in this work, correlation coefficient is used as the foundation of a cost function which is defined as: where and ′ represent reference window and matching window respectively， means the sum of the squares of the grey values of each pixel in the reference window ， ′ ′ is the sum of squares of the grey values of each pixel in the matching window，S ′ is the sum of the product of the grey value of the pixel at the corresponding position in the window to be matched and the corresponding position in the reference window， and ′ are the sum of the grey values of the pixels in the respective windows. The larger the correlation coefficient, the stronger the correlation between the images. A series of correlation coefficient values will be obtained by calculating the correlation coefficient between reference image windows and matching image windows. The window size and window-side length are predefined, e.g. window size can be set as 45 and window-side length can be set as 0.3 pixels, thus the length and width of the window equal 45*0.3=13.5 pixels. It helps to catch the small changes in image space caused by Z value change in object space by setting window-side length as a subpixel value. The grey values of the windows are interpolated using the bilinear interpolation method. Considering that there may be mismatches, results may be affected if the series of correlation coefficients are directly averaged. Therefore, the maximum value of the correlation coefficients between the reference image and matching images is chosen as the similarity metric for multiimages, denoted as ρ max and the matching cost of the voxel is further denoted as 1 -ρ max . So far, the matching cost of the voxels can be calculated, providing a data basis for the cost aggregation step.

Cost Aggregation in Object Space
The calculation of voxel cost aggregated from different paths in object space was modified as: The costs after aggregating from several different paths directly influence the final 3D point determination. Being similar to traditional SGM, the cost aggregation from several different paths can be fused as: The final 3D point position can be derived from Eq.7 by searching the minimum in ( , , ) for each voxel.

Height Calculation
The last step of object-based SGM for lunar terrain reconstruction leading to final DEM generation is calculating and refining the height value of each voxel. After cost aggregation, the struct ( , , ) was generated, which can be used to determine the final DEM. WTA was a choice to find the final height, it's also an easy and effective way to get the final 3D point position. However, fitting a spline function for each position ( , ) is helpful to find the subpixel Z value. Therefore, in this study, we fit the spline function for every planar position ( , ) to find its corresponding height value.

DEM Assessment
The performance of object-based SGM for lunar terrain reconstruction is assessed by evaluating the accuracy of the produced DEM. Since there is no control point on the lunar surface, it' s not possible to set several checkpoints to evaluate reconstruction precision. Also, it's hard to get a true lunar DEM as a reference. However, lunar orbiters have obtained a large number of high-resolution images in some areas, which have been used to derive high-resolution and high-precision DEM products. Those products can be used as reference DEM. The statistics of differences between produced DEM and reference DEM are applicable metrics to evaluate the produced DEM. Meanwhile, some subjective image quality evaluation metrics like Naturalness Image Quality Evaluator (NIQE) and Blind/Referenceless Image Spatial Quality Evaluator (BRISQUE) are widely used because they don't need a true reference (Mittal et al., 2012a;Mittal et al., 2012b). To evaluate the accuracy of DEM generated from multi-view images, the root-mean-square error (RMSE) is used between produced DEM and reference DEM, subjective image quality matrices NIQE and BRISQUE are also used as evaluation metrics.
On the other hand, to prove the accuracy of object-based semiglobal matching method for stereo pairs with small overlapping, the produced DEM will be compared with those produced by commercial photogrammetric software and widely used imagebased stereo matching method, respectively. Quantitative matrices may not be needed in this evaluation because visual effects are enough to show a difference.

EXPERIMENTAL EVALUATION
Two experimental datasets have been used to evaluate the proposed method in this study. The SLDEM (Barker et al., 2015) was selected to provide initial values. As for the reconstruction evaluation, for the areas where there exist high-precision and high-resolution DEM products, we use the DEM product which generated by images whose resolution is higher than experimental images as a reference. Moreover, the statistics of the differences between reference and produced DEM and image quality evaluation criteria, NIQE and BRISQUE, are used as evaluation criteria.

Evaluation Using Zhinyu Crater Data
Zhinyu Crater is a fresh crater located inside the Von Ká rmá n Crater, the landing site of Chang'e-4 probe. This evaluation aims to prove the advantages of the object-space semi-global matching method in handling small ratio conditions than commercial photogrammetry software and image-based matching method. A stereo pair of LROC NAC images involving four individual images (see Table 1) covering Zhinyu Crater are used as experimental data. These images can form three overlapping conditions, which are m184732140re and m184724991re, m184732140le and m184724991le, M184732140LE and M184724991RE. Among them, M184732140LE and M184724991RE form small overlapping stereo images, which may influence the completeness and correctness of generated DEM, thus becomes a problem in generating DEM. We regard these images as multi-view images rather than a stereo pair that require three stereo models according to the image overlap situation when we use the object-based semi-global matching method to generate DEM. And the resultant DEM is compared with the product generated from the four images by commercial photogrammetry software PCI and a DEM generated by a widely used matching method based on the triangle constraint, respectively (Peng et al., 2014). Table 1 lists the imaging conditions.  The DEM generated by a matching method based on the triangle constraint is shown in Figure 4(b). M184732140LE and M184724991RE are used to fill the gap area in figure 4(b). Figure  4(c) depicts the DEM produced by the object-space semi-global matching. Caused by orbit and attitude determination error, height inconsistency exists even when the matching result is correct, as shown in Figure 4(b). Nevertheless, applying semiglobal optimization in object space can guarantee the effectiveness of matching to some extent. Meanwhile, the height search range constrained by initial DEM helps to solve the small image overlapping ratio problem, as shown in Figure 4(c).

Evaluation Using Apollo-11 Landing Site Data
This investigation uses three LRO NAC images covering the Apollo-11 landing site as experimental data. To evaluate the advantages of multi-view images than a single stereo pair, we first choose a stereo pair with the best geopositioning precision from all the dual-image combinations, which are M1114021499R and M1126972080L. Meanwhile, the third image M1126986303R is added, which forms the best threeimage precision combination with M1114021499R and M1126972080L . Table 2 lists the detail information about these three images.
For the Apollo 11 landing site, the LROC team generated a 2m grid spacing DEM product using M150361817 L/R and M150368601 L/R, which can be used as a reference for precision evaluation due to that the pixel scales of M150361817 L/R and M150368601 L/R are about 0.5m, much higher than the images we used in the experiment. The product can be downloaded from http://wms.lroc.asu.edu/lroc/view_rdr/NAC_DTM_APOLLO11.  Figure 5 illustrates the shade relief maps of the DEMs generated by dual-image combination and three-image combination. Table  3 lists the experimental parameters, where "planar grid refinement times" equals the value that initial DEM's spatial resolution over the generated DEM's spatial resolution and thus controls the generated DEM's spatial resolution, "Correlation coefficient window (pixel)" is the window size when calculating the similarity between the image points projected by some voxel, "Window-side length (pixel)" is the size of the correlation coefficient window (the value of the window was resampled using bilinear interpolation), "The interval of Z" as the name suggests, controls how many meters the Z value changed each time, "Small penalty" and "Large penalty" are the parameters used to penalize small Z value changes and large Z value changes.
(a) (b) Figure 5. Shaded relief maps of the DEM generated from dualimage combination (a) and three-image combination (b), respectively.  Table 4 lists the RMSE of height differences, NIQE, BRISQUE and inliers of dual-image and three-image combinations. The statistical results show that the three-image combination is slightly better than the dual-image combination and thus further demonstrates the advantages of multi-view images to some extent.

CONCLUSIONS
In this study, we applied object-based semi-global matching method for lunar terrain reconstruction from multi-view LROC NAC images. The experimental result using Zhinyu Crater data demonstrates that the method could alleviate the matching problem caused by the uncertain conditions of multi-view image geometric positioning to some extent. The experimental result using Apollo-11 landing site data demonstrates the advantages of multi-view images relative to a single stereo pair.
Further research will be performed to improve the computational efficiency of our method based on GPU. Moreover, more planetary data will be used to test the effectiveness and accuracy of the method. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition)