# Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery

Cheng Jiang, Sylvain Douté, Bin Luo, Liangpei Zhang

► To cite this version:

Cheng Jiang, Sylvain Douté, Bin Luo, Liangpei Zhang. Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery. ISPRS Journal of Photogrammetry and Remote Sensing, 2017, 130, pp.418-430. 10.1016/j.isprsjprs.2017.06.010. insu-03692487

HAL Id: insu-03692487

https://insu.hal.science/insu-03692487v1

Submitted on 30 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery.

Cheng Jiang¹,†, Sylvain Douté*², Bing Luo³,⁴, Liangpei Zhang³,⁴

1 School of Geodesy and Geomatics, Wuhan University, Wuhan. China.

2 Université Grenoble Alpes/CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), UMR 5274, Grenoble, France.

3 State Key Laboratory of Information Engineering in Surveying Mapping and Remote Sensing, Wuhan University, Wuhan, China

4 Collaborative Innovation Center for Geospatial Information Technology, China.

*Corresponding author: Sylvain Douté, E-mail: sylvain.doute@univ-grenoble-alpes.fr
Phone: 33-476514171, Fax: 33-476448821.

**Abstract.** High-resolution Digital Elevation Models (DEMs) of the Martian surface are instrumental for studying the red planet. The Mars Orbiter Laser Altimeter (MOLA) instrument onboard the Mars Global Surveyor provided global DEM of high vertical resolution but with a limited spatial resolution that is not enough for characterizing small geological objects, normalizing illumination conditions across in-orbit optical images, modeling local meteorology, and other applications. This article addresses the problem of producing DEMs for regions of interest on Mars using available in-orbit imagery, typically ≈1000 km2 in area, while insuring a ≈10 meters vertical accuracy and a spatial accuracy which is comparable to that of the imagery. A method is proposed that combines photogrammetric and photoclinometric

<footnote>† The first author performed the work while at Université Grenoble Alpes/CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), Grenoble, France.</footnote>
approaches in order to retain their mutual advantages. According to experiments using
Mars Reconnaissance Orbiter Context Camera (CTX) images, the proposed method is
indeed able to produce DEMs satisfying the previous requirements, with less artifacts,
better surface continuity, and sharper details than the photogrammetric method when
it is used alone.

# 1. Introduction

Mars is the most similar Planet to Earth in the Solar System and offers
opportunities for fruitful comparative planetary studies as well as human exploration
and settlement. For a brief introduction to the planet knowledge, we refer readers to
Barlow (2008). Since the 60s of the last century, we have had over 40 space missions
sent to Mars allowing the acquisition of a wealth of data including extended
collections of images. To perform studies on Mars, we have to know in great details
the shape of the surface, namely the topographic information (Kirk et al., 2008). The
latter is not only mandatory for geology investigations, meteorology simulations, and
orthoimage generation, but is also critical for improving the precision of surface
reflectance retrieval, for spectral unmixing of hyperspectral images (Ceamanos et al.,
2011; Luo et al., 2013) and so on.

The topographic information can be represented by triangulated irregular
networks (TINs), which are a collection of vectorial objects, or digital elevation
models (DEMs) which are raster images. Since we can transform between these two
representations, and the raster form is more suitable for optimization problems in this
paper involving the gradient and orientation of surfaces, we will focus on DEMs. Two
types of techniques have been operated so far to produce DEMs on Mars. One is using
range measurements produced by instruments such as Lidar (light detection and
ranging) and Radar (radio detection and ranging). The Mars Orbiter Laser Altimeter
(MOLA) instrument onboard the Mars Global Surveyor operated for seven years until
2006 enabling the production of a global DEM of very high vertical resolution (less
than 1 meter), but of limited spatial resolution (nominally 168 meters/pixel) (Smith et
al., 2001). The individual data points are accessible to the public (PDS, 2009) and
their estimated absolute heights can be used as a reference for validating/calibrating
other methods.

The second type of techniques addresses the production of regional to local
DEMs (≈10⁵-10² km² in area) based on images. In the literature, there are two
approaches, namely photogrammetry (or stereo) and photoclinometry (or shape from
shading, SFS). Photogrammetric methods need at least two images capturing the same
region of interest but from a different viewpoint. It is based on the effect of parallax
that can be exploited by triangulation. The core part of photogrammetry is matching
pair of pixels, each belonging to one image, and solving a system of collinearity
equations involving the image coordinates of the matched pixels as observations, the
3D coordinates of the corresponding point on the scene, and camera parameters. We
refer readers to McGlone et al. (2013). Photoclinometric methods, can derive surface
gradients (i.e. slopes) from the intensity variation of an image, therefore allow to
produce a DEM from a single image if the surface can be reconstructed from these
gradients by integration for example. Photoclinometry can be considered as an inverse
problem of the image formation process with an emphasis on brightness. Good review
works are given by Horn (1990) and Zhang et al. (1999). According to the discussions
presented in Leclerc and Bobick (1991), the two approaches have the following own
advantages and drawbacks. Photogrammetric methods are fast, able to directly
produce DEMs with absolute heights – even though a calibration upon a reference is
often required - and have good comparative performances for large scale objects of
interest. However, the matching of pixels needs recognizable local intensity patterns
on both pair images that, consequently, need to be acquired with comparable surface
and illumination conditions. Furthermore the discriminability of height differences
decreases as the distance between the sensor and the scene increases. As for
photoclinometric methods, the discriminability does not depend on distance, and no
matching of pixels is needed, thus preserving small details in the DEM. However, the
method is relatively complex and slow to converge, requires integration of gradient
fields, and cannot produce absolute heights. Furthermore the modeling of intensity
spatial variations in the image entails a priori information on the bidirectional
reflectance properties of the surface or estimation of the latter from the images
themselves.

To benefit from as many advantages as possible, some pioneer works have
contributed to developing new methods by fusing photogrammetry and
photoclinometry. Such methods can be classified into three categories.

The methods of the first category directly combine the results by
photogrammetry and photoclinometry. In this category, Cryer et al. (1993) performed
an addition of components derived from spatial low pass filtering of the
photogrammetric result and the components derived from spatial high pass filtering of
the photoclinometric result; Samaras et al. (2000) first applied photogrammetric
method to the whole image and then selected some areas to perform photoclinometric
method independently. Kumar et al. (2009) adopted a complex strategy which builds
the relationship between the results of the two methods using a neural network. Where
photogrammetry has advantages its results are retained, whereas for areas lacking
textures, the photoclinometric results are transformed to be photogrammetry-like.

