An Integrated Photogrammetric and Photoclinometric Approach for Pixel-Resolution 3D Modelling of Lunar Surface

High-resolution 3D modelling of lunar surface is important for lunar scientific research and exploration missions. Photogrammetry is known for 3D mapping and modelling from a pair of stereo images based on dense image matching. However dense matching may fail in poorly textured areas and in situations when the image pair has large illumination differences. As a result, the actual achievable spatial resolution of the 3D model from photogrammetry is limited by the performance of dense image matching. On the other hand, photoclinometry (i.e., shape from shading) is characterised by its ability to recover pixel-wise surface shapes based on image intensity and imaging conditions such as illumination and viewing directions. More robust shape reconstruction through photoclinometry can be achieved by incorporating images acquired under different illumination conditions (i.e., photometric stereo). Introducing photoclinometry into photogrammetric processing can therefore effectively increase the achievable resolution of the mapping result while maintaining its overall accuracy. This research presents an integrated photogrammetric and photoclinometric approach for pixel-resolution 3D modelling of the lunar surface. First, photoclinometry is interacted with stereo image matching to create robust and spatially well distributed dense conjugate points. Then, based on the 3D point cloud derived from photogrammetric processing of the dense conjugate points, photoclinometry is further introduced to derive the 3D positions of the unmatched points and to refine the final point cloud. The approach is able to produce one 3D point for each image pixel within the overlapping area of the stereo pair so that to obtain pixel-resolution 3D models. Experiments using the Lunar Reconnaissance Orbiter Camera Narrow Angle Camera (LROC NAC) images show the superior performances of the approach compared with traditional photogrammetric technique. The results and findings from this research contribute to optimal exploitation of image information for high-resolution 3D modelling of the lunar surface, which is of significance for the advancement of lunar and planetary mapping.


INTRODUCTION
Lunar missions and researches such as rover explorations and geologic analysis require the support of high resolution and accurate 3D surface models, e.g., digital elevation models (DEMs).These surface models are usually obtained by processing stereo images through photogrammetry.With a pair of images acquired under suitable imaging geometry, dense conjugate points are matched through image matching algorithms.The 3D positions of these matched points can be accurately determined through photogrammetric space intersection (Wu et al., 2014).It is then apparent that the spatial resolution of the surface model is largely influenced by the number and distribution of conjugate points identified.In situations where the images are poorly textured or are acquired under significant illumination differences, dense image matching may fail to produce dense and reliable matches; hence leading to limited effective resolution and loss of surface details.On the other hand, photoclinometry (i.e., shape from shading) presented by Horn (1990) and Kirk (1987) is capable of generating pixel-level resolution surface models.When illumination (e.g., sunlight) hits the lunar surface, it interacts with the surface geometrically and physically and is partially reflected.The reflected radiance is then captured by the imaging sensor as pixel intensity.Photoclinometry makes use of this relationship to estimate surface slopes of each pixel and reconstructs 3D surface models based on pixel-wise slopes.Photoclinometry preserves local details while being unstable in coarse scale and it synergizes well with stereo photogrammetric processing where coarse scale accuracy can be ensured.

Algorithms
integrating stereo photogrammetry and photoclinometry for high-resolution lunar surface modelling have been explored and developed (Grumpe et al., 2014;Wu et al., 2017).These algorithms mainly focus on applying photoclinometry on the derived products (e.g., a DEM) from conventional approaches such as stereo photogrammetry and laser altimetry.The possibility and strategy to utilize photoclinometry in improving the performance of photogrammetric processing can be further explored.This paper presents an integrated photogrammetric and photoclinometric approach for high-resolution lunar surface modelling.In the presented approach, photoclinometry is involved in image matching by estimating the underlying slope and surface reflecting properties (i.e., intrinsic albedo) of the images.Intrinsic albedo is generally unaffected by illumination directions.As a result, more dense and robust matching can be achieved.3D positions of the matched conjugate points are photogrammetrically determined and an improved point cloud can be generated based on the increased number of conjugate points from photoclinometric processing.Improved pixelresolution 3D surface models can hence be obtained by using the aforementioned algorithms (Grumpe et al., 2014;Wu et al., 2017) with an improved photogrammetric 3D model from the proposed approach.This paper first describes the integrated photogrammetric and photoclinometric approach, followed by experimental analysis of the approach with real data using LROC NAC images.Finally, concluding remarks are presented and discussed.

