SHAPE AND ALBEDO FROM SHADING ( SAfS ) FOR PIXEL-LEVEL DEM GENERATION FROM MONOCULAR IMAGES CONSTRAINED BY LOW-RESOLUTION DEM

Lunar topographic information, e.g., lunar DEM (Digital Elevation Model), is very important for lunar exploration missions and scientific research. Lunar DEMs are typically generated from photogrammetric image processing or laser altimetry, of which photogrammetric methods require multiple stereo images of an area. DEMs generated from these methods are usually achieved by various interpolation techniques, leading to interpolation artifacts in the resulting DEM. On the other hand, photometric shape reconstruction, e.g., SfS (Shape from Shading), extensively studied in the field of Computer Vision has been introduced to pixellevel resolution DEM refinement. SfS methods have the ability to reconstruct pixel-wise terrain details that explain a given image of the terrain. If the terrain and its corresponding pixel-wise albedo were to be estimated simultaneously, this is a SAfS (Shape and Albedo from Shading) problem and it will be under-determined without additional information. Previous works show strong statistical regularities in albedo of natural objects, and this is even more logically valid in the case of lunar surface due to its lower surface albedo complexity than the Earth. In this paper we suggest a method that refines a lower-resolution DEM to pixel-level resolution given a monocular image of the coverage with known light source, at the same time we also estimate the corresponding pixel-wise albedo map. We regulate the behaviour of albedo and shape such that the optimized terrain and albedo are the likely solutions that explain the corresponding image. The parameters in the approach are optimized through a kernel-based relaxation framework to gain computational advantages. In this research we experimentally employ the Lunar-Lambertian model for reflectance modelling; the framework of the algorithm is expected to be independent of a specific reflectance model. Experiments are carried out using the monocular images from Lunar Reconnaissance Orbiter (LRO) Narrow Angle Camera (NAC) (0.5 m spatial resolution), constrained by the SELENE and LRO Elevation Model (SLDEM 2015) of 60 m spatial resolution. The results indicate that local details are largely recovered by the algorithm while low frequency topographic consistency is affected by the low-resolution DEM.


INTRODUCTION
Lunar DEMs are mainly obtained from stereo photogrammetry or laser altimetry.Stereo photogrammetric processing requires coverage of a target area with stereo images and is able to produce DEMs with best resolutions of about three times of the image resolution.This domain focuses mainly on the imaging geometry while the image photometric content is often used for conjugate matching and imaging pose determination.While stereo coverage is a prerequisite, for some lunar orbiters such as the LRO, high resolution stereo (i.e.LRO NAC) availability is limited due to its orbital and mission design (Tran et al., 2010) and hence hindering the production of high resolution DEMs through stereo photogrammetry.On the other hand, Laser altimetry such as the Lunar Orbiter Laser Altimeter (LOLA) (Smith et al., 2010) produces global and highly reliable topographic information of the target body (i.e. the Moon) for various scientific purposes.However it is also characterized by its large sample spacing, resulting in low resolution DEMs (e.g.1024 pixel-per-degree for LOLA DEMs), and utilizing the derived topographic products for high resolution purposes may introduce significant interpolation artifacts and thus limiting its applications.
Reflectance-based surface reconstruction methods estimate the 3-D surface shape solely based on the image photometric content; such methods are able to utilize the information captured in every pixel of one or more image(s) and reflect them on the resulted topography.Apart from its drawbacks of instability and solution ambiguity, it shows potential as a complement to the existing DEM production method (i.e.stereo photogrammetry and laser altimetry) which are remarkably reliable and accurate on large scales (Kirk et al., 2003).The relatively simple surface albedo complexity of major planetary bodies of interest also favours such techniques to be involved in lunar and planetary mapping.
There are studies and implementations utilizing photometric reconstruction algorithms in lunar and planetary mapping (e.g.Horn, 1977Horn, , 1990;;Kirk et al., 2003;Barron and Malik, 2011), where the terrain is reconstructed based solely on photometric information.Incorporating existing low resolution DEMs into SfS algorithms using single or multiple images (photometric stereo) has also been explored (Grumpe et al., 2014) and encouraging results have been achieved.
In this paper a method is described that refines an existing low resolution DEM (either from laser altimetry data or stereo photogrammetry) to image pixel-level resolution by reflectancebased surface reconstruction techniques (i.e.Shape and albedo from shading, SAfS) utilizing one single image.The resulted DEM details are desirable only if their pixel-wise albedo is also reasonable, therefore we complement the shape recovery algorithm with a pixel-wise albedo estimation algorithm based on regularity assumption.The following Section 2 describes the details of the approach.Experimental results using the LRO NAC imagery (0.5 m/pixel) constrained by the SLDEM2015 (resolution: 60 m) are presented in Section 3. Section 4 provides concluding remarks and a brief discussion.

