The topography of planetary surfaces provides essential information for a wide field of applications. For example, geomorphologic analysis requires high-resolution Digital Elevation Models (DEMs), but the analysis of hyper-spectral data also benefits greatly from accurate DEMs. To correct hyper-spectral data for photometric effects and thermal emission, a detailed DEM and, in particular, accurate slopes are vital. Furthermore, the planning of rover missions relies on DEMs to assess landing sites and possible hazards associated with steep terrain. In recent years the interest in rover missions to Mars, in particular, has again increased significantly. In 2021, both the National Aeronautics and Space Agency’s (NASA’s) Perseverance Rover and the Chinese National Space Agency’s (CNSA’s) Zhurong rover successfully landed on the Martian surface. The Rosalind Franklin rover of the European Space Agency (ESA) was planned to be launched to Mars in 2022, but has been postponed to a later date. The Martian surface is covered by a thin and non-negligible atmosphere that increases the model complexity of intensity-based reconstruction methods.

2. DEM Generation for Planetary Surfaces

In general, there are three major methodological approaches that are used to generate DEMs of planetary surfaces, i.e., ranging techniques, photogrammetric approaches and shading-based methods. Laser-altimetry samples the surface with laser pulses and the time-of-flight of the photons is translated into range measurements. The resulting DEMs have a high vertical fidelity and an extensive coverage of the planet is achieved. However, the lateral resolution is limited due to a comparatively coarse sampling of the surface (e.g., LOLA: 59 m/pix [2]^[1], MOLA: 463 m/pix [3]^[2]). Photogrammetric approaches, or stereo vision, use two or more images taken from different perspectives to infer a DEM, usually based on bundle adjustment. These methods do not require any physical reflectance model and are known to yield a good absolute height estimate. The state-of-the-art frameworks commonly employed in the planetary community are the Ames Stereo Pipeline [4]^[3] and BAE Systems’ SOCET SET^® [5]^[4]. Both rely on blockmatching to obtain tie-points for the bundle-adjustment procedure. Extensive regions of the planets are covered by regolith, which naturally appears textureless. The lack of texture may cause mismatching, which yields a variety of reconstruction artifacts such as spikes, holes where no matches could be generated and stair-like structures termed pixel-locking [6]^[5]. These stereo artifacts effectively lower the resolution of the DEM to several times below the pixel resolution. Shading-based methods require a reflectance model to connect radiance measurements and surface gradients. They also need proper initialization to ensure convergence. This additional effort is justified by obtaining a DEM of pixel level resolution with a very accurate reconstruction of slopes and heights and the elimination of stereo artifacts, especially in textureless areas (e.g., [7]^[6]). In order to enhance the quality of DEMs, recent approaches have successfully combined low resolution DEMs from photogrammetry and laser-altimetry with SfS, among others, on either the Moon (e.g., [1,8,9,10]^{[7][8][9][10]}), Mercury (e.g., ^[11]) or on Mars (e.g., [12,13,14,15,16,17,18]^{[12][13][14][15][16][17][18]}). The rationale is that the low resolution DEM provides a good absolute height estimate and SfS is used to refine the surface such that the whole procedure combines the advantages of both approaches. Recently, various works have explored machine-learning techniques to directly infer a relationship between measured image intensities and surface height (e.g., [19,20,21]^[19][20][21]). These approaches are purely data-driven and do not incorporate physical information. However, the training procedures implicitly learn the atmospheric conditions present in the training data, but the atmospheric influence has not been explicitly investigated in these works. Consequently, machine-learning approaches will also benefit from investigating the influence of the atmosphere on the reconstruction.

3. Shape from Shading in Planetary Remote Sensing

Galileo Galilei stated that surfaces that are tilted away from the sun appear darker and parts that face the sun appear brighter ^[22].This is probably the first time that surface shading was used to analyze the topography of a planetary body, and it has been used ever since. With the rise of more rigorous physical approaches and an increase in computational power, numerous approaches have been developed that use illumination geometry and shading to quantify the surface slope and topography, primarily of the Moon. In planetary science, these techniques are often termed photoclinometry, and the earliest approaches are given by Rindfleisch ^[23], Wildey ^[24], and Kirk ^[25]. A highly recognized technique from the computer vison community is the SfS method of Horn ^[26], which allows for the integrated recovery of slopes and heights. The surface reconstruction problem is encoded in terms of variational calculus and hence it is solved by minimizing a functional that penalizes the deviations between the shaded surface and the input image. To ensure integrability of the estimated surface, an additional regularization term is introduced. Further constraints were introduced by Shao et al. ^[27] to elastically tie the surface to a low frequency constraint surface, which improves the absolute vertical fidelity. Grumpe and Wöhler [1]^[7] generalized the variational approach and introduced a formalism to concurrently estimate the surface heights and the local reflectance properties in terms of the Hapke reflectance model, i.e., the albedo. Other recent approaches, which share many structural similarities, are, for example Wu et al. ^[9], Jiang et al. ^[18], and Alexandrov and Beyer ^[10]. Wu et al. ^[9] and Alexandrov and Beyer ^[10] use LRO-NAC images of the Moon and make use of the Lunar-Lambert model. The latter method additionally allows for multiple images. Jiang et al. ^[18] adopt the algorithm of Grumpe and Wöhler [1]^[7] and use the Mars-specific reflectance model and the Mars ReCo algorithm from Ceamanos et al. ^[28] to estimate the reflectance and atmospheric parameters.