The methods of the second category successively apply photogrammetry and
photoclinometry to the same data. For this category, the height field derived by
photogrammetry often serves as a constraint in the photoclinometric model (Dorrer et
al., 2004; Dorrer and Zhou, 1998; Fassold et al., 2004; Grumpe et al., 2014; Pien and
Gauch, 1995).

For the last category, photogrammetry and photoclinometry are combined into a
single model. The representative work is done by Fua and Leclerc (1995), Lohse et al.
(2006)_ENREF_18, and Gehrke (2006). In their work, the solution is found in an object
space where the height of the point would affect both the position and intensity value
of the pixel in every observed image, but the numerical stability of the algorithms are
only insured if the images were taken according to oblique views. The three
categories of methods follow an increasing fusion level, as well as increasing
complexity.

Based on the pioneer works, we propose a new fusion algorithm that belongs to
the second category. Our choice is motivated by the existence of an independent,
well-known and capable photogrammetric tool, which is open to public and keeps
updating, i.e. NASA Ames Stereo Pipeline (ASP) (NASA, 2009; Shean et al., 2016).

In the proposed algorithm, the photoclinometric scheme subsequently takes this
"original" DEM as an input along with the image and refines it at small scales by
inverting a radiative transfer model that depicts the image intensity field for a given
topography, atmospheric, illumination and viewing conditions. Another specific
aspect of the algorithm is that the model is inverted with an optimization procedure
containing two regularization terms. The first one insures that variables describing the
height field on the one hand and the gradients on the other hand are separated for a
better numerical stability. The second regularization term insures that the algorithm
retains the properties of the original DEM at large scales. It also speeds up the
convergence and constrains the solution space. As a result, the complexity of the
optimization scheme is moderate which implies flexibility if one needs further
adjustment. Furthermore our strategy avoids considering complex rules for fusion,
and the method does not need oblique views, i.e. is more suitable for images produced
by most Martian sensors. Finally, the novelty of the proposed method is multiple. First
the intensity model integrates an innovative radiative transfer scheme that takes into
account atmospheric effects such as the diffuse scattered radiance from the sky. In
current Martian conditions, the latter constitutes a significant source of illumination
modulated locally depending on the fraction of sky visible in addition to the direct
illumination by the sun. Second the model also uses a realistic description of the
bidirectional reflectance distribution function (BRDF) of the surface, as an anisotropic semi-empirical kernel-based model, namely the Ross-Thick Li-Sparse (RTLS) model. The kernel weights of the RTLS model for a given scene are retrieved using the Mars-Reco algorithm (Ceamanos et al., 2013) applied to multi-angular sequences of hyperspectral images by the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) if they are available. Finally, this work is the first successive application of photogrammetry and photoclinometry applied to CTX (Context Camera) images (ASU, 2006). All the data treated by the experiments is downloaded through the Mars Orbital Data Explorer (http://ode.rsl.wustl.edu/mars/index.aspx; last accessed 7 June 2017).

The article is organized as follows. The proposed model and the optimization scheme are detailed in Section 2. In Section 3 the method is tested with increasingly complex data. Eventually, the advantages and problems of the proposed method are discussed in Section 4.

## 2. Methods

The proposed method exploits photogrammetric and photoclinometric information in a successive way, as shown in Fig. 1.
<figure><img src="image_2.png" /><figcaption>Fig. 1 The flowchart of the proposed method</figcaption></figure>

## 2.1 DEMs derived by photogrammetric method

The DEMs derived by photogrammetric method are used as initial DEMs in the iterative procedure during optimization and also for setting one of the regularization terms. For the photogrammetric method, we adopt version 2.5.1 of ASP, which comes with a very detailed guide book available online (NASA, 2009). For generating the initial DEM, we follow the procedure recommended by Öhman (2013). In this paper, the procedure mainly follows three steps. First, the CTX images are selected and downloaded (ASU, 2006). Then the images are preprocessed, such as applying a radiometric calibration and map projection, using the Integrated Software for Imagers and Spectrometers (ISIS) (USGS, 1992). Finally, ASP is used to derive the DEM. The default parameters of ASP as defined by Öhman (2013) are used for all the experiments in this paper.
## 2.2 Expression for the bottom of the atmosphere reflectance

In order to implement the photoclinometric scheme, the imaging process must be modeled. For simplicity we assume a linear relationship between image intensity and the bottom-of-atmosphere (BOA) reflectance implying that the image must first be corrected for the atmospheric path radiance – an additive term corresponding to the photons scattered by the atmosphere alone towards the sensor and atmospheric attenuation between the surface and the sensor (Ceamanos et al., 2013). According to Eq. 2 in Ceamanos et al. (2013), the BOA reflectance factor can be expressed as:

$$ \mathbf{R}_{BOA} = S\mu_0 e^{-\tau_0/\mu_0}\rho(s_0, s, p, q) + S\mu_0 e^{-\tau_0/\mu_0}\alpha c_0 \rho_1(\mu)\rho_2(\mu_0) \\ + \frac{\alpha}{\pi} \int_{\Omega^+} D_s(s_0, s')\rho(s', s, p, q)\mu' ds' \quad (1) $$

where S is the extraterrestrial solar spectral irradiance, $s_0$ and S are the sun and view directions defined by $(\theta_0, \varphi_0)$ and $(\theta, \varphi)$ respectively ($\theta_0$ and $\theta$ are the corresponding zenith angles, and $\varphi_0$ and $\varphi$ are the corresponding azimuth angles), $\mu_0$ and $\mu$ are the cosine of solar and view zenith angle respectively, $ds'$ is an integration variable representing all possible incoming directions of light from the sky, $\tau_0$ is the atmospheric opacity, $c_0$ is the spherical albedo of the atmosphere, $\alpha$ is the multiple reflection factor that depends on the surface albedo $q(\mu_0)$ such that $\alpha = (1-q(\mu_0)c_0)^{-1}$, $\rho = \pi f$ (where f is the BRDF), p and q denote the gradient of the surface height according to the horizontal and vertical axis of the image, and $D_s(s_0, s)$ is the path radiance incident on the surface. $\rho_1(s)$ and $\rho_2(s_0)$ are defined as:
$$ \rho_1(s) = \frac{1}{2\pi} \int_{\Omega^+} \rho(s', s, p, q) ds', \rho_2(s_0) = \frac{1}{2\pi} \int_{\Omega^-} \rho(s_0, s', p, q) ds' \quad (2) $$

Here, the sign (+) and the sign (-) represent downward and upward directions, respectively. Eq. (1) consists of three terms that represent the reflection of the anisotropic sources (sun and its aureole), the multiple scattering between surface and atmosphere and the reflection of the isotropic component of the sky illumination, respectively. Note that this formulation is adapted to integrate the parameterization of (Spiga and Forget, 2008) for the fast and accurate estimation of solar and sky irradiance for any Martian slope. In addition, in this chapter we omit the dependency on wavelength in the equations. For CTX images, we used the central band of 0.6µm for all the experiments.

When adopting the RTLS model, the reflectance factor $\rho$ can be expressed as:

$$ \rho(s_0, s, p, q) = k^L + k^G f_G(s_0, s, p, q) + k^V f_V(s_0, s, p, q) \quad (3) $$

where the three terms are considered as Lambertian (L), geometric (G) and volumetric (V) components, respectively. The function $f_G$ and $f_V$ are predefined generic kernels, and $k^L$, $k^G$ and $k^V$ are their weights which are defined by the surface reflective properties. Since the components of $f_G$ and $f_V$ are usually relatively smaller than the first term, we can assume $p=0$ and $q=0$ in Eq. (2) implying a flat surface for the surface-atmosphere multiple reflections, to make the equation simpler. For efficiency purposes, besides the inputs of $\tau_0$, $s_0$, $s$, $p$, $q$, $k^L$, $k^G$ and $k^V$, all the other components in Eq. (1) can be retrieved or calculated using a look-up table which is pre-computed using the radiative transfer program named Discrete Ordinates Radiative Transfer Program for a Multi-layered Plane-Parallel
Medium (DISORT) (Stamnes et al., 1988). Since this paper focuses on getting high-quality DEMs, interested readers are referred to Ceamanos et al. (2013) for more details.

## 2.3 The photoclinometric model

The core part of the photoclinometric methods is the brightness constraint or the so-called “image irradiance equation” (Horn, 1990; Zhang et al., 1999) which can be expressed as:

$$ \beta \cdot \mathbf{R} = \mathbf{I} \qquad (4) $$

where **I** is the measured image intensity, **R** is the modeled reflectance map, and $\beta$ is a scalar factor. For Mars, we choose to express **R** using $\mathbf{R}_{BOA}$ in Eq. (1). For a certain image, **R** is usually modeled as a function of surface gradient vectors (lexicographically stacked) **p** and **q** by omitting the other variables. Therefore, by solving Eq. (4), we can derive the surface gradients, and this is the basic idea of photoclinometric methods. However, the number of unknowns is twice the number of equations, so extra constraints are needed to narrow the solution space. Such as the smoothness constraint, the integrability constraint, the gradient constraint, the unit normal constraint and so on (Zhang et al., 1999).

In the literature dedicated to photoclinometric methods, the height information $z$ or specifically the DEM can be acquired by integrating the solution **p** and **q** (Grumpe et al., 2014; Grumpe and Wöhler, 2014) (two step strategy) or by including the height field $z$ in the optimization scheme (Fassold et al., 2004; Leclerc and Bobick, 1991) (one step strategy). To be able to produce DEMs from images directly, while
avoiding the accumulation of errors, the one step strategy is employed in our method. Nevertheless **p** and **q** are still used in our model in the form of independent variables, so as to prevent the optimization problem from being a complicated variational problem. This consideration is accompanied by adopting the integrability constraint. Moreover, the initial DEM **z**₀ derived by the ASP photogrammetric method (Section 2.1) is used to constrain the solution by retaining its properties at large scales. So, the proposed model that satisfies all these requirements can be expressed as:

$$ 
\begin{gathered}
\mathbf{z} = \operatorname*{arg\,min}_{\mathbf{p,q,z}, \beta} \mathbf{E} \\
\text{with } \mathbf{E} = \frac{1}{2} (\boldsymbol{\beta} \cdot \mathbf{R}(\mathbf{p}, \mathbf{q}) - \mathbf{I})^2 + \frac{\lambda_1}{2} [(\mathbf{p} - \mathbf{z}_x)^2 + (\mathbf{q} - \mathbf{z}_y)^2] + \frac{\lambda_2}{2} (\mathbf{W} \cdot \mathbf{z} - \mathbf{W} \cdot \mathbf{z}_0)^2
\end{gathered}
\quad (5) $$

where **W** is a mask to filter out the points where ASP fails to find the matched points, $\lambda_1$ and $\lambda_2$ are regularization parameters to balance the contributions of the different terms, while the third term represents the photogrammetric constraint. $\mathbf{z}_x$ and $\mathbf{z}_y$ are the calculated gradients derived by convolution of **z** with filters $f_x$ and $f_y$ respectively:

$$ f_x = \frac{1}{8c} \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{bmatrix}, f_y = \frac{1}{8c} \begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{bmatrix} \quad (6) $$

where c is the cell size. Invoking the same reason as in Grumpe et al. (2014), we do not consider the empirical term that enforces the smoothness of the surface for better results. Note that, unlike the work by Grumpe et al. (2014) where $\frac{\partial \mathbf{z}_0}{\partial x}$ and $\frac{\partial \mathbf{z}_0}{\partial y}$ are used to construct the photogrammetric constraint, we use $\mathbf{z}_0$ directly, in order to not
amplify the noise contained in. $\mathbf{z}_0$.

## 2.4 Numerical solutions

Since Eq. (5) is an energy function to be minimized with more than one variable, the most direct way is to alternately minimize the function with each variable successively.

For $\mathbf{p}$ and $\mathbf{q}$, the updating equations are obtained by using classic gradient descent:

$$ \mathbf{p}^{n+1} = \mathbf{p}^n - \alpha_1^n \frac{\partial E}{\partial \mathbf{p}} = \mathbf{p}^n - \alpha_1^n \left[ \left( \beta^n \mathbf{R}(\mathbf{p}^n, \mathbf{q}^n) - \mathbf{I} \right) \beta^n \frac{\partial \mathbf{R}(\mathbf{p}^n, \mathbf{q}^n)}{\partial \mathbf{p}} + \lambda_1 (\mathbf{p}^n - \mathbf{z}_x^n) \right] \quad (7) $$

$$ \mathbf{q}^{n+1} = \mathbf{q}^n - \alpha_2^n \frac{\partial E}{\partial \mathbf{q}} = \mathbf{q}^n - \alpha_2^n \left[ \left( \beta^n \mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^n) - \mathbf{I} \right) \beta^n \frac{\partial \mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^n)}{\partial \mathbf{q}} + \lambda_1 (\mathbf{q}^n - \mathbf{z}_y^n) \right] \quad (8) $$

where $\alpha_1^n$ and $\alpha_2^n$ are the step sizes that can be decided according to line-search algorithms. In this work, the inexact line search algorithm proposed by Fletcher (2013) is adopted. Since $\mathbf{R}$ is a very complicated function, the partial derivatives of $\mathbf{R}$ with respect to $\mathbf{p}$ and $\mathbf{q}$ are calculated using numerical approximations:

$$ \frac{\partial \mathbf{R}(\mathbf{p}^n, \mathbf{q}^n)}{\partial \mathbf{p}} = \frac{\mathbf{R}(\mathbf{p}^n + \varepsilon, \mathbf{q}^n) - \mathbf{R}(\mathbf{p}^n, \mathbf{q}^n)}{\varepsilon} \quad (9) $$

$$ \frac{\partial \mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^n)}{\partial \mathbf{q}} = \frac{\mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^n + \varepsilon) - \mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^n)}{\varepsilon} \quad (10) $$