Framework of the approach
The general structure of the proposed method is outlined in Figure 1.Our algorithm adopts an iterative pyramid (hierarchical) structure throughout the optimization process starting from lowest resolution (i.e.resolution of the input gross DEM) and subsequently increases the resolution by a factor of 2 until it approaches pixel resolution of the image.At each level of the pyramid, the image is resampled to fit to the grid resolution of the DEM and the reflectance can be estimated by utilizing the albedo regularization strategy.Afterwards a least squares optimization is performed to adjust the surface heights (shape) to approach the estimated reflectance, constraints from the refined DEM of previous pyramid level is also introduced to regulate the optimization.Similar to other works (Horn, 1990;Grumpe et al., 2014), we optimize the DEM by a relaxation strategy where the grid node heights are adjusted one at a time assuming other heights are constants; the adjusted height will be updated immediately to the DEM.The final output of the hierarchical processing is a refined DEM and an albedo map with resolutions identical to the image resolution.Detailed procedures are described in the following sections.

Fundamentals of reflectance-based surface reconstruction
Reflectance-based surface reconstruction (e.g., SfS) has been extensively studied for more than half a century (Kirk et al., 2003).Such techniques recover object details based on the relation between the surface gradients (e.g.defined by slope and aspect), the illumination conditions and the recorded image intensity; in other words, reflectance-based surface reconstruction searches for a shape that is able to explain or produce the captured image.
The relationship between photometric image content and topology can be written as (1) where I = image intensity A = albedo G(p,q) = topological reflectance p = surface normal in x direction q = surface normal in y direction The reflectance model G(p,q) builds up the relationship between image photometric content and 3-D topography; these models relate surface gradients (slope and aspect), illumination conditions, surface photometric properties (i.e., albedo) and the resulting intensity captured in the image.Various models can be employed depending on applications and data, such as the sophisticated Hapke model (Hapke, 2002(Hapke, , 2012) ) if measurements in physical units are available (Grumpe et al., 2014) or the empirical Lunar-Lambert model (McEwen, 1991) if only image intensities are available.In this paper the Lunar-Lambert model is used.
Photoclinometry is a 1-D profile-based slope and height estimation algorithm, this algorithm extracts profiles along the illumination direction named "characteristics" (Horn, 1977) and estimates the surface slopes and heights along profiles based on image intensities of the profiles, where the assumption is made that, amongst all combinations of slope and aspect, the slope along the illumination direction has the dominating effect on the resulting image reflectance.There are developments in extending photoclinometry to 2D such as (Kirk et al., 2003), but in this paper we keep the meaning of photoclinometry as 1-D profile-based.
Shape from Shading (SfS) is a 2-D based photometric surface reconstruction approach.Different from profile-based photoclinometry which controls the aspect to be sun-facing, SfS algorithms solve for both slope and aspect to explain the image intensity, thus having a high ambiguity but are more flexible.