Overview
The principle of the proposed approach is to use photoclinometry to increase matching points and to use stereo photogrammetry to derive their 3D positions.Figure 1 presents the overall workflow of the proposed approach: The input image pair is forwarded to photoclinometric processing separately to generate new image pairs.The new image pairs represent the albedo and slope (i.e., tangent) of the underlying terrain.Together with the original image pairs, three image pairs are made available.The three pairs are then matched separately to produce conjugate points.Repeated and invalid conjugate points are filtered out during the process.The combined matching results can be treated as ordinary conjugate points; their 3D locations are then computed by space intersection which is a standard photogrammetric approach.As a result, all the 3D points derived at this stage are photogrammetrically determined; doing so ensures a consistent accuracy and uncertainty pattern for all computed 3D points.Afterwards, the 3D model is used as constraints in photoclinometric reconstruction algorithms such as Grumpe et al. (2014) and Wu et al. (2017) for pixel-resolution modelling.This strategy derives additional information (i.e., conjugate points) from photoclinometric processing and extends its application spectrum.

Photoclinometric processing for deriving albedo information
Image pixel intensity can be modelled as the reflected energy from the illumination source due to surface geometry and albedo of the surface as follows: Where I is the intensity of an image pixel, A refers to albedo and G(p,q) is the geometric reflectance based on a reflectance model and surface gradients.p and q represent the surface gradients in x-and y-direction (i.e.,   and

𝜕𝑦
) respectively.In the proposed approach, photoclinometry is used to estimate the intrinsic albedo map as well as the slope map of each image separately in order to facilitate matching.This is an under constrained problem and is achieved by first estimating the albedo by kernel operations and use the estimated albedo to estimate the underlying slope.
Kernel-based albedo estimation has been proposed and implemented by Lee and Rosenfeld (1983) and Kirk (2003).
The former (Lee and Rosenfeld, 1983) assumes a smooth and implicitly convex surface and iteratively estimates image albedo based on a lambert reflectance model and equation ( 1); it has a characteristic that the albedo of an object will gradually converge towards the brightest image value of the object and is unstable when the underlying shape is non-convex.The later (Kirk, 2003) uses a smoothed version of the image as albedo for photoclinometric processing which is much more stable and straightforward but leaving coarse scale traces of shading.The proposed albedo estimation in this research is an integrated extension of the aforementioned methods.As image intensity is a ratio of incoming and reflected energy ranging from 0 to 1, assuming a pixel must be lit up by illumination (i.e., I > 0), there is a physical limit for both albedo and geometric reflectance (Wu et al., 2017): Equation ( 2) implies that both albedo and geometric reflectance must be at least equal to the image intensity, hence for a dark pixel, the valid range of albedo is larger than that of a bright pixel.In other words, a bright pixel has lower uncertainty and hence is more believable.Based on above, a kernel can be developed which favours brighter surrounding pixels and suppresses darker surrounding pixels.This resembles the characteristic observed in Lee and Rosenfeld (1983)'s approach while keeping the stability of Kirk ( 2003)'s approach.Note that the image I mentioned here does not refer to the input image but an initial estimate of albedo image.
The estimation is performed on multiple resolutions.First of all the highpass component between two consecutive resampling levels is extracted according to Figure 2: The input image is down-sampled to two consecutive levels n and n-1, where level n has a lower resolution than level n-1.The image downsampled at level n is up-sampled again to match the resolution of level n-1; although both resampled images have the same resolution, the one directly down-sampled from input image contains more details than the other.The image highpass component is extracted by the difference between the downsampled image of level n-1 and the up-sampled image of level n.
The process of hierarchical albedo estimation is described in Figure 3: The albedo map of previous level n is first up-sampled to level n-1; the highpass component is added to it and initialized the albedo map of the current level.Albedo map at the current level is then estimated by the proposed kernel processing.This process requires an initial albedo estimate at the coarsest level and it is achieved by the following equation: (5) From Equation ( 5), the initial albedo A0 is obtained by dividing the mean image intensity value I0 by the reflectance value G(p0,q0) of mean surface gradients p0 and q0.p0 and q0 are predefined by assuming certain mean underlying surface gradients such as a flat plane (i.e., p=q=0).The resulted albedo map needs to be scaled in order to be used for slope estimation and currently, this is performed by manually selecting an arbitrary scalar.
Figure 2. Workflow for computing the highpass component between two resampling levels Figure 3. Algorithm for hierarchical albedo estimation