where $\varepsilon$ is a small constant, and we set it to 0.0001 in this paper.

For $\mathbf{z}$, we deal with a simple variational problem of which solution can be found by solving the classical Euler-Lagrange equation:
$$ \lambda_1 \left( \Delta \mathbf{z} - \frac{\partial \mathbf{p}}{\partial x} - \frac{\partial \mathbf{q}}{\partial y} \right) = \lambda_2 \mathbf{W} (\mathbf{z} - \mathbf{z}_0) \quad (11) $$

that can be optimized using Gauss-Seidel method. Similar to the method in Jiang et al.

(2016), the numerical updating equation is:

$$ \mathbf{z}_{i,j}^{n+1} = \frac{[c^2 \mathbf{W}_{i,j} \lambda_2 (\mathbf{z}_0)_{i,j} / \lambda_1 + c (\mathbf{p}_{i-1,j}^{n+1} - \mathbf{p}_{i,j}^{n+1} + \mathbf{q}_{i+1,j}^{n+1} - \mathbf{q}_{i,j}^{n+1}) + \mathbf{z}_{i-1,j}^{n} + \mathbf{z}_{i+1,j}^{n} + \mathbf{z}_{i,j-1}^{n} + \mathbf{z}_{i,j+1}^{n}]}{(c^2 \mathbf{W}_{i,j} \lambda_2 / \lambda_1 + 4)} \quad (12) $$

