Shape-from-shading Refinement of LOLA and LROC NAC Digital Elevation Models: Applications to Upcoming Human and Robotic Exploration of the Moon

The Lunar Reconnaissance Orbiter (LRO) has returned a wealth of remotely sensed data of the Moon over the past 15 years. As preparations are under way to return humans to the lunar surface with the Artemis campaign, LRO data have become a cornerstone for the characterization of potential sites of scientific and exploration interest on the Moon's surface. One critical aspect of landing site selection is knowledge of topography, slope, and surface hazards. Digital elevation models derived from the Lunar Orbiter Laser Altimeter (LOLA) and Lunar Reconnaissance Orbiter Camera (LROC) instruments can provide this information at scales of meters to decameters. Shape-from-shading (SfS), or photoclinometry, is a technique for independently deriving surface height information by correlating surface reflectance with incidence angle and can theoretically approach an effective resolution equivalent to the input images themselves, typically better than 1 m per pixel with the LROC Narrow Angle Camera (NAC). We present a high-level, semiautomated pipeline that utilizes preexisting Ames Stereo Pipeline tools along with image alignment and parallel processing routines to generate SfS-refined digital elevation models using LRO data. In addition to the present focus on the lunar south pole with Artemis, we also demonstrate the usefulness of SfS for characterizing meter-scale lunar topography at lower equatorial latitudes.


Background
The Lunar Reconnaissance Orbiter (LRO; Chin et al. 2007) has been in continuous operation around the Moon since 2009 and is currently in its fifth extended science mission (ESM5).Data returned from LRO have revolutionized our understanding of the Moon with near-global laser altimetry, radar, thermal emissivity, and high-resolution visible imaging coverage.In anticipation of the first crewed lunar landing in over 50 years with the Artemis campaign, interest in high-priority lunar landing sites has increased dramatically (e.g., Lunar Exploration Analysis Group 2016Group , 2018;;National Research Council 2007, 2011), particularly those in the South Polar Region (SPR; ³84°S) where the first Artemis missions are targeted (National Space Council 2020).Detailed knowledge of landing site geology, topography, and terrain hazards will be critical to ensuring safe operations for future human and robotic missions to the lunar surface (Artemis III Science Definition Team 2020; Lunar Critical Data Products Specific Action Team 2021).The topography of the Moon is currently best understood through altimetric ranging with the LRO Lunar Orbiter Laser Altimeter (LOLA; Smith et al. 2010) and targeted digital terrain models derived from Lunar Reconnaissance Orbiter Camera (LROC; Robinson et al. 2010) stereo imaging.
("DEM," digital elevation model, and "DTM," digital terrain model, are interchangeable; we use "DEM" as a general term and "DTM" to refer specifically to NAC DTM products.)LOLA ranging provides highly accurate vertical height information but relies on interpolation to fill gaps between point shots (Barker et al. 2021).LROC Narrow Angle Camera (NAC) DTM products use stereo photogrammetry, which is limited to ∼3-4x coarser resolution than their source image pairs and does not perform as well in areas of low image contrast or shadows (e.g., Henriksen et al. 2017).
Shape-from-shading (SfS), also known as photoclinometry, is a technique used to interpret surface relief from the intensity of light in an image.SfS can theoretically approach a resolution in topography equivalent to the input images themselves, which is typically better than 1 m per pixel (mpp) using LROC NAC.SfS is also much more tolerant of the extreme and highly variable illumination conditions at the poles, which makes it a particularly powerful tool in cases when stereo photogrammetry is impractical.The fundamental assumption of SfS is that for a surface with known reflectance, the intensity of reflected light should scale as a function of the incidence angle between the light source and the surface (e.g., Van Diggelen 1951;Ikeuchi & Horn 1981;Kirk 1987;Horn 1990;McEwen 1991;Kirk et al. 2003).Thus, a surface will theoretically be brightest when the incidence angle is perpendicular and darkest when the incidence is parallel (i.e., in shadow).While the concept of SfS for the Moon has been treated several times in the literature (Lohse et al. 2006;Archinal et al. 2011;Grumpe & Wöhler 2011;Grumpe et al. 2014;Wu et al. 2018;Chen et al. 2022), we use the built-in SfS functionality within Ames Stereo Pipeline (ASP), an open-source toolset developed by NASA for the generation of DEMs from planetary image data (Alexandrov et al. 2023) and the current state of the art for ongoing SPR landing region characterization.The ASP SfS algorithm recovers topography by performing a least-squares minimization of a cost function over multiple camera images k (Alexandrov & Beyer 2018): where h(x, y) is the height of the terrain at map-projected coordinates (x, y); I k is the kth camera image obtained by projecting 3D points from the terrain through a camera with known position and orientation relative to the ground; and T k , A, and R k are the kth image exposure, terrain-dependent albedo, and reflectance, respectively.The first term collectively constrains the brightness of each pixel; the second two terms are additional constraints on terrain smoothness (μ2 with the sum of squares of all second-order partial derivatives ∇ 2 h) and deviation from the reference DEM height (λ 2 with the reference DEM height h 0 ).The ASP SfS implementation uses a lunar Lambertian reflectance model (McEwen 1991;Lohse et al. 2006), which is modified from a simple Lambertian model that assumes reflectance is a function of the cosine of the incidence angle.Shadow thresholds are introduced to avoid the misinterpretation of shadowed areas as parallel to the incidence angle and thus extremely steep; these "holes" are filled by replacing them with values from the reference DEM (Alexandrov et al. 2023).The minimization problem is discretized using a finite difference method and employs Googleʼs Ceres solver (Agarwal et al. 2022) for the optimization step.
The ASP SfS implementation continues to undergo iterative updates as its usefulness for improving lunar topography models is explored in greater detail.In this regard, we present a high-level, semiautomated processing pipeline for refining digital elevation models of the lunar surface with SfS that can be implemented using only existing tools from the most recent version of ASP (3.3 as of this writing; Alexandrov et al. 2023) and USGS Integrated Software for Imagers and Spectrometers (ISIS; Laura et al. 2023).Taking full advantage of the parallel processing that we describe requires hardware that may only be available through high-performance computing clusters, but the basic pipeline can be implemented on any modern desktop computer.Applications of SfS to the Moon using ASP and other tools have thus far focused on human exploration of the SPR, mostly as part of NASA Artemis internal science team activities (e.g., Bertone et al. 2023;Kirk et al. 2023;Lawrence et al. 2023;Petro et al. 2023).Thus, we aim to highlight the usefulness of LRO-derived SfS data products for the characterization of meter-scale topography of human and robotic landing sites both in the SPR and at lower equatorial latitudes.Below, we present our shape-from-shading methodology and two examples of SfS applied to different sites of scientific and exploration interest across the Moon.