Hierarchical shape reconstruction constrained by lowresolution DEM
The regularization of shape reconstruction is achieved by two distinct strategies: implicit regularization by computational and system design and explicit regularization by adding constraints into the optimization.(2) Where [ ] [ ] is the element-wise multiplication.
This means there would be four surface normals and the subsequent reflectance (i.e. the grids in Figure 2) surrounding the center node; holding the surrounding heights as constants by the relaxation scheme retains local topography.The center height is optimized such that the four surrounding reflectances comply with the expected values computed from the piece-wise constant albedo assumption.
The algorithm utilizes the Gauss-Seidel relaxation scheme which sequentially replaces the surface heights with newly computed values.As a result the sequence of computation within the relaxation algorithm has significant effects on the resulting DEM, especially at highly shaded or illuminated regions (e.g.craters, boulders).The refined DEM may apparently deviate from the desired topography due to the selection of the rendering sequence, which thus becomes another implicit shape constraint.In this paper we choose the rendering sequence to be in the same quadrant of the illumination direction: for example if the illumination is from southeast, the rendering sequence would start at bottom-right corner and end at upper-left corner.
We adopt a similar explicit DEM constraint as Grumpe et al., (2014) on the relative depth which requires the low-pass component of a group of surface normals to comply with their counterparts in the gross DEM: where is the low-pass operator over the surface normal.
The algorithm adopts an iterative hierarchical structure throughout the optimization process, starting from lowest resolution.We include the DEM constraint into the algorithm by encouraging the low-pass component of a group of surface normals to approach their counterparts in the previous pyramid level.As we up-sample the DEM by a factor of 2, a DEM cell is split into 2x2 cells in the next level; the distances between the new and old centers are thus equal, and therefore the low-pass operator would become an arithmetic mean of the four surface normals derived from the same cell of the previous level.We deliberately avoid the use of absolute heights to constrain the optimization because we are aware of topographic inaccuracies of the gross DEM (see Figure 5a, the crater is not captured by the gross DEM) that may further flatten the results.
As we base every shape constraint on the refined DEM of the previous pyramid level, we are aware that the result may slowly deviate from the original topography, which is not considered apparent in our experiments; on the other hand, constraining shape using original gross DEM may increase optimization stability while sacrificing some extent of local details.

Locally varying albedo
A terrain cannot be properly refined if the albedo is not reasonably estimated pixel-wise, however such an optimization system with only one image would be under-determined without regularization; therefore we did not bundle albedo optimization with shape refinement, instead, we separately estimate the albedo.While the constant albedo assumption is often plausible in lunar and planetary mapping (Horn, 1990;Kirk et al., 2003), previous works on natural image statistics (Huang and Mumford, 1999) and Shape from Shading (Alldrin et al., 2007;Barron and Malik, 2011) assume a specific statistical behaviour of the albedo, that albedo values tends to be highly concentrated and the difference between adjacent albedo is very small, effectively suggesting a piece-wise constant albedo assumption.
In our algorithm we adopt the piece-wise constant albedo assumption which maximizes surface detail recovery by a reasonable extent.Such assumption is implemented by first taking the log of Equation ( 1) which linearizes the relationship: (4) If we consider the log intensity difference between adjacent pixels, the relationship becomes (5) Applying the piece-wise constant albedo assumption, equation ( 6) can be simplified to: (6) The log reflectance is estimated to satisfy equation ( 7), and we extend the idea to a 2x2 kernel leading to the relationship ) of the reformulated equation ( 1) (8) which may cause numerical instabilities.

Optimization
Assuming a known light source and a known viewer position, we optimize and refine the DEM by the image.The optimization starts by first estimating the albedo and the subsequent estimated reflectance.
At the start of each iteration we compute the reflectance map (i.e.G) of the up-sampled DEM and perform a smoothing over the resulting reflectance, leading to a suppression of possible invalid reflectance values (i.e.negative reflectance) which may appear due to relaxation optimization; afterwards we extract a 2x2 reflectance patch per DEM node and estimate the likely reflectance by satisfying equation ( 8) and the expected reflectance (i.e.E[G]) for the whole region.
The optimization of the DEM is achieved by utilizing the estimated reflectance map and minimizing the error term where the first term refers to reflectance and the second term refers to the DEM shape constraint, while are the weights for each term.
A refined DEM is generated after optimization; it is then upsampled for the next level of the hierarchy until the DEM resolution reaches the resolution of the image.The subsequent new albedo map is then computed from the image and the refined DEM.

Datasets
We tested our algorithm on a typical area on lunar surface (refer to Figure 3).The area was selected by its topological complexity.It contains:


An edge of a rima which splits the region into highland and lowland  A comparatively larger crater located right at the edge of the rima  Groups of small boulders scattered around the larger crater  Groups of small craters all over the region We used the LRO NAC image (0.5 m/pixel) (image pair no.M173246166, downloaded from http://lroc.sese.asu.edu/) and ortho-rectified by the software ISIS3 as the sole input image.
The image is projected to a Mercator system and the region is covered by 1382x1262 pixels.For the low-resolution DEM, the co-registered SLDEM2015 of the area is utilized (http://pdsgeosciences.wustl.edu/missions/lro/lola.htm).SLDEM2015 is the DEM created from LOLA data and images of the Japanese Selenological and Engineering Explorer (SELENE) terrain camera (TC) (Barker et al., 2015).The SLDEM2015 has a resolution of 60 m and covers the testing region with only 12x11 pixels (Figure 3b).The crater spans across a few cells of the SLDEM2015 and its morphology was not properly captured by the SLDEM2015, solving these issues might help to increase the performance of SAfS for low resolution DEM refinement.
Figure 3c shows another DEM (referred to as NAC DEM) with a resolution of 3 m, which was generated from stereo NAC images through photogrammetric process (Wu et al., 2014) and is employed in this study for comparison purposes.

Experimental results
After applying our SAfS algorithm to the testing data, the final refined DEM and albedo map were obtained, which are shown in Figure 4. Figure 5 provides a side-by-side comparison of our refined DEM with the gross SLDEM2015 and the NAC DEM in 3-D visualization.
From Figures 4 and 5, it can be noticed that the refined DEM from our SAfS algorithm shows promising local detail recovery compared with the other two DEMs.We also noticed that the refined albedo map to some extent contains topographic details and slight shades, which implies that a small extent of topographic detail is suppressed by shape constraints and is eventually attributed to the albedo.Figure 6 shows the result of a detailed profile comparison.Two reference lines crossing the crater area were selected, and for each reference line three elevation profiles were interpolated from the three DEMs.The profile comparison presents a shape that deviates from the NAC DEM, especially in Figure 6a where the depth of the crater in the refined DEM is around 10 m shallower than its counterpart in the NAC DEM.Such deviations are likely due to the large resolution differences between the input image and the gross DEM (i.e.0.5 m resolution for the image and 60 m for the SLDEM).
For recovery of topographic details, our SAfS algorithm shows very promising performance.The profile in Figure 6a crosses a boulder at the bottom of the crater.The profile in Figure 6b crosses boulders on the slope of the crater.They are marked with boxes in the figures.These detailed small features are apparent in the input image and the refined DEM obtained with our algorithm, but not in the NAC DEM.

CONCLUSIONS AND DISCUSSION
In this paper we have presented a SAfS algorithm.With a single image and a low-resolution DEM, a refined DEM at image pixel level as well as its corresponding albedo map have produced by the proposed algorithm.
Encouraging results have been shown, especially regarding the recovery of local detail.Areas where very small boulders and craters captured by the image are not presented in both the low resolution SLDEM (60 m) and the higher resolution NAC DEM (3 m) have been reconstructed and preserved in our refined DEM (0.5 m).Slight shadows and shades are incorporated in the albedo map, which will be our next step to improve.
On the other hand, geometric deviations between the refined DEM and the NAC DEM have been spotted.These areas are mainly craters that are covered by only a few cells in the input DEM.We conclude that inaccuracies in the low resolution DEM will cause visible effects in the large scale topography of the refined DEM, as the large scale terrain is largely determined in the initial stages of the hierarchical process.The importance of self-correction at the initial stages is thus emphasized.
Future work will be to further assess the robustness of the algorithm and to allow the algorithm to be able to correct the large scale components of the input DEM before SAfS refinement as well as reducing shading and shadow components that remain in the final refined albedo.

Figure 1 .
Figure 1.Framework of the SAfS algorithm

Figure 2
Figure 2. 3-by-3 Kernel over the DEM and its enclosed surface reflectance (i.e.G 1 to G 4 )At every height node of the relaxation scheme, we extract a 3x3 kernel and only the center node (i.e. the star in Figure2) is 7) with N([*]) as the normal vector formed by the matrix [*].computed using the operator in Equation (3).This strategy avoids the reciprocal component (i.e. Figure 3. (a) The LRO NAC image (0.5 m/pixel), (b) SLDEM2015 (60-m resolution); and (c) NAC DEM (3-m resolution) for comparison purpose.
(a) The refined DEM obtained with our SAfS algorithm and (b) the refined albedo map of the region.3-D views of (a) the refined DEM (0.5 m/pixel) obtained with our SAfS algorithm, (b) the SLDEM2015 (60-m resolution) used as constraint, and (c) the NAC DEM (3-m resolution) used for comparison.
Figure 6.(a) Profile comparison for reference 1, (b) profile comparison for reference 2, (c) local regions of the image, SAfS DEM, and NAC DEM corresponding to the box mark in (a) showing a boulder at the bottom of the crater, and (d) local regions of the image, shaded SAfS DEM, and shaded NAC DEM corresponding to the box mark in (b) showing boulders on the crater slope.