where the subscripts i and j denote the position of the pixel. Then $\mathbf{z}_x^{n+1}$ and $\mathbf{z}_y^{n+1}$ are also be updated by convolving $\mathbf{z}^{n+1}$ with filters defined in Eq. (6).

For $\beta$, by setting the derivative of $\mathbf{E}$ with respect to $\beta$ to 0, the closed form of the solution is expressed as:

$$ \beta^{n+1} = \frac{(\mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^{n+1}))^T \mathbf{I}}{(\mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^{n+1}))^T (\mathbf{R}(\mathbf{p}^{n+1}, \mathbf{q}^{n+1})))} \quad (13) $$

Finally, the iteration number $n$ is increased by 1, and a new iteration is started.

To sum up, the algorithm aimed at finding the solution of Eq. (5) consists of updating $\mathbf{p}$, $\mathbf{q}$, $\mathbf{z}$ and $\beta$ successively using Eq. (7), (8), (12), and (13) until the stop criteria is satisfied. As for the stop criteria, we adopt a similar idea as in Grumpe and Wöhler (2014): the algorithm will stop when the maximum iteration number $N_{iter}$ is reached, or a smaller value of the energy function $E$ is no longer found for $N_{step}$ iterations since the last successful decrease.

### 3. Experiments and discussion

To test the proposed method, we have performed experiments with different CTX datasets of increasing complexity. By complexity we mean the fraction of missing values and the density of artifacts in the original DEM produced by ASP as
well as the degree of spatial inhomogeneity and contrast affecting the intrinsic albedo
of the surface. The intrinsic albedo is defined as the fraction of incident irradiance
reflected by materials in any direction. The test datasets are mostly composed of CTX
images that have been cropped according to regions of interest 800×800 pixels
wide (Fig. 2, first row). Nevertheless the fact that our algorithm is capable of
processing larger regions (5000x5000 pixels wide corresponding to ≈1000km2 in area)
will be illustrated in the first experiment. For the subsequent experiments the reasons
for presenting regions of limited extent are as follows: first, the results are allowed to
be presented in greater details for giving a better description of the proposed method's
characteristics. Second, the solution for every pixel of the proposed method only
relies on a limited neighborhood. This means that the size of the refined DEM has
limited influence on the quality of the reconstruction. In addition we have verified that
the meta-parameter values of the optimization scheme stay relevant for a large area
regarding possible changes of topographical, geometrical and physical conditions.
Applying the method to large areas is mainly a matter of computing time, of the order
of a few hours for a ≈1000km2 area.

The second row of Fig. 2 shows the shaded images of the initial DEMs acquired
by using ASP (The sun position is the same as in the corresponding images). The three
datasets imply an increased complexity. Because of an extensive mantling by dust
deposited from the atmosphere, the intrinsic albedo of the surface is known to be
fairly homogenous for the first two datasets contrary to the third one. As for the initial
DEMs, the first one is relatively good with no missing values (i.e. holes) and not too
many artifacts. The second one has a high density of artifacts that are featured as abrupt changes and discontinuities due to irregularities with sensor scanning and due to intrinsic limitation of photogrammetric approach. Nevertheless it has only a few small holes. The last one has some big holes meaning no initial height information is available in those areas. For efficiency, only the left image was used in model of eq. (5). For the stopping criteria, $N_{iter}$ was set to 100000, and $N_{step}$ was set to 300.

<figure><img src="image_3.png" /><figcaption>Fig. 2. The tested CTX cropped images (first row) and the initial DEMs acquired by ASP (second row).</figcaption></figure>

All the results have been analyzed qualitatively and quantitatively. For visual inspections, the shaded images of the DEMs will be shown as well as the simulated BOA reflectance image since their intensity variations are very sensitive to the topographic solution. For quantitative assessment, the DEMs were compared with the available MOLA points, and the Root Mean Square Error (RMSE) between both
products were calculated to evaluate the height accuracy. In addition, the original left image was compared with the simulated BOA reflectance image (defined in Eq. (1)) by means of the Gradient Magnitude Similarity Deviation (GMSD) index (Xue et al., 2014) which constitutes an additional indicator to evaluate the accuracy of the surface derived from the DEM. The GMSD index was invented to evaluate the similarity of two images in the gradient domain, the standard deviation of the Gradient magnitude similarity being defined as:

$$GMS(i) = \frac{2\mathbf{m}_r(i)\mathbf{m}_d(i) + \xi}{\mathbf{m}_r^2(i) + \mathbf{m}_d^2(i) + \xi} \quad (14)$$

where *i* denotes the pixel index, **m**<sub>r</sub> and **m**<sub>d</sub> are the gradient magnitude images which are acquired by convolving the images with Prewitt filters, and ξ is a positive constant for numerical stability. Smaller RMSE and GMSD imply better results.

The atmosphere opacity τ<sub>0</sub> in Eq. (1) was assumed to be at usual levels, i.e. 0.5 on Mars. The regularization parameters λ<sub>1</sub> and λ<sub>2</sub> in Eq. (5) were empirically chosen from the series of values [0.1,0.01,0.001,0.0001] and [1.e-5, 1.e-6,1.e-7,1.e-8,1.e-9,1.e-10,1.e-11], respectively.

## 3.1 Results of the first CTX image pair.

For the first experiment, the information of the two images is listed in Table 1. Since there is no Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) multi-angular acquisition in this area, the kernel weights of the RTLS model were set using the globally averaged values for mineral dusty terrains on Mars, i.e. k<sup>L</sup> = 0.3966, k<sup>G</sup> = 0.0330, k<sup>V</sup> = 0.1277 (Vincendon, 2013). The λ<sub>1</sub> and λ<sub>2</sub> were set to 0.01 and 1e-6 respectively, and there are 88 MOLA points available for the
whole area. First Fig. 3 illustrates the capability of our algorithm to process large CTX
left and right images 5000x5000 pixels wide corresponding to ≈1000km2 in area.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Left: screenshot of the 3D visualization of the DEM produced by the proposed method (the vertical exaggeration factor is x2). Right: corresponding shaded image. In both cases, the sun azimuth is 256.63° and the sun elevation is 50.20°, same as for the left component of the CTX image pair.</figcaption></figure>

Then the results obtained for the cropped area respectively using only ASP and
the proposed method are summarized in Fig. 4. Compared with the result using ASP,
the result of the proposed method has less artifacts, its surface become more
continuous, its edges become sharper and more details are revealed. The result of the
proposed method also has smaller RMSE and GMSD that stresses the consistency of
its good performance.

<table><caption>Table 1
Image information of the first CTX pair.</caption><thead><tr><th></th><th>Left image</th><th>Right image</th></tr></thead><tbody><tr><td>Image ID</td><td>G20_025904_2209_XN_40N102W</td><td>G20_025970_2217_XN_41N102W</td></tr><tr><td>Sun elevation</td><td>50.20°</td><td>48.74°</td></tr><tr><td>Sensor elevation</td><td>82.53°</td><td>80.10°</td></tr><tr><td>Sun azimuth angle</td><td>256.63°</td><td>258.06°</td></tr><tr><td>Phase angle</td><td>47.25°</td><td>31.45°</td></tr><tr><td>Used Region</td><td colspan="2">Upper-left corner, lat: 42.9316°, lon: -102.8187°;<br>Lower-right corner, lat: 42.8572°, lon: -102.7443°.</td></tr></tbody></table>
<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>(a) ASP
RMSE=7.2342, GMSD=0.1653</figcaption></figure>

<caption>(b) The proposed method
RMSE=7.0928, GMSD=0.0621</caption>

<caption>Fig. 4. Results obtained using ASP only (first column) and the proposed method (second column). The shaded images in the first row are produced when the sun azimuth is 256.63° and the sun elevation is 50.20°, same as for the left component of the CTX image pair. For the</caption>
shaded images in the second line, the sun azimuth is 180° and the sun elevation is 45°. In the third line, the BOA reflectance image (eq. 1) simulated using the corresponding DEMs is presented and ought to be compared with image (row1, column1) of Fig. 2. Note also the respective values of the RMSE and GMSD indicators.

## 3.2 Results of the second CTX image pair

For the second experiment, the information of the two images is listed in Table 2. The kernel weights of the RTLS model were supplied by applying Mars-Reco on the CRISM acquisition sequence FRT9E0D and by averaging the obtained surface properties over the whole region of interest leading to $k^L = 0.1028, k^G = 0.0042, k^V = 0.0627$. The $\lambda_1$ and $\lambda_2$ were set to 0.001 and 1e-9 respectively, and there are 25 MOLA points available in the selected area. The results are shown in Fig. 5. By using the proposed method better results were also obtained. However, the shaded image of the refined DEM when the sun position is very different to that of the original image reveals that artifacts of a different type have been introduced by the photoclinometric scheme. They appeared as a wavy pattern of weak amplitude whose main direction is aligned with the original direction of the sun for the left component of the CTX image pair. Similar observations have also been reported in the work of O’Hara and Barnes (2012). The reason will be discussed in Section 3.6.

<table><caption>Table 2
Image information of the second CTX pair.</caption><thead><tr><th></th><th>Left image</th><th>Right image</th></tr></thead><tbody><tr><td>Image ID</td><td>B20_017600_1538_XN_26S183W</td><td>B21_017745_1538_XN_26S183W</td></tr><tr><td>Sun elevation</td><td>19.57°</td><td>19.95°</td></tr><tr><td>Sensor elevation</td><td>89.91°</td><td>75.52°</td></tr><tr><td>Sun azimuth angle</td><td>312.70°</td><td>313.32°</td></tr><tr><td>Phase angle</td><td>70.36°</td><td>79.72°</td></tr><tr><td>Used Region</td><td colspan="2">Upper-left corner, lat: -25.8450°, lon: 176.7125°; Lower-right corner, lat: -25.9100°, lon: 176.7774°.</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Results obtained using ASP only (first column) and the proposed method (second column). The shaded images in the first row are produced when the sun azimuth is 312.70° and</figcaption></figure>
the sun elevation is 19.57°, same as for the left component of the CTX image pair. For the shaded images in the second line, the sun azimuth is 180° and the sun elevation is 45°. In the third line, the BOA reflectance image (eq. 1) simulated using the corresponding DEMs is presented and ought to be compared with image (row1, column2) of Fig. 2. Note also the respective values of the RMSE and GMSD indicators.