Shape-from-shading Method
The SfS method requires a reference DEM for registration to a geodetic datum and thus map projection and incorporation into a geographic information system (GIS).We use publicly available LOLA LDEM (Barker et al. 2021) and NAC DTM (Henriksen et al. 2017) products.The latest 5 mpp LOLA LDEM data sets of the SPR are available from the Planetary Geodesy Data Archive (PGDA)1 ; we browse and select NAC DTM data sets from the Planetary Data System (PDS) Orbital Data Explorer. 2Full URLs for the data are provided at the end of the paper.We use the LROC Image Search Tool3 to select all raw (experimental data record, EDR) NAC image products that fall within the lat-lon extent of the reference DEM.For sites in the SPR, the total number of overlapping images is significantly greater, so we limit our search to those data acquired after LROʼs periapsis was lowered in 2015 to allow better access to the south pole (corresponding approximately to data release 24 onward).The DEM products are upsampled using cubic convolution in Geospatial Data Abstraction Library (GDAL)4 to match either the resolution of the associated orthoimage, in the case of NAC DTMs, or the highest resolution among NAC images that overlap with the LOLA LDEM, resulting in upsampled reference DEMs on the order of ∼0.5-1 mpp.The ASP SfS algorithm is capable of parallel processing by splitting large DEMs into tiles.We choose instead to manually parallelize our processing to have greater control over the exact size and extent of tiles, typically 600 × 600 pixels with a 100 pixel overlap, which are later mosaicked into a single SfS-refined DEM using dem_mosaic in ASP.In most cases, we crop the original DEM to a region of interest of size m × n tiles before upsampling.Each NAC EDR product is converted to a USGS ISIS-readable image cube and radiometrically calibrated using the standard ISIS tools lronac2isis, spiceinit, lronaccal, and lronacecho.The full set of overlapping NAC image cubes is then bundle adjusted using the ASP bundle_adjust tool; this corrects camera position and pointing discrepancies between images.
SfS requires exact pixel-level coregistration between input image frames and the reference DEM.Bundle adjustment, while necessary, often leaves residual alignment errors that require correction.After the initial bundle adjustment step, we run the ASP mapproject tool on the NAC images to generate cropped, map-projected versions of the NAC images that match the extent of each smaller DEM tile.The recently implemented ASP image_align tool performs 2D interest point matching to align a source image to a reference image.The tool is capable of performing alignment transforms with 2 (translation), 4 (rigid), 6 (affine), or 8 (homography) degrees of freedom.However, in order to be georeferenced, the transforms are limited to 4 degrees of freedom that describe 3D movements about the bodyʼs center.In ASP, this type of reference frame is known as "Earth-centered Earth-fixed," or ECEF, but it is interchangeable with geocentric coordinate systems for any planetary body, including the Moon.All of the image alignments described here are referenced to the standard lunar spheroid with radius 1737 km.We found that adding degrees of freedom beyond simple translation made little difference in the accuracy of the final alignment, so we use x − y translation only to prevent the introduction of unnecessary image distortion that could affect the accuracy of the SfS solution.
In addition to aiding in parallel processing, we found that image_align performs much better with individual tiles than across whole images, where in the former case the required adjustments are relatively small.Even though the same image may be adjusted slightly differently for each tile, the SfS end product will still be self-consistent as long as the individual NAC images are aligned to the tiled DEM.We run a translation transform on the cropped, map-projected NAC images for each tile separately.For the NAC DTMs, we align each input NAC image to the representative orthoimage tile, since this is already coregistered with the DTM that will be used for SfS refinement.
For LOLA, we use controlled NAC mosaics of the sites of interest in place of an orthoimage.We found that adjusting the interest point inlier threshold upward was necessary in cases when the default value failed to identify a sufficient number of matches.The alignment is applied to the NAC images using the ECEF transform in an additional bundle adjustment and map projection step.
Due to varying illumination conditions, frames with bad data, and other factors, the alignment transform occasionally outputs spurious results.To correct for this, we implement a second automated filtering step using a structural similarity index (SSIM; Wang et al. 2004) to determine which aligned images match the orthoimage or controlled mosaic within a given threshold.For each SfS region, we manually separated correctly and incorrectly aligned images for certain tiles and tested their SSIM index distribution against a filtered subset with an arbitrary SSIM threshold to find the lowest SSIM index that still captured the correctly aligned images.In the final SfS product, the difference between outputs using manually filtered versus SSIM-filtered source images was virtually indistinguishable.
We perform an initial SfS run to calculate exposures for each tile followed by the actual SfS processing run.While SfS can run the Ceres solver to full convergence, we found that the first run solution is often sufficient and that further runs resulted in relatively little qualitative difference in the output; a single run for each tile also significantly decreases processing times.We also experimented with varying the μ and λ parameters in Equation (1).A lower value of μ results in a qualitatively sharper solution with less smoothing but tends to emphasize processing artifacts; a value of μ = 0.08 strikes a balance between the potential loss of legitimate height information with artifact removal, where they exist.Varying the λ parameter made no discernible difference in the solution, so we maintain the default value of λ = 0.001.We use the default shadow threshold (∼floating-point pixel brightness value) of 0.005.
SfS outputs a final DEM along with a solution for terraindependent albedo.In addition to these output products, we generate the following ancillary data products: (1) hillshaded DEMs using the Sun azimuth and elevation of the respective orthoimage or controlled mosaic; (2) absolute height difference maps, which show how the SfS solution differs from the reference; and (3) slope and aspect maps.