4. Shape from Shading Applied to Mars

SfS has been applied to the Martian surface for forty years, and over time, different approaches of an increasing level of sophistication have been presented. This section provides a brief review of previous approaches for SfS applied to Mars and points out the methodical challenges and open questions addressed by this paper. Previous approaches primarily differ in the reflectance model, the specific implementation of the SfS algorithm and the dataset that determines the actual resolution. Early works applied photoclinometry to Mariner 9 images ^[24] and SfS to Viking imagery [12,13]^[12][13]. Dorrer et al. ^[29] and Dorrer et al. ^[14] employ SfS to refine stereo DEMs, which are derived from the High-Resolution Stereo Camera (HRSC) onboard the Mars Express orbiter ^[30]. O’Hara and Barnes ^[16] propose the Large Deformation Optimization Shape from Shading (LDO-SfS) technique for recovering the surface shape without initialization. All approaches use simple reflectance models such as Lambert or Oren–Nayar ^[31], and mostly assume a constant albedo, if any, and do not model any atmospheric effects (e.g., ^[16]). Gehrke ^[15] is the first to address atmospheric effects on Mars in the context of surface reconstruction. The work combines facet stereo vision with SfS, applied to HRSC imagery. The Lunar-Lambert model is used for radiometric modeling, and the optical depth of the atmosphere is estimated from two HRSC images acquired under different observation angles. Jiang et al. ^[18] were the first to propose a scheme for an integrated stereo and SfS approach on CTX images (up to 5 m/pix) with atmospheric compensation. They employ their previous reflectance model ^[28] to build a thorough physical reflectance and atmospheric model based on additional multi-angle CRISM measurements to estimate the model parameters. Compared to the photogrammetric reconstruction, the DEM results show improvement and are consistent with a sparse set of MOLA points. In a subsequent publication, Douté et al. ^[32] extend their approach to work with HiRISE imagery; due to the lack of validation data, the Isotropic Undecimated Wavelet Transform (IUWT) is employed for a consistency-check. Despite the overall viability and adequate results of SfS on Mars, wresearchers identify two methodological challenges and one open question related to atmospheric compensation. First, all previous works, and most recently the approach of Jiang et al. ^[18], assume constant albedo throughout the scene. If a more sophisticated area with locally varying albedo is examined, Jiang et al. ^[18] perform clustering and divide the scene into multiple regions with constant albedo. Even though algorithms for SfS with locally varying albedo exist (e.g., Wu et al. ^[9] and Grumpe and Wöhler [1]^[7]), they have not yet been fully explored for Mars. Secondly, modeling the atmosphere remains challenging and largely unexplored. Many studies have simply neglected atmospheric effects and Gehrke ^[15] only estimated the optical depth. The approach of Jiang et al. ^[18] requires external data to estimate the parameters of the reflectance model and the atmospheric model. Significant additional processing is required for obtaining the image parameters. Due to the lack of such data, Jiang et al. ^[18] assume globally averaged parameters. In general, it would be beneficial to directly estimate the atmospheric parameters from the image data such that no external information is necessary. Thirdly, Jiang et al. ^[18] employed an atmospheric model for their SfS implementation, but did not provide any further analysis of the influence of atmospheric conditions on the reconstruction results. In fact, they utilized images acquired under extremely clear atmospheric conditions with an optical depth of τ≈0.16 (G20_025904_2209_XN_40N102W), τ≈0.08 (B20_017600_1538_XN_26S183W), and τ≈0.19 (B21_017786_1746_XN_05S222W), as can be retrieved from the maps of Montabone et al. ^[33]. Hess et al. [7]^[6] combined an atmospheric model with the SfS procedure from Grumpe and Wöhler [1]^[7] and Grumpe et al. ^[8]. They provided a short comparison with and without an atmospheric model for the first scene from Jiang et al. ^[18] (G20_025904_2209_XN_40N102W). Even the model without atmospheric compensation produced good results.