We now analyze the impact of the RTLS kernel weights supplied by Mars-Reco algorithm on the results. For that purpose we also present in Fig. 6 the best result obtained using the globally averaged kernel weights mentioned in Section 3.1 and from varying the regularization parameters, $\lambda_1$ and $\lambda_2$ being set to 0.01 and 1e-10 respectively. The result was compared with the result derived from the Mars-Reco supplied kernel weights. We can see that adopting the Mars-Reco kernel weights leads to better results with sharper details and significantly lower RMSE and GMSD. The improvement (going from the right part of Fig. 6 to the left part) is significant since it represents respectively 37% (13%) of the RMSE (GMSD) difference between the initial DEM and the refined DEM (going from the left part of Fig. 5 to the right part). Therefore, the performance of the proposed method depends on the accuracy of the initial DEM as well as the accuracy of the RTLS kernel weights.

<figure><img src="image_8.png" /><figcaption>(a) Using the kernel weights supplied by (b) Using the kernel weights of globally</figcaption></figure>
<table><tr><td>Mars-Reco</td><td>averaged values for mineral dusty terrains on Mars</td></tr><tr><td>RMSE=5.7416, GMSD=0.0457</td><td>RMSE=5.8841, GMSD=0.0769</td></tr></table>

Fig. 6. Comparison of the shaded images generated by the proposed method using two different sets of RTLS kernel weights. The sun azimuth is 312.70° and the sun elevation is 19.57° according to the left component of the CTX image pair.

## 3.3 Results of the third CTX image pair

For the third experiment, the information of the two images is listed in Table 3. Since the test area contains at least three distinct classes of intrinsic surface albedo, the left CTX image was first segmented using the mean shift algorithm (Comaniciu and Meer, 2002), as shown in Fig. 7. Then, the RTLS kernel weights supplied by applying Mars-Reco on the CRISM acquisition sequence FRT17327 were averaged accordingly and assigned to each class. If we sort the classes in the order they appear in the CTX image from left to right, the weights can be presented in the form of vectors with each component corresponding to one class:

$$k^L = [0.0749, 0.1105, 0.1568]$$
$$k^G = [-0.0244, -0.0186, -0.0004]$$
$$k^V = [0.1605, 0.1579, 0.1433]$$

The parameters $\lambda_1$ and $\lambda_2$ were set to 0.01 and 1e-8 respectively, and there are only 16 MOLA points available in the selected area, specifically in the area of the first class.

<table><caption>Table 3
Image information of the third CTX pair.</caption><thead><tr><th></th><th>Left image</th><th>Right image</th></tr></thead><tbody><tr><td>Image ID</td><td>B21_017786_1746_XN_05s222w</td><td>B21_017931_1746_XN_05s222w</td></tr><tr><td>Sun elevation</td><td>32.04°</td><td>32.36°</td></tr><tr><td>Sensor elevation</td><td>75.12°</td><td>80.10°</td></tr><tr><td>Sun azimuth angle</td><td>304.64°</td><td>305.23°</td></tr><tr><td>Phase angle</td><td>47.38°</td><td>58.48°</td></tr><tr><td>Used Region</td><td colspan="2">Upper-left corner, lat: -5.0812°, lon: 137.0650°; Lower-right corner, lat: -5.1503°, lon: 137.1341°.</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Segmenting the original left image using the mean shift algorithm.</figcaption></figure>

The results are shown in Fig. 8. As we can see from the latter, using the proposed method rather than simply applying the photogrammetric method allows filling the holes, avoiding discontinuities and getting sharper details. However, it comes at the expense of introducing wavy artifacts in the DEM, which are obvious in the shaded image when the sun direction is very different from the one of the original image. The RMSE value of the result of the proposed method is worse than that of ASP although such comparison may be biased because there are only 16 MOLA points available for the evaluation of a DEM of size 800×800. This is probably not statistically meaningful.
<figure><img src="image_10.png" /><figcaption>Fig. 8. Results obtained using ASP only (first column) and the proposed method (second column). The shaded images in the first row are produced when the sun azimuth is 304.64° and the sun elevation is 32.04°, same as for the left component of the CTX image pair. For the</figcaption></figure>
shaded images in the second line, the sun azimuth is 180° and the sun elevation is 45°. In the third line, the BOA reflectance image (eq. 1) simulated using the corresponding DEMs is presented and ought to be compared with image (row1, column3) of Fig. 2. Note also the respective values of the RMSE and GMSD indicators.

### 3.4 The influence of the meta-parameters

The proposed method contains two meta-parameters $\lambda_1$ and $\lambda_2$ to control the relative contribution of the different regularization terms. Take the first dataset for example, the quantitative results using different combinations of these parameters are shown in Table 4. The best combination that insures the minimum RMSE and GMSD is marked bold and the second best combination is marked with an underline. As we can see from the results by correlation, parameter $\lambda_1$ mainly influences the GMSD index while parameter $\lambda_2$ mainly influences the RMSE index. This phenomenon is not surprising. For smaller $\lambda_1$, the contribution of the image intensity constraint (first term of Eq. 5) becomes more dominant, therefore better GMSD can be expected as less differences between the simulated BOA image and the original image are achieved. For parameter $\lambda_2$, values exceeding a certain threshold will cause the result being too close to the noisy DEM obtained by ASP, and too small values $\lambda_2$ will cause the lost of the absolute height information.

<table><caption>Table 4
Influence of the regularization parameters on the first dataset.</caption><thead><tr><th></th><th>λ<sub>1</sub></th><th>0.1</th><th>0.1</th><th>0.1</th><th>0.1</th><th>0.01</th><th>0.01</th><th>0.01</th><th>0.01</th><th>0.001</th><th>0.001</th><th>0.001</th><th>0.001</th></tr><tr><th></th><th>λ<sub>2</sub></th><th>1e-5</th><th>1e-6</th><th>1e-7</th><th>1e-8</th><th>1e-5</th><th>1e-6</th><th>1e-7</th><th>1e-8</th><th>1e-5</th><th>1e-6</th><th>1e-7</th><th>1e-8</th></tr></thead><tbody><tr><th>RMSE</th><td>7.1534</td><td>7.1097</td><td><strong>7.0910</strong></td><td>7.0967</td><td>7.1753</td><td><u>7.0928</u></td><td>7.3096</td><td>7.8577</td><td>7.2225</td><td>7.1618</td><td>7.1472</td><td>8.0543</td></tr><tr><th>GMSD</th><td>0.1079</td><td>0.1032</td><td>0.1041</td><td>0.1044</td><td>0.0793</td><td><u>0.0621</u></td><td>0.0618</td><td>0.0624</td><td>0.1289</td><td>0.0623</td><td><strong>0.0494</strong></td><td><strong>0.0494</strong></td></tr></tbody></table>
### 3.5 The convergence property and the computation time