Example Sites with Applications to Science and Exploration
Below we present two SfS processing examples, one each using a LOLA LDEM and NAC DTM as reference.The LOLA LDEM site is of interest for upcoming Artemis landings in the SPR, while the NAC DTM site showcases the possibility of deriving SfS solutions for regions of scientific and exploration interest at equatorial latitudes.

LOLA LDEM: Malapert Massif (Artemis)
Thirteen candidate landing regions were announced in 2022 for Artemis III and future human lunar landings in the SPR (NASA 2022).A number of the 13 regions and other sites within the SPR now have coverage with 5 mpp LOLA LDEM gridded topography mosaics (Barker et al. 2021).For an example of SfS refinement of an LDEM, we highlight the Malapert Massif landing region that was chosen for the Artemis III Geology Team EVA planning exercise and has also garnered significant interest as a candidate for future human landings (e.g., Basilevsky et al. 2019Basilevsky et al. , 2023;;Longo 2023).
The nominal landing site for the EVA planning exercise is near the point of maximum solar illumination at 85.964°S, 357.681°E on a ridge near the summit of Mons Malapert.The elongated profile of Mons Malapert is thought to be an uplifted remnant of the South Pole-Aitken (SPA) basin rim (Garrick-Bethell & Zuber 2009; Krasilnikov et al. 2023) with an age of ∼4.2 Ga (Basilevsky et al. 2019).The Malapert Massif region would be an ideal location for understanding the earliest geologic history of the Moon and likely samples ancient crustal material from the SPA impact (Krasilnikov et al. 2023).The terrain on either side of the ridge is dominated by steep slopes and the presence of elephant hide terrain (EHT; e.g., Kreslavsky et al. 2021), which could present trafficability challenges for walking traverses or piloted rovers at the site.
We based our processing extent upon the required 2 km EVA radius from the landing site as described in the A3GT NOFO document (NASA ROSES 2022).Since no predefined orthoimage exists, we generated a 4 × 4 km clipped version of the 5 mpp Site23 LDEM product resampled to a target resolution of 1 mpp, or an improvement of ∼5x.We then used the controlled NAC mosaic of Malapert Massif (Henriksen et al. 2023) for the image alignment step.We split SfS processing into 8 × 8 tiles for a processing extent of 4000 × 4000 pixels (16 km 2 ).The automatic alignment and filtering steps resulted in between 4 and 25 usable NAC images per tile that were used to compute the SfS solution.
Both the reference LOLA LDEM (Figures 1(a ).Knowledge of these hazards will be critical for Artemis Human Landing System (HLS) design and testing as well as planning for extravehicular activity.Both solutions indicate that slopes off the ridgeline are convex in profile, increasing slowly at first (∼10°-15°) and becoming increasingly steep away from the ridge (∼20°-25°).The SfS solution also removes LOLA interpolation artifacts that appear as diagonal lines in the LDEM product.The absolute height difference map shows that six of the 64 tiles failed to generate an SfS solution (Figure 3(a)).While these could be easily reprocessed manually, it demonstrates a potential shortcoming of the automated alignment and filtering steps in these instances.
One advantage of LOLA is coverage even in permanently shadowed regions (PSRs), which are inaccessible to any topographic data product that relies on visible imaging such as photogrammetry or photoclinometry.The SfS solution smooths the underlying LOLA data to fill gaps in the SfS solution (e.g., Figure 3(a)), and the seams between the layers often result in the introduction of slope artifacts; future planning efforts should take this limitation into account.This would significantly extend the known lifetime of volcanic activity on the Moon beyond previously established estimates (e.g., Hiesinger et al. 2011).Other hypotheses suggest that a highly porous substrate has masked impact craters and other indicators of weathering and that Ina is closer in age to the ancient maria surrounding it, ∼3.5 Ga (Qiao et al. 2017;Wilson & Head 2017;Qiao et al. 2020).Thus, the determination of the true age of the Ina deposits would be key to unraveling the history of lunar mare volcanism and the Moonʼs thermal evolution.
In order to address these questions, the Dating an Irregular Mare Patch with a Lunar Explorer  launch no earlier than 2027 (Anderson et al. 2023).The nominal landing site for DIMPLE is Mons Agnes, a large mound located near the eastern edge of Ina.Mons Agnes will allow access to both mound and floor material for a small rover that will return samples to the lander for radiometric dating.The mission profile is similar to those previously proposed for Ina (Stopar et al. 2019;Qiao et al. 2021) that would seek to determine the age and origin of the Ina deposits through in situ analysis.
We generated a 2 × 3 tile SfS solution for Mons Agnes and its immediate surroundings (centered at 18.657°N, 5.337°E) using the 2 mpp NAC_DTM_INACALDERA1 product as reference with an orthoimage and SfS processing resolution of 0.6 mpp, or an improvement of ∼3x.A size of 2 × 3 tiles results in a total processing extent of 1100 × 1600 pixels or 660 × 960 m (∼0.63 km 2 ).Between 18 and 28 individual NAC images per tile were used to compute the SfS solution after automated alignment and filtering.
Both the reference NAC DTM (Figure 4    The data products presented above can be used for both qualitative and semiquantitative calculations of DEM quality, which we demonstrate for Mons Agnes in Figure 7.Such techniques include (but are not limited to) plotting pixelwise slope or height differences between the two solutions and looking for outliers (Figure 7   .Other outliers can be explained by a slight horizontal offset between the two solutions that is visible in the bias along westward-facing slopes in the absolute height difference map (Figure 6(a)).This is due to an underlying difference in alignment between the NAC DTM and the orthoimage that we use as a basis for image alignment in the SfS solution.Different data products will always have some inherent spatial variability at such fine scales, and users should be aware of these discrepancies when choosing source data sets from which to derive topographic solutions.We chose to use NAC orthoimages as our ground truth since SfS relies upon shading information from individual image frames, as opposed to stereo photogrammetry that uses interest point matching across multiple images to derive 3D match points that may not necessarily correspond to the same point on the ground as indicated by either image in a stereo pair.

Conclusions
Shape-from-shading is a powerful tool for refining digital elevation models up to an effective resolution that approaches that of the best image resolution of the surface.Such highresolution DEMs will be critical components of the site selection and characterization necessary to perform safe landings and mission operations on the Moon.As work on the computational technique behind SfS has progressed in preparation for returning humans to the Moon with the Artemis campaign, we have taken advantage of preexisting tools and processes to create a high-level, semiautomated pipeline for generating SfS-refined DEMs of the Moon.Three-dimensional image alignment transforms combined with structural similarity (SSIM) filtering result in highly robust input data sets for SfS with minimal human involvement required, and parallel processing and image tiling aid in reducing processing times.Using this pipeline, we have demonstrated the usefulness of SfS for any site where reference DEMs can be found; beyond targeted NAC DTMs, gridded global topography mosaics from LOLA and LROC Wide Angle Camera (WAC) exist at a range of resolutions for most of the lunar surface.Continuing to scale up topographic data processing efforts will be key to developing a robust lunar data infrastructure as part of a sustained human presence on the Moon.

ORCID iDs
Benjamin D. Boatwright https://orcid.org/0000-0002-0370-2599James W. Head https:/ /orcid.org/0000-0003-2013-560X ), (c)) and SfS solution (Figures 1(b), (d)) show the overall topographic signature of the site: a central east-west ridgeline with primarily north-and south-facing slopes.(NAC images of the site additionally show that only the northern or southern slopes are ever illuminated at one time while the top of the ridge remains mostly in sunlight.)The LOLA LDEM gives a general sense of surface slopes (Figures 2(a), (c)), but the SfS solution (Figure 2(b)) reveals areas within 100 m of the proposed landing site that are dominated by short-wavelength roughness and small impact craters (Figure 2(d)) that correspond to visible hazards in the NAC mosaic (Figure 2(e) 3.2.NAC DTM: Ina Irregular Mare Patch (PRISM) Irregular mare patches (IMPs) represent a lunar geologic enigma.On the basis of crater size-frequency distributions, optical immaturity, and sharp surface relief, these features appear to indicate extremely young mare volcanism (Schultz et al. 2006; Braden et al. 2014; Fassett & Thomson 2014; Basilevsky & Michael 2021) inside of collapsed summit calderas (pit craters) together with solidified mounds of latestage lava extrusions (Wilson & Head 2017), both of which have absolute model ages of <100 Ma (Braden et al. 2014).
(DIMPLE) investigation will land at Ina as the third Payloads and Research Investigations on the Surface of the Moon (PRISM) selection through NASAʼs Commercial Lunar Payload Services (CLPS) program, due to

Figure 1 .
Figure 1.Cropped LOLA LDEM (a), (c) and SfS solution (b), (d) for the Malapert Massif candidate landing region, centered at 85.964°S, 357.681°E on a ridge near the summit of Mons Malapert.Both products show a central east-west ridgeline with primarily north-and south-facing slopes.Two hillshade images match illumination conditions of the low-Sun controlled NAC mosaic with subsolar longitude 315°((a)-(b), Sun from top left) and 235°((c)-(d), Sun from bottom left), elevation 5°above the horizon.
(a)) and SfS solution (Figure 4(b)) highlight the height differences between Mons Agnes, the surrounding floor, and the exterior of Ina itself.The SfS solution improves upon the textural contrast between Mons Agnes and the surrounding floor.The characteristic roughness length scale of the floor approaches or exceeds the NAC DTM resolution, leading to obscuration from noise in the reference DTM.The same can be observed for the mounds, where an artificial surface roughness is introduced in the NAC DTM (Figure 4(c)); the hillshaded SfS solution (Figure 4(d)) comes much closer to approximating the actual appearance of the mound as viewed in NAC images of the notional DIMPLE landing site (Figure 4(e)).The slope distribution between the NAC DTM and SfS solution is similar, with a mean slope of 5°( Figure 5).Both correctly represent the steep margins of Mons Agnes, which are in excess of 20°, as well as the much gentler westward slope toward the center of Ina (also apparent in the aspect map, Figure 6(b)).The SfS solution inherits a few small errors from the reference NAC DTM, such as the unnatural squared off margin near the southwest corner of Mons Agnes (Figure 5(b), arrow; Figure 7(b)), but makes noticeable improvements on submeter topographic variations that the

Figure 2 .
Figure 2. Slope maps for the LOLA LDEM (a) and SfS solution (b) for Malapert Massif classified in increments of 2°with a cutoff at the +2σ value in the slope distribution (∼24°, mean 14°).A narrower classification of the SfS solution (cutoff 10°, increments of 1°) reveals areas near the top of the relatively flat ridge (insets (c)-(e); boxes in (a)-(b)) that are dominated by short-wavelength roughness and small impact craters that are missed in the LDEM (c) but visible in the SfS (d) and corresponding NAC mosaic (e).Stars in (c)-(e) show the location of the proposed Artemis landing site.

Figure 3 .
Figure 3. Absolute height difference map (a) and aspect map (b) for Malapert Massif.Blank tiles in the absolute difference map indicate areas where the algorithm filled gaps in the SfS solution with underlying LOLA data.

Figure 4 .
Figure 4. Cropped NAC DTM (a) and SfS solution (b) for Mons Agnes, Ina irregular mare patch, centered at 18.657°N, 5.337°E.The hillshaded SfS solution (d) shows improved textural contrast between Mons Agnes and the surrounding floor ((e), NAC orthoimage) by removing photogrammetry processing artifacts present in the NAC DTM (c).Hillshades (c)-(e) match illumination conditions of the NAC DTM orthoimage, subsolar longitude 257°(Sun from bottom left) and elevation 35°.
(a)) and visually comparing SfS hillshades with the same solar illumination conditions as the original source image(s) to identify artifacts (Figure 7(b)).The highest density of pixelwise slope points (Figure 7(a)) falls

Figure 5 .
Figure 5. Slope maps for the NAC DTM (a) and SfS solution (b) for Mons Agnes classified in increments of 2°with a cutoff at the +2σ value of the slope distribution (14°, mean 5°).Inset panels show an area near the eastern margin of Ina that contains high-albedo blocky material (e, NAC orthoimage); the SfS solution (d) resolves the short-wavelength slope variation of this terrain while the NAC DTM (c) shows generally flat slopes.

Figure 6 .
Figure 6.Absolute height difference map (a) and aspect map (b) for Mons Agnes.Large height differences (>2 m; dark red) are due to the slight misregistration between the NAC DTM and SfS solution, which is caused by an underlying offset between the NAC DTM (as viewed via hillshade) and the orthoimage.More modest height differences in the 1-2 m range reflect detections of short-wavelength roughness in the SfS solution, e.g., Figure 5(d).Comparison of the two hillshaded results to the orthoimage (Figure 4(c)-(e)) shows that the SfS solution is qualitatively more accurate to ground truth.