Photoclinometric processing for deriving slope information
The underlying slope can be estimated by using the resulted albedo map using Equation (1).Although there are methods to solve for both p and q in one observation (Grumpe et al., 2014;Wu et al., 2017), they usually require additional information or constraints that may lose optimal image details for matching.Instead, this paper solves for the normal (i.e., negative tangent) along the illumination direction which plausible solutions can be determined (Horn, 1990).Equation (1) will be rewritten as: The reflectance value G(s) is now solely determined by the normal s along illumination direction and can be optimized iteratively or approximately non-iteratively.For example, if the Lambert reflectance model is used such that: where i is the incidence angle between the surface normal and illumination vector as shown in Figure 4.The surface normal s can be computed by the following equation: where i0 is the incidence angle between the Z-axis and the illumination vector.Equation ( 8) implicitly prefers surface normal closer to Z-axis by subtracting i0.As surface normal is the negative of surface tangent.Slope maps (i.e., tangent maps) can be obtained by computing -s in Equation ( 8).By this stage, both albedo maps and slope maps are generated and can be used for photogrammetric processing.

Integrated photogrammetric and photoclinometric processing
Photogrammetric processing for 3D surface modelling in the proposed approach is straightforward, i.e., through image matching and space intersection of matched conjugate points.
As the albedo and slope maps are now treated as new image pairs, there are now a total of three pairs of images referring to the same scene for matching: (1) original input image pair which describes the combined albedo and shading information; (2) albedo image pair which is unlikely to be affected by illumination differences; and (3) surface slope image pair which describes the underlying 3D information of the scene.Currently, they are separately matched and the matching results are combined as a whole set; this strategy has an advantage of flexibility where ordinary matching algorithms and tools such as SIFT can be directly used with least modifications.The combined conjugate point set requires a filtering process to remove repeated points and is then ready for photogrammetric processes.The combined point set will have more and denser points than the result from original images alone.When the exterior orientation of the images are known, the X, Y and Z coordinates of each of the conjugate point pairs can be determined by space intersection using the well-known collinearity equations.The derived photogrammetric 3D model will then be forwarded to shape-from-shading algorithms (Grumpe et al., 2014;Wu et al., 2017) for pixel-resolution refinement.

EXPERIMENTAL ANALYSIS
The performance of the proposed approach is validated using real lunar surface images.As presented in Figure 5, the stereo image pair will undergo the proposed routine described in section 2 to produce conjugate point sets.In this paper, the popular SIFT matching algorithm is employed for image matching of the stereo pair.In order to facilitate a more quantitative analysis, a reference DEM is employed to filter out undesirable points.This is achieved by comparing the vertical error between the 3D points computed by space intersection of conjugate points and that of the reference DEM; vertical error larger than a predefined threshold is discarded.The experimental analysis will present and analyse the number of conjugate points and the corresponding DEM from the proposed approach and compare them with the conventional approach where only the original image pair is matched.
Figure 5.General procedure for experimental analysis