Still considering the first dataset, the convergence speed of our algorithm is shown in Fig. 9. As we can see from the curve, the optimization algorithm described in Section 2.4 can effectively reduce the energy of the proposed method, and the convergence rate is satisfactory. Similar results are found for other datasets and are omitted to save space.

<figure><img src="image_11.png" /><figcaption>Fig. 9. Reduction of the energy on the first dataset during optimization.</figcaption></figure>

The computation time for each iteration with respect to varying number of pixels (i.e. image size) is evaluated from experiments on the first CTX dataset and is shown in Fig. 10. Note that the tests were performed using Matlab with Intel CPU of 2.60 GHz. The figure implies that the complexity of the algorithm is approximately linear with the image size.
<figure><img src="image_12.png" /><figcaption>Fig. 10. Computation time cost for each iteration with respect to varying image size</figcaption></figure>

3.6 Analysis of the remaining artifacts

As can be seen from the results in Fig. 5 (row 2, column 2) and Fig. 8 (row 2, column 2), the photoclinometric scheme introduces in some cases artifacts in the DEM in the form of “plane waves” of limited amplitude affecting the field z with a “direction of propagation” approximately orthogonal to the original direction of the sun in the left component of the CTX image pair. It is worth noting that, contrary to experiment 1, the incidence angle of experiments 2 and 3 is large (≈70°). Fig. 11 illustrates the fact that this phenomenon is likely due to a relative invariance of the intensity **R**(*p*,*q*) for some specific variations of *p* and *q*. When the orientation of a planar facet of terrain is changed such that its normal sweeps out a cone centered at the sun direction, the local incidence angle (which controls the angular aperture of the cone) is invariant by this transformation. The phase angle also remains unchanged. In some cases, for example if the surface behaves like a Lambertian reflector, the local
incidence angle is the most influent parameter that controls the intensity of the
corresponding pixel through the first term of the photometric model (Eq. (1))
regardless of the local emergence angle and azimuth. Then the intensity will be
invariant too. As a consequence, the brightness constraint alone in the
photoclinometric scheme of Eq. (5) sometimes entails an indetermination for any
"precession" of the facet normal around the sun direction. This "precession" should
be limited since the surface continuity needs to be preserved. Besides, in the studied
cases, the cone angular aperture is large so the planar facet behaves as if it would
rotate around an axis that shares the same azimuth as the sun. Consequently the
surface gradient is somewhat undetermined in the direction orthogonal to the original
direction of the sun. In the presence of noise in the original DEM produced by ASP,
one component of the slope will thus remain uncertain and likely biased creating an
"excitation" that can numerically propagate across the scene. The wavy artifact thus
introduced in the refined DEM has no influence on the modeling of the CTX left
image by the BOA reflectance image because of the photometric invariance. How to
deal with these artifacts either by filtering or by breaking the photometric
indetermination with additional images will be the subject of future work.
<figure><img src="image_13.png" /><figcaption>Fig. 11. The invariance of the intensity for a Lambertian facet rotated such that its normal sweeps out a cone centered at the sun direction because of the invariance of the incidence angle. $\theta_0$ and $\theta$ denote the local incidence and emergence angles, respectively. $\varphi_0$ and $\varphi$ denote the azimuth angles of the sun and sensor, respectively.</figcaption></figure>

4. Conclusions and future work

This paper proposes a new method that merges photogrammetric and photoclinometric information for producing high resolution DEMs over ≈1000 km2 region of interest of Mars from orbital imagery. The DEM produced by ASP is used as an initial input of an iterative optimization procedure minimizing a photoclinometric energy function. The latter contains three terms expressing respectively the adequacy of the DEM as respect to the image intensity variation through a surface-atmosphere radiative transfer model, the integrability of the resulting surface and common characteristics for the original and refined DEMs at large scales. From this point of view the photogrammetric information is also used as a constraint in the optimization function. In the radiative transfer scheme, the RTLS surface reflectance model is
adopted and its kernel weights are supplied by Mars-Reco algorithm. The experiments
were performed using pairs of CTX images with increasing complexity. The results of
the different datasets imply that, the proposed method can produce DEMs with fewer
artifacts and better details than the photogrammetric method alone.

The results of the proposed method depend on the accuracy of the initial DEM produced by ASP and the accuracy of the RTLS kernel weights. Therefore, the future work can be two-folds. First, the patterns that appear in the shaded image derived from the refined version of the DEM when the sun direction is very different from that of the original image are to be further studied. The exact reason of the artifacts, the dependence on the incidence angle, the influence of the surface reflectance model, and how to reduce these artifacts remain an interesting problem. Second, the proposed method aims at producing DEMs for areas where the reflectance properties of the surface are supplied by independent methods as a priori information. We envision estimating conjointly and coherently the RTLS weights as well as the topography if series of images with different geometric conditions are available.

## Acknowledgments

This research was carried out under the project "I2- Mars" conjointly funded by the
Agence Nationale de la Recherche (ANR) (grant number ANR-12-IS05-0001-01) and
the National Science Foundation of China (NSFC).

## References

ASU, 2006. Mars Reconnaissance Orbiter Context Camera (CTX) Image Map, <http://global-data.mars.asu.edu/bin/ctx.pl>. (accessed 31.10.2016)

Barlow, N., 2008. Mars: An introduction to its interior, surface and atmosphere. Cambridge University Press.

Ceamanos, X., Douté, S., Fernando, J., Schmidt, F., Pinet, P., Lyapustin, A., 2013. Surface reflectance of Mars observed by CRISM/MRO: 1. Multi-angle Approach for Retrieval of Surface Reflectance from
CRISM observations (MARS-ReCO). Journal of Geophysical Research: Planets 118 (3), 514–533.

Ceamanos, X., Douté, S., Luo, B., Schmidt, F., Jouannic, G., Chanussot, J., 2011. Intercomparison and validation of techniques for spectral unmixing of hyperspectral images: A planetary case study. IEEE Transactions on Geoscience and Remote Sensing 49 (11), 4341–4358.

Comaniciu, D., Meer, P., 2002. Mean shift: A robust approach toward feature space analysis. IEEE Transactions on pattern analysis and machine intelligence 24 (5), 603–619.

Cryer, J.E., Tsai, P.-S., Shah, M., 1993. Integration of shape from x modules: Combining stereo and shading, Computer Vision and Pattern Recognition, 1993. Proceedings CVPR'93., 1993 IEEE Computer Society Conference on. IEEE, pp. 720–721.

Dorrer, E., Mayer, H., Ostrovskiy, A., Reznik, S., Neukum, G., 2004. De-and re-shading of Mars Express HRSC image data for homogenization of map relief shading. Int. Arch. Photogramm. Remote Sensing Spatial Inf. Sci 35, 1304–1308.

Dorrer, E., Zhou, X., 1998. Towards optimal relief representation from Mars imagery by combination of DEM and shape from shading. International Archives of Photogrammetry and Remote Sensing 32 (4), 156–161.

Fassold, H., Danzl, R., Schindler, K., Bischof, H., 2004. Reconstruction of archaeological finds using shape from stereo and shape from shading, Proc. 9th Computer Vision Winter Workshop, Piran, Slovenia, pp. 21–30.

Fletcher, R., 2013. Practical methods of optimization. John Wiley & Sons.

Fua, P., Leclerc, Y.G., 1995. Object-centered surface reconstruction: Combining multi-image stereo and shading. International Journal of Computer Vision 16 (1), 35–56.

Gehrke, S., 2006. Geometric and radiometric modeling of the martian surface based on object space matching and photoclinometry. ISPRS Archives, XXXVI (B4).

Grumpe, A., Belkhir, F., Wöhler, C., 2014. Construction of lunar DEMs based on reflectance modelling. Advances in Space Research 53 (12), 1735–1767.

Grumpe, A., Wöhler, C., 2014. Recovery of elevation from estimated gradient fields constrained by digital elevation maps of lower lateral resolution. ISPRS Journal of Photogrammetry and Remote Sensing 94, 37–54.

Horn, B.K., 1990. Height and gradient from shading. International journal of computer vision 5 (1), 37–75.

Jiang, C., Zhang, H., Zhang, L., Shen, H., Yuan, Q., 2016. Hyperspectral Image Denoising with a Combined Spatial and Spectral Weighted Hyperspectral Total Variation Model. Canadian Journal of Remote Sensing 42 (1), 53–72.

Kirk, R., Howington-Kraus, E., Rosiek, M., Anderson, J., Archinal, B., Becker, K., Cook, D., Galuszka, D., Geissler, P., Hare, T., 2008. Ultrahigh resolution topographic mapping of Mars with MRO HiRISE stereo images: Meter-scale slopes of candidate Phoenix landing sites. Journal of Geophysical Research: Planets 113 (E3).

Kumar, S., Kumary, M., Sukavanamy, N., Balasubramanian, R., Bhargava, R., 2009. Depth Recovery of Complex Surfaces from Texture-less Pairs of Stereo Images. ELCVIA: electronic letters on computer vision and image analysis 8 (1), 44–56.

Leclerc, Y.G., Bobick, A.F., 1991. The direct computation of height from shading, Computer Vision and Pattern Recognition, 1991. Proceedings CVPR'91., IEEE Computer Society Conference on. IEEE, pp. 552–558.

Lohse, V., Heipke, C., Kirk, R.L., 2006. Derivation of planetary topography using multi-image
shape-from-shading. *Planetary and Space Science* 54 (7), 661–674.

Luo, B., Chanussot, J., Douté, S., Zhang, L., 2013. Empirical automatic estimation of the number of endmembers in hyperspectral images. *IEEE Geoscience and Remote Sensing Letters* 10 (1), 24–28.

McGlone, C., Mikhail, E., Bethel, J., 2013. Manual of photogrammetry, Sixth ed. American Society of Photogrammetry and Remote Sensing, Bethesda, USA.

NASA, 2009. NASA Ames Stereo Pipeline, <https://ti.arc.nasa.gov/tech/asr/intelligent-robotics/ngt/stereo/>. (accessed 31.10.2016)

PDS Geosciences Node, 2009. NASA Planetary Data System MOLA PEDR Query Tool, <http://oderest.rsl.wustl.edu/GDSWeb/GDSMOLAPEDR.html>. (accessed 26.10.2016)

O'Hara, R., Barnes, D., 2012. A new shape from shading technique with application to Mars Express HRSC images. *ISPRS journal of photogrammetry and remote sensing* 67, 27–34.

Öhman, T., 2013. A beginer's guide to stereo-derived DEM production and analysis using ISIS, ASP, and ArcMap. Lunar and Planetary Institute, internal report, 30 pp.,doi: 10.13140/RG.2.1.1743.7609., <http://www.lpi.usra.edu/lunar/tools/dems/Ohman_2013_ISIS-ASP-ArcMap workflow.pdf>. (accessed 31.10.2016)

Pien, H.H., Gauch, J.M., 1995. A variational approach to multi-sensor fusion of images. Applied Intelligence 5 (3), 217–235.

Samaras, D., Metaxas, D., Fua, P., Leclerc, Y.G., 2000. Variable albedo surface reconstruction from stereo and shape from shading, Computer Vision and Pattern Recognition, 2000. Proceedings. IEEE Conference on. IEEE, pp. 480–487.

Shean, D.E., Alexandrov, O., Moratto, Z.M., Smith, B.E., Joughin, I.R., Porter, C., Morin, P., 2016. An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery. *ISPRS Journal of Photogrammetry and Remote Sensing* 116, 101–117.

Smith, D.E., Zuber, M.T., Frey, H.V., Garvin, J.B., Head, J.W., Muhleman, D.O., Pettengill, G.H., Phillips, R.J., Solomon, S.C., Zwally, H.J., 2001. Mars Orbiter Laser Altimeter: Experiment summary after the first year of global mapping of Mars. *Journal of Geophysical Research: Planets* 106 (E10), 23689–23722.

Spiga, A., Forget, F., 2008. Fast and accurate estimation of solar irradiance on Martian slopes. *Geophysical Research Letters* 35 (15).

Stamnes, K., Tsay, S.-C., Wiscombe, W., Jayaweera, K., 1988. Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Applied optics 27 (12), 2502–2509.

USGS, 1992. Integrated Software for Imagers and Spectrometers, <https://isis.astrogeology.usgs.gov/index.html>. (accessed 02.11.2016)

Vincendon, M., 2013. Mars surface phase function constrained by orbital observations. *Planetary and Space Science* 76, 87–95.

Xue, W., Zhang, L., Mou, X., Bovik, A.C., 2014. Gradient magnitude similarity deviation: a highly efficient perceptual image quality index. *IEEE Transactions on Image Processing* 23 (2), 684–695.

Zhang, R., Tsai, P.-S., Cryer, J.E., Shah, M., 1999. Shape-from-shading: a survey. *IEEE transactions on pattern analysis and machine intelligence* 21 (8), 690–706.