Dataset
A stereo pair of LROC NAC CDR (calibrated data record) images as shown in Figure 6 is used for experimental analysis.The two images were acquired within about two hours and hence their illumination directions are very similar; the resolution is about 0.5 m/pixel.The reference DEM of the same region is obtained from the LROC archive (http://lroc.sese.asu.edu/archive).The reference DEM is produced according to Burns et al. (2012) with the same stereo pair presented here; the reference DEM has a spatial resolution of 2 m/pixel.The threshold for this dataset is set to 2 m, which means points with vertical error larger than 2 m with respect to the reference DEM will be regarded as invalid.

Experimental results
Figure 7 presents the albedo maps (Figure 7a) and slope maps (Figure 7b) derived from the original image pair in Figure 6.Bright values in slope maps mean surface facing away and vice versa.Because the slope maps represent surface gradients along illumination direction, they have apparent similarity with the original image pair; however, they possess different pattern in fine detail which allows image matching algorithms to obtain a different result.The albedo maps in Figure 7a show no apparent signs of shading.Bright regions are usually concentrated around the craters; this is because of the ejecta that is associated with the craters and the existence of small rocks and boulders inside them.For example, Figure 7c reveals a bright ring within the largest crater; they are mostly small rocks similar to those in Figure 7d which produces relatively strong albedo changes.Figure 8 presents the estimated albedo over a region with ejecta.
It is apparent in Figure 8 that the bright regions radiating outwards of a crater are likely to be crater ejecta.Although its corresponding albedo map does not visually resembles the outline of the ejecta, the intensity profile (lower dashed line in Figure 8) and its corresponding albedo profile (upper solid line in Figure 8) reveals that the albedo change is actually captured in the albedo map.The albedo profile presented here is scaled by 32767 to match the data type of the image which is signed 16 bit integer.9 summarises the number of conjugate points matched on each image pair.It is noted that the number of points found in original image pair is similar to that in slope maps while the number of points found in albedo maps is about half of the others.This is because the original image and slope maps preserve most of the details while albedo maps generally preserve features with apparent albedo changes.The absence of shading in albedo maps also suppresses interest point extraction; however, the points found in albedo maps will be better associated with features.The total number of unique matched points presented in the last row of Table 1 is more than a double of matching results from only the original image pair; this indicates that the combined approach significantly increases the population of conjugate points.For all pairs, the number of valid points (i.e., vertical error within 2m) is about 60% of total matched points; this implies that the matching accuracy is consistent throughout all combinations.Figure 10 presents the 3D views of DEMs produced by the conventional approach (i.e., matching on the original image only) and that from the proposed approach.Although they are visually similar, the DEM from proposed approach possesses more topographic details; for example, the close up of the DEMs in Figure 10a reveals that the proposed approach generates better shapes along the edge of the Rima and for the concentrated groups of craters.Figure 10b presents a profile comparison of a crater that is not recovered by the conventional approach but preserved in the proposed approach.This is because very few points were matched from the original image pair within the crater due to problems such as shadows, while there are much more matched points within the same crater from the photoclinometrically derived images.Hence the proposed approach is able to obtain more accurate points where conventional approach may fail to achieve.

CONCLUSIONS AND DISCUSSION
This paper presents an integrated photogrammetric and photoclinometric approach for pixel-resolution 3D modelling of the lunar surface.Albedo and slope maps derived from the original images through photoclinometry are used in image matching together with the original images.The experimental results indicate a significant increase in unique and accurate conjugate points where conventional approach may produce unsatisfactory results.This leads to a more detailed surface model and can be further processed to pixel-resolution through photoclinometric reconstruction algorithms such as Grumpe et al. (2014) and Wu et al. (2017).Currently, the slope maps describe slopes along illumination direction which makes them work best when the illumination condition of the stereo pair is similar.They can be more robust by introducing additional constraints of photometric stereo (i.e., photoclinometry from multiple images and illumination conditions) such as Woodham (1980), Heipke (1992) and Liu et al. (2018) to minimize the influence of illumination differences.A more sophisticated matching approach where image, albedo and slope are simultaneously utilized can be explored based on the strategy presented here.This paper demonstrated that photoclinometry and its concepts can be extended into other applications and it is believed that a wider spectrum of applications can be explored and developed with the use of photoclinometry.

Figure 1 .
Figure 1.Workflow for the proposed approach 3)Equation set (3) describes the mathematical process of the albedo estimation kernel.First of all, a 3-by-3 kernel slides over the image; matrix D records the intensity difference  , between the centre pixel  , ;  = 2and it's surrounding eight neighbours  , .The weighting function (, ) enlarges the difference by a power of k if the difference is positive (i.e., , >  , ) and returns zero otherwise.() is the mean of surrounding intensity difference weighted by (, ), in the case of local maxima, sum of all (, ) is zero and () will directly return zero.Albedo is then computed by the following:  , =  , + () (4)

Figure 4 .
Figure 4. Illustration of the relationship between the surface slope and illumination source

Figure 6 .
Figure 6.LROC NAC CDR stereo pair for experimental analysis and their brief information Figure 7. Albedo maps and slope maps of the image pair derived by the proposed photoclinometric process and particular close up of the experimental dataset b) Profile comparison of a selected crater Figure 10.3D views of DEMs of the experimental dataset

Table 9 .
Number of conjugate points for each image pair