Nathalie Bertrand

Concurrent Games with Arbitrarily Many Players

Keywords: 2012 ACM Subject Classification Theory of computation → Verification concurrent games, parameterized verification Digital Object Identifier 10, 4230/LIPIcs, MFCS, 2020, 1 Category Invited Talk

Traditional concurrent games on graphs involve a fixed number of players, who take decisions simultaneously, determining the next state of the game. With Anirban Majumdar and Patricia Bouyer, we introduced a parameterized variant of concurrent games on graphs, where the parameter is precisely the number of players. Parameterized concurrent games are described by finite graphs, in which the transitions bear finite-word languages to describe the possible move combinations that lead from one vertex to another.

We report on results on two problems for such concurrent games with arbitrary many players. To start with, we studied the problem of determining whether the first player, say Eve, has a strategy to ensure a reachability objective against any strategy profile of her opponents as a coalition. In particular Eve's strategy should be independent of the number of opponents she actually has. We establish the precise complexities of the problem for reachability objectives. Second, we considered a synthesis problem, where one aims at designing a strategy for each of the (arbitrarily many) players so as to achieve a common objective. For safety objectives, we show that this kind of distributed synthesis problem is decidable.

Motivation

We introduce and study concurrent games in which the number of players is a priori unknown. Games with arbitrarily many players seem particularly relevant to model modern distributed systems. A first typical situation is the one of a global server, answering requests from an arbitrary number of clients. One can also think of a fleet of drones trying to cooperate to achieve a common goal. Wireless sensors networks and ant colonies are more application examples of games with arbitrarily many players.

Games with arbitrarily many players

Starting from concurrent games with a fixed number of players [START_REF] Luca De Alfaro | Concurrent reachability games[END_REF][START_REF] Alur | Alternating-time temporal logic[END_REF], a natural idea is to define concurrent with arbitrarily many players by equipping edges of the arena with languages of finite words. For instance, an edge from vertex v to vertex v can be labelled with a language L, representing the situation where, if there are k players, and in v, player i chooses action a i , and a 1 • • • a k ∈ L, then the next vertex will be v . As an example, L can be the regular language described by the regular expression a(Σ 2 ) * + (bb) * ; with six players that all choose b, or with seven players if the first one chooses a while the choices of all others are arbitrary, will lead to v .

1:2 Concurrent Games with Arbitrarily Many Players

Note that the number of players k is unknown to them, but fixed all along the play. Choosing k, and resolving the nondeterminism is performed by the (adversarial) environment. Since the number of players is a parameter, we refer to the arenas as parameterized arenas, that we now formally define.

Definition 1 (Parameterized arena). A parameterized arena is a tuple

A = (V, Σ, ∆) with V is a finite set of vertices; Σ is a finite alphabet of actions; ∆ : V × V → 2 Σ + is a partial transition function.
Example 2. The notion of parameterized arena is illustrated on an example, depicted in Figure 1. Here, the alphabet of actions is Σ = {a, b}, and for instance, the language that labels the edge from v 3 to v 4 is ba + abab + . Thus, when in vertex v 3 , if either there are two players and Player 1 plays b while Player 2 plays a, or if there are at least 4 players and their actions form a word in abab + , then the game moves to v 4 . Note that this arena is nondeterministic: for instance, from vertex v 2 , any word in a(aa) * can lead to v 0 , v 3 and v 5 . Also, in this example, all languages are regular, and are thus denoted by regular expressions. In the whole paper, we assume that arenas are complete: from any vertex v, for any non-empty word w ∈ Σ + , there exists an edge v L -→ v with w ∈ L. For conciseness, the examples -as above-might depict incomplete arenas; a sink state can be added so as to obtain complete arenas.

Definition 3 (Strategies and induced plays). A strategy for Player i in the arena

A = (V, Σ, ∆) is a function σ i : V + → Σ.
An infinite strategy profile π = (σ 1 , σ 2 , • • • ) induces the plays:

Plays A (π) = k Plays A (σ 1 , σ 2 , • • • , σ k ) = k {v 0 v 1 v 2 • • • | ∀j ≥ 0, σ 1 (v 0 • • • v j ) • • • σ k (v 0 • • • v j ) ∈ ∆(v j , v j+1 )}
In words, a strategy dictates which action to play depending on the sequence of vertices seen so far. The plays induced by a strategy profile are formed of induced k-plays for each possible number of players k. An induced k-play satisfies that at each step j, the word, obtained by concatenating the actions prescribed by the strategies for players from 1 to k, belongs to the language labelling the edge from v j to v j+1 . The initial choice of k, and the resolution of nondeterminism during the play are taken care of by an adversarial environment.

1: Example 4. Let us illustrate the notions of strategies and plays on Figure 2. Notice that in this arena, the actions of all players but Player 1 are irrelevant: indeed the languages are particularly simple aΣ, aΣ + , etc.

3 v0 v1 v2 v3 v4 v5 v6 a Σ a Σ ≥ 2 a Σ + a Σ + a Σ + b Σ ≥ 2 b Σ + a Σ ≥ 2 Σ ≥2 Σ ≥2
Assuming the game starts in v 0 , an example of strategy for Player 1 is the following:

σ 1 (v 0 ) = σ 1 (v 0 v 1 ) = σ 1 (v 0 v 2 ) = σ 1 (v 0 v 1 v 3 ) = a σ 1 (v 0 v 2 v 3 ) = b Examples of plays consistent with σ 1 are v 0 v 1 v 3 v 4 if k = 2 and v 0 v 2 v 3 v 4 if k = 3.
Indeed, if we annotate plays with actions of each players, these plays are for instance obtained by

v 0 aa -→ v 1 ab -→ v 3 aa -→ v 4 and v 0 aab --→ v 2 abb --→ v 3 baa --→ v 4 , respectively.
Even if the actions of players 2 to k are irrelevant, the arena is nondeterministic: in the very first step, under σ 1 , the environment resolves the nondeterminism between going to v 6 , or progressing to v 2 or v 3 (depending on the number of opponents). For a fixed number of players k, there are thus two plays induced by σ 1 .

Definition 5 (Game). A game G = (A, Win) is an arena equipped with a set of infinite plays: Win

⊆ V ω .

Typical examples of winning conditions that we use in this paper are

Reachability: For a target set

T ⊆ V , Win = {v 0 v 1 • • • |∃i : v i ∈ T }; Safety: For a safe set S ⊆ V , Win = {v 0 v 1 • • • |∀i : v i ∈ S}.

Outline

Parameterized arenas raise many interesting problems, and in this article we focus on two of them. In the first problem, Player 1 is distinguished, and she aims at achieving an objective independently of the number of opponents she has and of the strategies they play. In the second problem, all players try to achieve an objective as a coalition, not knowing however a priori how many they are. In the next two sections, we formalize these problems, and give decidability and complexity results.

One player against all

We first consider a setting in which player 1 aims at achieving an objective independently on the number of her adversaries, and whatever strategy they play. This situation can be motivated for instance by a scenario in which a server aims at answering requests by arbitrarily many clients. Formally

Eve vs rest of the world Input: A parameterized arena A and a winning condition Win. Question: The PSPACE-hardness (see [START_REF] Bertrand | Concurrent parameterized games[END_REF]) is obtained via a natural reduction from QBF [START_REF] Stockmeyer | Word problems requiring exponential time (preliminary report)[END_REF], and we focus now on sketching the proof of PSPACE membership.

∃σ 1 ∀k∀σ 2 • • • σ k Plays A (σ 1 , σ 2 , • • • , σ k ) ⊆ Win? M F C S
Recall that Eve must win against her opponents playing as a coalition and also against the environment that chooses the number of players before the play starts, and resolves potential nondeterminism. It seems thus quite natural to reduce to a 2-player game (see for instance [START_REF] Grädel | Automata, Logics, and Infinite Games: A Guide to Current Research[END_REF]Chapter 12] for a gentle introduction to such games). As we have seen in Example 4, in order to win, Eve must gain information on the number of opponents she actually has. We thus build the knowledge game, a 2-player turn-based game tracking Eve's knowledge. The knowledge game starts in a state (v 0 , N) reflecting that the parameterized game is in v 0 , and Eve has no information on the number of her opponents. The construction of the knowledge game is depicted in Figure 3: on the left-hand side is a subgame of the parameterized arena, and on the right-hand side part of the corresponding knowledge game. If K ⊆ N is the current knowledge, from a vertex (v, K) of the existential player (represented with rounded rectangles), she chooses an action in Σ, leading to a vertex of the universal player (represented with diamonds). The universal player resolves nondeterminism (if any), leading to a vertex (v i , K i ), where K i represents the updated knowledge. Here L denotes the set of lengths of words in L: L = {|w| | w ∈ L}. Assuming that she played a, and the play moves to v i , Eve updates the actual number of opponents to K ∩ a -1 L i . While building the knowledge game, only existential vertices (v i , K i ) with non-empty K i are constructed. Also, only universal vertices corresponding to a feasible action of Eve are built.

v v1 v2 v,K v1,K∩ a -1 L1 v2,K∩ a -1 L2 v1,K∩ b -1 L1 v2,K∩ b -1 L2 L1 L2 a b
Figure 4 provides the knowledge game for the example from Figure 2.

v0,N v1,{1} v2,[2,∞[ v6,N v3,{1} v3,[2,∞[ v4,{1} v5,{1} v4,[2,∞[ v5,[2,

N. Bertrand 1:5

If we transfer the winning condition Win to the knowledge game, one can show that Eve has a winning strategy in the parameterized arena if and only if the existential player wins the knowledge game.

Example 7. Back to the parameterized arena of Figure 2, consider the objective to reach the set {v 4 , v 6 }. On the corresponding knowledge game (see Figure 4) the objective of the existential player is thus to reach any vertex (v, K) with v ∈ {v 4 , v 6 }. An obvious winning strategy to do so is to play b in the first step. However, even if a is the first action, the existential player can ensure reaching the target set: the chosen action should be a after going through v 1 , and b after going through v 2 . Thus, the strategy σ 1 described in Example 4 ensures Eve to reach {v 4 , v 6 } independently of the number of her opponents.

To solve Eve vs rest of the world, it thus suffices to solve the 2-player turn-based knowledge game. Starting with a parameterized arena A with regular languages on the edges, the corresponding knowledge game is at most exponential in |A|. This is due to the fact that vertices in the knowledge game encode the knowledge Eve has, which is obtained by intersecting lengths of words in the initially given languages. These subsets of N are all semilinear [START_REF] Parikh | On context-free languages[END_REF], and there can be at most exponentially many, corresponding to the possible combinations of intersections. Constructing the knowledge game, and solving it for a reachability objective would thus yield an exponential time algorithm.

To obtain the announced PSPACE complexity upper bound, we show that storing the whole knowledge game is not necessary. In constrast, taking a dynamic programming approach, it is sufficient to store only subgames that are polynomial in the size of the parameterized arena. Each subgame is rooted at some existential vertex (v, K) and stops as soon as, either the target set T is reached (with arbitrary knowledge), or the knowledge changes to some K K. In such a subgame, there are at most polynomially many vertices, and the objective of the existential player is to reach vertices (v , K ) that are winning. Such winning vertices are computing recursively, and the recursion depth is polynomially bounded. The interested reader can find more details on this polynomial space procedure in [START_REF] Bertrand | Concurrent parameterized games[END_REF].

Strategy synthesis for a coalition

After the one player against all setting, we now consider the case where agents want to collectively achieve a goal, independently of the number they actually are. Formally

Synthesis for arbitrarily-large coalition

Input: A parameterized arena A and a winning condition Win.

Question: ∃σ 1 , σ 2 • • • ∀k Plays A (σ 1 , σ 2 , • • • , σ k ) ⊆ Win?
Theorem 8. When Win is a safety objective, Synthesis for arbitrarily-large coalition is in EXPSPACE and PSPACE-hard.

The PSPACE-hardness proof for Eve vs rest of the world can be adapted to obtain the same lower-bound here. In the sequel, we explain how we establish the EXPSPACE complexity upper bound. Closing the complexity gap is currently on our agenda. Example 9. Consider the nondeterministic parameterized arena from Figure 5, and assume the winning objective Win is to avoid the sink vertices v 2 and v 3 . Assuming the game starts at v 0 , to ensure this safety condition as a coalition, the players can apply the following memoryless strategy profile: in v 0 , all players but Player 1 play a, and in v 1 , all players play a. Under this strategy profile, for any number of players k ≥ 1, all consistent plays avoid v 2 and v 3 . In the sequel, we focus on safety objectives specified by a safe set of vertices S ⊆ V , and sketch how to decide whether there exists a strategy profile π = (σ 1 , σ 2 , • • • ) that ensures to stay within S independently of the number of players.

Observe that a strategy profile for arbitrarily many players can equivalently be seen as a map π : V + → Σ ω from sequences of vertices to infinite words. Such a profile operates as follows: after observing a history

v 0 v 1 v 2 • • • v j , if there are k players, the effect of profile π is the prefix of length k of π(v 0 v 1 v 2 • • • v j ).
As a first step towards deciding the Synthesis for arbitrarily-large coalition problem for safety objectives, we first unfold the arena into a tree T , and stop a branch as soon as either it reaches a node labelled with an unsafe vertex v / ∈ S, or it reaches a node labelled with a vertex v that has an ancestor with same label. This construction is illustrated in Figure 6 on the example of Figure 5. The two left-most branches are stopped because of the repetition of the label v 0 , and the two right-most branches are stopped because of the labelling by an unsafe vertex (v 2 or v 3 ). The size (i.e., number of nodes) of the tree unfolding T can be exponential in the number of vertices of the parameterized arena A, however, its branching degree and height are linear in the size of A.

The unfolding T can itself be seen as a parameterized arena, in which the objective of the coalition is to avoid unsafe branches that end in a node labelled with an unsafe vertex. Strategies of the coalition in T map inner nodes (nodes that are not leaves) to ω-words. One can show that the coalition has a winning strategy profile in A for the safety objective defined by S, if and only if it has a winning strategy profile in T to avoid unsafe branches. Intuitively, from a winning strategy profile in the tree, one can build a winning strategy profile in the arena with finite-memory bounded by the height of the tree. On our example, once we show there is a winning strategy profile π T in the tree, one can define a winning strategy profile π A in the arena by: π A (V + v 0 ) = π T (v 0 ) and π A (V + v 1 ) = π T (v 0 v 1 ). Clearly enough, this justifies stopping a branch as soon as there is a repetition.

The second and most involved step is to characterize, at the tree unfolding level, the winning strategy profiles. If m is the number of inner nodes of T , we show that one can effectively build a deterministic safety automaton B over Σ m (thus reading one letter of the prescribed strategy at each inner node simultaneously) that accepts all infinite words in (Σ m ) ω that correspond to winning strategies in T . Altogether, the Synthesis for arbitrarily-large coalition problem reduces to checking non-emptiness of B. The latter being at most doubly exponential in the size of A, we obtain an EXPSPACE complexity upper-bound.

Discussion

In this article, we reported on recent results about a new model of concurrent games, where the number of players is arbitrary [START_REF] Bertrand | Concurrent parameterized games[END_REF]. Concurrent games with arbitrarily many players extend 2-player concurrent games, and more generally concurrent games with a fixed number of players. The edges in the arenas are equipped with languages over finite words. Such parameterized arenas can represent at once a denumerable number of standard arenas, each with a fixed number of players. They enable the definition and study of a number of parameterized game-theoretic problems.

We first considered a setting in which one player, say Eve, wants to achieve an objective independently of how many opponents she has, and whatever strategy they choose. We show that, for reachability objectives, deciding the existence of a uniform winning strategy for Eve against the rest of the world, is a PSPACE-complete problem. Second, we started to explore a synthesis problem, in which all players want to achieve an objective as a coalition. The difficulty here lies in the fact that they do not know a priori how many they are. For safety objectives, the existence of a coalition strategy is in EXPSPACE and PSPACE-hard.

We believe our preliminary work on parameterized arenas opens up many research paths. On the theoretical side, we currently put our effort on the coalition synthesis problem for reachability objectives. As an example, consider the parameterized arena from Figure 7, and assume the objective is to reach vertex v 1 . Without knowing a priori how many they are, the players can collectively achieve this objective with the following profile: as long as the play is in v 0 , at step i, Player i plays b while all other players play a. Under this strategy for the coalition, if there are k players, at step k the play moves from v 0 to v 1 . In constrast to safety objectives, synthesizing such a symbolic strategy profile (or even deciding its existence) calls for more involved techniques. Studying other solution concepts [START_REF] Grädel | Solution concepts and algorithms for infinite multiplayer games[END_REF] such as Nash equilibria [START_REF] Ummels | The complexity of Nash equilibria in limit-average games[END_REF][START_REF] Bouyer | Pure Nash equilibria in concurrent games[END_REF], secure equilibria [START_REF] Chatterjee | Games with secure equilibria[END_REF] and subgame perfect equilibria [START_REF] Selten | Spieltheoretische Behandlung eines Oligopolmodells mit Nachfrageträgheit[END_REF][START_REF] Brihaye | On subgame perfection in quantitative reachability games[END_REF] is also on our agenda. On the practical side, we believe parameterized arenas could be used to represent a variety of distributed systems. We would be quite interested in exploring potential applications of this model that we find fascinating.

)Figure 1

 1 Figure 1 An example of a parameterized arena.

Figure 2 A

 2 Figure 2 A simple concurrent parameterized arena.

Figure 3

 3 Figure 3 Construction of the knowledge game.

Figure 4

 4 Figure 4 Knowledge game for the arena from Figure 2.

ΣFigure 5

 5 Figure 5 Parameterized arena for the synthesis of a coalition strategy.

aFigure 6

 6 Figure 6 Finite unfolding of the arena from Figure 5 with S = {v0, v1}.

Figure 7

 7 Figure 7 Synthesis for a reachability objective.

2 0 2 0 1:4 Concurrent Games with Arbitrarily Many Players Theorem

6. When Win is a reachability objective, Eve vs rest of the world is a PSPACEcomplete problem.

Acknowledgements This paper is based on partly published results obtained in a collaboration with Patricia Bouyer and Anirban Majumdar.


Keywords: Mars, surface processes, topography, image processing

I. INTRODUCTION

s one of the four terrestrial planets, Mars has drawn people's attention since ancient times [START_REF] Barlow | Mars: An introduction to its interior, surface and atmosphere[END_REF]. There have been over 40 space missions to Mars since 1960s, and the exploration enthusiasm continues with ESA's Trace Gas Orbiter and NASA's Insight missions. In the last decade, the large datasets acquired by NASA's Mars Reconnaissance Orbiter (MRO) and ESA's Mars Express have allowed us to better understand the atmosphere, surface, and sub-surface of Mars. Producing 3D surface models, the so-called digital elevation models (DEMs), from remote sensing measurements is fundamental for geology investigations, meteorology simulations, radar sounding interpretation, and image radiometric-geometric correction before analysis. Ongoing research aims at reaching a high horizontal and vertical precision (meter or sub-meter scale) for the DEM because this is critical for certain tasks, e.g. choosing a landing site for exploration rovers [START_REF] Beyer | Meter-Scale Slopes of Candidate InSight Landing Sites from Point Photoclinometry[END_REF].

So far, for planetary DEM construction, two technologies are available. One is based on range measurements, using instruments such as Laser Altimeter, e.g. Lunar Orbiter Laser Altimeter (LOLA) and Mars Orbiter Laser Altimeter (MOLA). Usually they can provide global DEM at very high vertical resolution but at very low horizontal resolution due to the interpolation process of a limited number of tracks. For the DEM produced by MOLA, the vertical resolution is less than 1m, while the spatial resolution is approximately 168m.pixel -1 [START_REF] Smith | Mars Orbiter Laser Altimeter: Experiment summary after the first year of global mapping of Mars[END_REF].

The second technique adopts images' geometrical and intensity information to derive the DEM, namely photogrammetry (or stereo) and photoclinometry (or shape from shading, SFS). As an effective and cost-effective way, photogrammetry is widely adopted for producing DEMs on non-desert areas on Earth [START_REF] Eckert | Comparison of automatic DSM generation modules by processing IKONOS stereo data of an urban area[END_REF][START_REF] Mcglone | Manual of photogrammetry[END_REF]. The technique searches a given scene point in at least two images (matching process) and then inversely solves the collinearity equations which relates the point in the scene and the conjugated pixels in the images taking into account a camera model with its parameters (triangulation process). Usually ground control points, which coordinates are precisely known by GPS measurements, are integrated in the process so as to better constrain the image relative orientations and the internal camera parameters ("bundle adjustment"). Ideally, photogrammetry methods can yield DEMs at lateral resolutions two times lower than the image resolution. However, for missions on Mars, the resolution of the orbital images and the camera parameters' precision are not as high as those on Earth. In addition the density of ground control points is low and there are less areas containing features or recognizable intensity patterns for the matching process. This makes photogrammetry method alone not able to produce high resolution DEMs devoid of various artifacts and noise. For fine details, photoclinometry methods [START_REF] Zhang | Shape-from-shading: a survey[END_REF][START_REF] Horn | Height and gradient from shading[END_REF] seem to be a good solution which is reported to be able to refine the DEM with comparable resolution as the image [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF]. The first principle of photoclinometry is the sensitivity of the image intensity variations to the altitude gradients of the surface especially in the sun direction. Therefore inverting the image intensity formation process allows us to derive the gradients as well as the relative height. In this way, photoclinometry does not need multiple images, which avoids the problems for image matching.

Photoclinometry alone also has its drawbacks, as it cannot produce absolute height fields and the model is complex involving the bidirectional reflectance properties of the surface, a radiative transfer model, and regularizations. To alleviate this, photoclinometry can be combined with Laser Altimeter data [START_REF] Soderblom | Accurate fine-scale topography for the Martian south polar region from combining MOLA profiles and MOC NA images[END_REF][START_REF] Soderblom | Meter-Scale 3-D models of the martian surface from combining MOC and MOLA data[END_REF] or photogrammetry which represents the trend in the literature. Generally, previous studies combine photoclinometry with photogrammetry at three levels, i.e. data-level, model-level and regularization-level methods.

For data-level methods, photogrammetry and photoclinometry are individually performed and then the Small-scale Topographical Characterization of the Martian Surface with In-orbit Imagery.

A results are combined in the Fourier domain [START_REF] Cryer | Integration of shape from x modules: combining stereo and shading[END_REF], the original object space [START_REF] Samaras | Variable albedo surface reconstruction from stereo and shape from shading[END_REF], and even via a neural network [START_REF] Kumar | Depth Recovery of Complex Surfaces from Texture-less Pairs of Stereo Images[END_REF]. The model-level methods [START_REF] Lohse | Derivation of planetary topography using multi-image shape-from-shading[END_REF][START_REF] Gehrke | Geometric and radiometric modeling of the martian surface based on object space matching and photoclinometry[END_REF][START_REF] Fua | Object-centered surface reconstruction: Combining multi-image stereo and shading[END_REF] combine photogrammetry and photoclinometry into a single model, in which the position and intensity of a certain pixel in the image are determined, among other parameters, by the height of the point in the object space.

For regularization-level methods [START_REF] Pien | A variational approach to multi-sensor fusion of images[END_REF][START_REF] Dorrer | Towards optimal relief representation from Mars imagery by combination of DEM and shape from shading[END_REF][START_REF] Dorrer | De-and re-shading of Mars Express HRSC image data for homogenization of map relief shading[END_REF][START_REF] Fassold | Reconstruction of archaeological finds using shape from stereo and shape from shading[END_REF][START_REF] Grumpe | Recovery of elevation from estimated gradient fields constrained by digital elevation maps of lower lateral resolution[END_REF], the low-resolution DEM produced by photogrammetry is integrated in a regularization term in the photoclinometry model. In addition the low-resolution DEM constitutes an a priori solution of the latter. Should the difference of resolution between the low-resolution DEM and the image exceeds a certain threshold, the photoclinometic workflow can be integrated in a hierarchical iterative strategy during which the effective spatial and vertical resolutions of the reconstituted topography increase progressively [START_REF] Wu | Construction of pixel-level resolution DEMs from monocular images by shape and albedo from shading constrained with low-resolution DEM[END_REF]. Existing hybrid methods often produce better results than those by using photogrammery or photoclinometry alone. However, a single model that can generate a DEM at the resolution of the image without strong requirements while not being restricted to some type of terrains and being free of artifacts is still lacking. These requirements include: use of multi-images, oblique-views of the camera [START_REF] Lohse | Derivation of planetary topography using multi-image shape-from-shading[END_REF][START_REF] Gehrke | Geometric and radiometric modeling of the martian surface based on object space matching and photoclinometry[END_REF][START_REF] Fua | Object-centered surface reconstruction: Combining multi-image stereo and shading[END_REF], and a Lambertial assumption for the BRDF of the surface [START_REF] Dorrer | De-and re-shading of Mars Express HRSC image data for homogenization of map relief shading[END_REF][START_REF] Barron | High-frequency shape and albedo from shading using natural image statistics[END_REF]. As a matter of fact sensors around Mars, e.g., the Context Camera (CTX) or the High Resolution Imaging Science Experiment (HiRISE) both aboard the Mars Reconnaissance Orbiter (MRO) do not allow presently the near simultaneous acquisition of more than two images. When multiple images are compiled, one has to face phenomena of diachronism (different shooting times, variation of brightness, movement of shadows) and strong affine transformations (very different angles of view). In a previous work [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF], we proposed a method combining photogrammetric coarse DEMs and photoclinometry at the regularization-level to produce large high resolution DEMs from individual Mars images. The method integrates an intensity model of the image based on a novel radiative transfer scheme and a realistic bidirectional reflectance distribution function (BRDF) model of the surface, namely Ross-Thick Li-Sparse (RTLS) model [START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF]. These characteristics make the method being able to recover details missing or degraded in the photogrammetrically generated DEMs such as the ones produced by Ames Stereo Pipeline (ASP) [START_REF] Shean | An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery[END_REF]. However, at the time of first release, the method still suffered from wave-like artifacts approximately orthogonal to the direction of the sun because of intrinsic limitations of the photoclinometry method. Some experiments also showed that the quality of the initial DEM has some effect on the final refined DEM precluding the use of the method for certain areas with very poor initial DEMs. Besides, the validation of the method was restricted to a minimum and we had no investigation about the characteristic spatial scales controlled by the photoclinometric process. The previous problems are largely alleviated by the novel approach proposed in this work. The latter is much more robust to the initial DEM conditions and is nearly free of artifacts with better vertical precisions. The contributions of the present paper can be summarized as follows: 1. We solve the artifact problem by using two strategies simultaneously. First we propose a cartesian coordinate system aligned with the directions parallel and perpendicular to the sun azimuth instead of using the initial frame of the image [START_REF] Kirk | A fast finite-element algorithm for two-dimensional photoclinometry[END_REF]. In addition, we split the integration regularization term into two terms with a weighting factor to pose more penalization on the direction normal to the sun azimuth. Second, a Gaussian convolution scheme is integrated into the photogrammetric constraint in order to take into account any difference of spatial resolution between the initial and refined DEMs. Indeed the former is often filtered in order to eliminate the noise at high-frequencies at the expense of the resolution before it can be used as the initial solution of the photoclinometric processing. 2. The proposed method is validated and carefully characterized on images acquired above Gusev Crater by both CTX and HiRISE. The Gusev Crater, 103 miles across and lying just south of the Martian equator, was the landing site for the Mars Exploration Rover Spirit in January 2004. The validation and the characterization of our products are based on inter-comparison and multi-scale analysis. 3. A multi-resolution strategy is proposed when applying the method to poor initial DEMs containing large-area with strong artifacts, making the method even more robust to the initial DEM conditions. Different from other iterative works [START_REF] Gehrke | Geometric and radiometric modeling of the martian surface based on object space matching and photoclinometry[END_REF][START_REF] Wu | Construction of pixel-level resolution DEMs from monocular images by shape and albedo from shading constrained with low-resolution DEM[END_REF][START_REF] Capanna | Three-dimensional reconstruction using multiresolution photoclinometry by deformation[END_REF], a simplification and optimization effort is made to reach maximum numerical efficiency while maintaining a high precision based on analyzing spatial scales controlled by the photoclinometric process.

For that purpose we consider images covering the Russell Dune, an interesting and challenging region. 4. Finally to illustrate another added value of our method, we present topographical products derived from the refined DEMs aimed at understanding planetary processes. The article is organized as follows. The proposed method is described in Section II. In Section III, various experiments and analysis are presented and discussed. Eventually, our results are summarized and the conclusion is drawn in Section IV.

II. METHODS

The steps followed by the proposed method are illustrated in Fig. 1. The inputs include the surface photometric parameters that can be derived from multi-angular in-orbit observations [START_REF] Ceamanos | Surface reflectance of Mars observed by CRISM/MRO: 1. Multi-angle Approach for Retrieval of Surface Reflectance from CRISM observations (MARS-ReCO)[END_REF] and the initial DEM usually produced by photogrammetry tools. The core part of the method is the optimization model and the numerical algorithms in the sun coordinates. The model is composed of three terms, i.e., the image irradiance equation, the improved integrability constraint and the improved photogrammetry constraint, which will be discussed in detail as follows.

A. Mathematical and Physical Basis of the Photoclinometric Model

Our goal is to recover a refined DEM z from a low-resolution one 𝐳 " , using a single image I. Suppose we have a reflectance function 𝑅(𝐩, 𝐪) depending on the surface orientations which is represented by the gradient of the surface in the horizontal and vertical axis of the image as p and q. The basic equation in the photoclinometry scheme is the image irradiance equation expressed as: 𝛽𝑅(𝐩, 𝐪) = 𝐈, where 𝛽 is a scalar factor relying on the camera amplifications and other factors. The equation is underdetermined which means we have more unknowns than the number of the equations. To solve the underdetermined problem, regularization terms are needed, such as the smoothness constraint, the integrability constraint, the gradient constraint, the unit normal constraint and so on [START_REF] Zhang | Shape-from-shading: a survey[END_REF]. For previous studies in photoclinometry, some methods [START_REF] Grumpe | Recovery of elevation from estimated gradient fields constrained by digital elevation maps of lower lateral resolution[END_REF][START_REF] Grumpe | Construction of lunar DEMs based on reflectance modelling[END_REF] get p and q first, then acquire z by integrating p and q, while other methods [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF][START_REF] Fassold | Reconstruction of archaeological finds using shape from stereo and shape from shading[END_REF][START_REF] Leclerc | The direct computation of height from shading[END_REF] include z in the optimization scheme and solve for z directly. We prefer the second choice since it avoids the accumulation of errors and is more concise in the formulation. Therefore, the proposed energy function is expressed as: 𝟐 (1) where the first regularization term with parameter 𝜆 A is the improved integrability constraint in the sun coordinates that will be detailed in Section II-B. The second regularization term with parameter 𝜆 C is the improved photogrammetry constraint that will be detailed in Section II-C. Here, v and w are just p and q expressed in the sun coordinates (see Section II-B). Note that zx and zy (respectively zv and zw) are defined as the derivative of z along the horizontal (x) and vertical (y) axis of the image (respectively along the directions parallel (v) and perpendicular (w) to the sun).

𝑬 = 𝟏 𝟐 (𝜷 • 𝑹(𝐯, 𝐰) -𝐈) 𝟐 + 𝝀 𝟏 𝟐 [(𝟏 -𝜹)(𝐯 -𝐳 𝒗 ) 𝟐 + 𝜹(𝐰 - 𝐳 𝒘 ) 𝟐 ] + 𝝀 𝟐 𝟐 (𝐁 • 𝐅 * 𝐳 -𝐁 • 𝐅 * 𝐳 𝟎 )
For the reflectance function R, we employ the bottom-of-atmosphere (BOA) reflectance for a balance between complexity and accuracy. The BOA reflectance factor [START_REF] Ceamanos | Surface reflectance of Mars observed by CRISM/MRO: 1. Multi-angle Approach for Retrieval of Surface Reflectance from CRISM observations (MARS-ReCO)[END_REF] can be expressed as: 𝑹 𝑩𝑶𝑨 = 𝑺𝝁 𝟎 𝒆 J𝝉 𝟎 /𝝁 𝟎 𝝆(𝒔 𝟎 , 𝒔, 𝒗, 𝒘) + 𝑺𝝁 𝟎 𝒆 J𝝉 𝟎 /𝝁 𝟎 𝜶𝒄 𝟎 𝝆 𝟏 (𝝁)𝝆 𝟏 (𝝁 𝟎 ) + 𝜶 𝝅 × ∫ 𝑫 𝒔 𝛀 V W𝒔 𝟎 , 𝒔 ′ X 𝝆 W𝒔 ′ , 𝒔, 𝒗, 𝒘X 𝝁 ′ 𝒅𝒔 ′ where the dependency on wavelength is omitted for simplicity. The equation contains three components that are the reflection of direct sunlight, the multiple scattering between surface and atmosphere and the reflection of the sky diffuse illumination (both the sun aureole and the isotropic component), respectively. The parameter S is the extraterrestrial solar spectral irradiance. The vector 𝑠 " represents the sun direction that is depicted by the zenith angle 𝜃 " and the azimuth angle 𝜑 " , so as the view direction s depicted by 𝜃 and 𝜑. The variables 𝜇 " and 𝜇 are the cosine of 𝜃 " and 𝜃 respectively. The quantity 𝜌 is defined as π𝑓 (where f is the BRDF function that will be defined later). The parameter 𝜏 " is the atmospheric opacity, c " is the spherical albedo of the atmosphere, 𝛼 is the multiple reflection factor which is defined as (1 -𝑞(𝜇 " )𝑐 " ) JA (where 𝑞(𝜇 " ) is the surface albedo). The integration functions 𝑑𝑠 ′ (where the sign (+) and (-) represent downward and upward directions, respectively), respectively. The path radiance incident on the surface is noted 𝐷 l .

For the BRDF function, we adopt the Ross-Thick Li-Sparse (RTLS) model [START_REF] Lucht | An algorithm for the retrieval of albedo from space using semiempirical BRDF models[END_REF] that is not only realistic but also linear as it is defined as: 𝝆(𝒔 𝟎 , 𝒔, 𝒗, 𝒘) = 𝒌 𝑳 + 𝒌 𝑮 𝒇 𝑮 (𝒔 𝟎 , 𝒔, 𝒗, 𝒘) + 𝒌 𝑽 𝒇 𝒗 (𝒔 𝟎 , 𝒔, 𝒗, 𝒘) where the three terms play the roles of Lambertian (L), geometric (G), and volumetric (V) components, respectively. The weights 𝑘 s , 𝑘 t , and 𝑘 x of the three components are determined by the surface reflectance properties. The functions 𝑓 t and 𝑓 y are predefined generic kernels. The RTLS model comes from the processing of multi-angular sequences of CRISM images (EPF observations) acquired at or close to the Martian scenes of interest [START_REF] Ceamanos | Surface reflectance of Mars observed by CRISM/MRO: 1. Multi-angle Approach for Retrieval of Surface Reflectance from CRISM observations (MARS-ReCO)[END_REF]. Estimations of atmospheric opacities are extracted from Martian climatic databases or from the processing of EPF observations simultaneously acquired with the CTX or HiRISE images. Note that the equation giving RBOA assumes a flat surface for the surface-atmosphere multiple reflections to make the calculation more tractable. For efficiency, all the components in (3) except 𝜏 " , 𝑠 " , s, p, q, 𝑘 s , 𝑘 t , and 𝑘 x are retrieved using a look-up table which is precomputed using the Discrete Ordinates Radiative Transfer Program for a Multi-layered Plane Parallel Medium (DISORT) [START_REF] Stamnes | Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media[END_REF].

B. Integrability Constraint

According to [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF], 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 processed image. These effects originate from the photometric indetermination. It is a relative invariance of the intensity reflected by a planar facet of terrain when it behaves like a Lambertian reflector and when its orientation is changed such that its normal sweeps out a cone centered at the sun direction (constant local incidence angle). Experiments conducted in [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF] indicate that the waves are excited by the presence of noise or other defects in the original DEM. In order to minimize these artifacts we propose two modifications in the optimization scheme: introducing a new coordinate system and splitting the integrability constraint.

Instead of using the initial coordinate system of the image (x, y) we shall prefer a rotated coordinate system (x', y') with axis x' in the direction of the sun and axis y' in the direction perpendicular to the sun. In that case we deal with two variables v and w that are respectively the derivative of the surface height z along the direction of the sun and the counterpart derivative in the direction perpendicular to the sun. They are related to the slopes in the horizontal and vertical dimension of the image p and q by the following equations: 𝒗 = 𝐬𝐢𝐧(𝝋 𝟎 ) 𝒑 -𝐜𝐨𝐬(𝝋 𝟎 ) 𝒒 𝒘 = 𝐜𝐨𝐬(𝝋 𝟎 ) 𝒑 + 𝐬𝐢𝐧(𝝋 𝟎 ) 𝒒 (4) where 𝜑 " is the azimuth of the sun direction. Inversely, p and q can be recovered from v and w by: 𝒑 = 𝐬𝐢𝐧(𝝋 𝟎 ) 𝒗 + 𝐜𝐨𝐬(𝝋 𝟎 ) 𝒘 𝒒 = -𝐜𝐨𝐬 (𝝋 𝟎 )𝒗 + 𝐬𝐢𝐧 (𝝋 𝟎 )𝒘 Following the change of coordinate system, the formulae that give the local incidence, emergence, and azimuth angles as a function of the sun elevation/azimuth, and as a function of the facet slope/aspect are modified to take into account the new slope components resulting in an updated definition of R(v,w).

When present, the wavy pattern is mostly concentrated in the image w whereas there is little contribution in the image of the conjugated variable v (First row of images in Fig. 2). Consequently errors in the y' direction should be more penalized than errors in the x' direction. The integrability constraint (first regularization term of ( 1)) is split in a first term: 𝜆 A (1 -𝛿)(𝑣 -𝑧 y ) C for the direction parallel to the sun and a second term: 𝜆 A 𝛿(𝑤 -𝑧 " ) C for the direction perpendicular to the sun. The factor 𝛿≥0.5 allows to put more emphasis on the penalization of the error in the y' direction. Fig. 2 (second row) illustrates that applying a weight three times higher for the second term than for the first term (𝛿=0.75) in the integrability constraint diminishes very significantly the wavy pattern in the direction perpendicular to the sun. Note that the images shown in Fig. 2 may be stretched differently, which will not influence the discussions.

C. Photogrammetry Constraint

The photogrammetry constraint penalizes directly solution height fields that diverge too much from the initial DEM and is defined as:

𝟏 𝟐 (𝐁 • 𝐅 * 𝐳 -𝐁 • 𝐅 * 𝐳 𝟎 ) 𝟐
where B is a binary mask filtering the points where the photogrammetry method fails to find the matched points. F denotes a low-pass filter, accounting for the resolution difference between the initial and refined DEMs. In this paper we use a Gaussian function with fixed standard deviation expressed in pixels for the experiments. 𝐳 " is acquired from released DEMs for HIRISE or using ASP for CTX images.

D. Numerical Solutions

To minimize (1) with respect to v, w, z, and 𝛽, we employ the alternate minimization method. where 𝛼 A " and 𝛼 C " are the step sizes which is solved using the inexact line search algorithm [START_REF] Fletcher | Practical methods of optimization[END_REF]. For the partial derivative of R(𝐯, 𝐰) with respect v and w, numerical approximations can be obtained using: [START_REF] Soderblom | Meter-Scale 3-D models of the martian surface from combining MOC and MOLA data[END_REF] where 𝜖 is a small constant which is empirically set as 0.0001 For z, the problem becomes minimizing the cost function:

𝑬(𝐳) = 𝝀 𝟏 𝟐 •(𝐩 -𝐳 𝒙 ) 𝟐 + '𝐪 -𝐳 𝒚 ' 𝟐 " + 𝝀 𝟐 𝟐 (𝐁 • 𝐅 * 𝐳 -𝐁 • 𝐅 * 𝐳 𝟎 ) 𝟐 (11
) This Euler-Lagrange equation expresses a variational problem that admits well know solutions. Because the equation is expressed in the original image coordinate system, one needs to transform the vectors v and w back to (p, q): 𝐩 = 𝐬𝐢𝐧(𝝋 𝟎 ) 𝐯 + 𝐜𝐨𝐬(𝝋 𝟎 ) 𝐰 𝐪 = -𝐜𝐨𝐬(𝝋 𝟎 ) 𝐯 + 𝐬𝐢𝐧(𝝋 𝟎 ) 𝐰 (12) Here, we omit the deriving process which is complicated and lengthy as, interested readers may refer to Grumpe and Wöhler [START_REF] Grumpe | Recovery of elevation from estimated gradient fields constrained by digital elevation maps of lower lateral resolution[END_REF] for a similar detailed derivation. The numerical solution is given as:

𝐳 𝒏OE𝟏 = '𝐳 𝒉 𝒏 OE𝐳 𝒗 𝒏 'J 𝒄 𝟐 𝐁𝝀 𝟐 𝝀 𝟏 𝒄 𝟏 𝒄 𝟐 J𝒄 𝟐 '𝐩 𝒙 𝒏V𝟏 OE𝐪 𝒚 𝒏V𝟏 ' 𝒄 𝟐 𝐁𝝀 𝟐 𝝀 𝟏 𝐅 𝟎 𝒄 𝟏 OE𝟒 ( 13
)
where c is the cell size,

p … "OEA = oe •,žVŸ VŸ Joe •,žjŸ VŸ C¡ , q † "OEA = £ •VŸ,ž VŸ J£ •jŸ,ž VŸ C¡ , z ¥ " = z ¦, §JA " + z ¦, §OEA " , z y " = z ¦JA, § " + z ¦OEA, § " , F " = F ¦ © , § © denotes
F ¯= F ¦ © , § © JA , F °= F ¦ © , § © OEA , F ± = F ¦ © JA, § © , F ² = F ¦ © OEA, § © . 𝐅 𝒇𝟎 = 𝐅 ¦ © ³", § © ³"
denotes a filter with the same weights as F but with center weight being zero. Similarly,

𝐅 𝒇𝒉 = 𝐅 (¦ © ³", § © JA³"),(¦ © ³", § © OEA³") , 𝐅 ´y = 𝐅 (¦ © JA³", § © ³"),(¦ © OEA³", § © ³") . z … " = µ •,žVŸ Jµ •,žjŸ C¡ , z † " = µ •VŸ,ž Jµ •jŸ,ž C¡ , z C°" = z ¦, §OEC " , z C" = z ¦, §JC " , z C± " = z ¦JC, § " , z C² " = z ¦OEC
The algorithm, identified by the acronym HDEM, is summarized as follows: First, for a new 𝐳 " (𝐳 " will be used for the first iteration), 𝐳 … " and 𝐳 † " are calculated, then transformed to 𝐳 " " and 𝐳 y " in the sun coordinate system. Second, we get 𝐯 "OEA and 𝐰 "OEA by decreasing 𝐸(𝐯, 𝐰), then transform them back to 𝐩 "OEA and 𝐪 "OEA in the image coordinate. Third, a newer 𝐳 "OEA is updated by decreasing 𝐸(𝐳) and a newer 𝛽 "OEA is updated by decreasing 𝐸(𝛽). Finally, increasing the number n starts a new iteration. The algorithm will stop when a maximum iteration number is reached or when the energy function of (1) is no longer decreased for a certain number of iterations.

The proposed method contains two meta-parameters 𝜆 A and 𝜆 𝟐 to control the relative contribution of the regularization terms. As in [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF], the optimal values for 𝜆 A and 𝜆 C are those that insure the minimization of two quality criteria after convergence of the algorithm. First we have the Root Mean Square Error (RMSE) between available MOLA reference points and the corresponding points in the refined DEM. Second we define the Gradient Magnitude Similarity Deviation (GMSD) between the ortho-image and the corresponding simulated BOA reflectance image. It was noted in [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF] that parameter 𝜆 A mainly influences the GMSD index through the weighting of the image intensity constraint while parameter 𝜆 C mainly influences the RMSE index through the weighting of the photogrammetry constraint.

E. Multi-scale Analysis of Digital Elevation Models

The Isotropic Undecimated Wavelet Transform of DEMs allows us to conduct their multi-scale analysis. For that purpose we use the IUWT algorithm proposed by [Starck and Murtagh [33]] which is adapted to discrete 2D fields and is very convenient as it can be applied to a matrix of values regardless of its dimensions. Furthermore the reconstruction of the signal is straightforward and is written as the sum:

𝒄 𝑱 = 𝒄 𝟎 + ½ 𝒅 𝒋 𝑱 𝒋³𝟏
where c0 represents the smoothest version of the signal that is considered and 𝑑 § the details of the signal at decreasing spatial scales j. The J scales follow a dyadic behavior as the spatial resolution of the j th scale is 𝑟 § = 2 ÁJ § 𝑟 Á , rJ being the initial resolution of the signal. In this paper our spatial 2D fields (elevation or image intensity) are usually decomposed into J+1=9 scales, leading to a pyramidal image cube of depth 9.

III. EXPERIMENTS AND DISCUSSION

A. Validation of the HDEM Algorithm: the Columbia Hills Imaged by CTX and HiRISE

In this section we aim at validating and characterizing the products generated by the HDEM algorithm by inter-comparison and multi-scale analysis. We consider two sources for that purpose: (i) a photogrammetrically generated medium resolution CTX DEM (6 m.pixel-1) covering the region of the Columbia Hills at Gusev Crater (Mars) where the rover Spirit operated from 2004 to 2010. The hills were named for the lost Columbia space shuttle astronauts. The highest peak, Husband Hill, rises about 90 meters above the crater's floor. The DEM is accompanied by the corresponding orthorectified image. (ii) a photogrammetrically generated high resolution HiRISE DEM (1 m.pixel-1) covering a limited sector of the previous region. The DEM is also accompanied by the corresponding orthorectified image.

1) Producing the Refined DEMs

The procedure described in Section II-D is first applied to a digital elevation model and orthorectified image generated by photogrammetry from CTX imagery. These products are respectively described at URL1 * and at URL2 † . Fergason et al. (2017) produced the DEM as part of an effort to characterize topography within candidate landing ellipses for the Mars2020 mission and to aid entry, descent, and landing simulations for the mission. The DEM was constructed from a pair of CTX images (see TABLE I) using the stereo photogrammetry software SOCET SET (BAE Systems). The initial DEM (5535x9480) presents significant noise and artifacts. The best combination of meta-parameters for the retrieval is 𝜆 A =1.0e-03 and 𝜆 C =1.0e-07 based on the conjugated minimization of the cost function (1), the root mean square error between z1 and zMOLA (RMSE= 7.12 m instead of 7.20m for z0) at the available MOLA points in the scene (color-coded in Fig. 7), and the value of the GMSD indicator (0.0755) between the reconstituted image I (𝑅(𝐯, 𝐰)) and the original IL. There are other parameters used in the calculation i.e. the kernel values of the photometric model (a geographical mixture of 70 % standard dust with k L =0.40, k G =0.03, and k V =0.13, and 30% basaltic sand with k L =0.17, k G =-0.012, and k V =0.087), the geometrical conditions (see TABLE I, left column), and the Aerosol Optical Depth (AOD=0.5). Upon convergence (800 iterations), the refined DEM generated by HDEM is satisfactory as illustrated in Fig. 3 for a representative region of interest. One shaded version of the initial CTX DEM z0 (𝜃 " =59°, 𝜑 " =311°) and two shaded versions of the refined DEM z1 -one according to the initial geometrical conditions and the other with (𝜃 " =60°, 𝜑 " =180°) to reveal any possible artefact linked to photometry -are displayed. The defects and noise initially present in the DEM are removed while increasing the level of details as shown by the improved sharpness of the shaded image. The HDEM algorithm is then applied to a digital elevation model and orthorectified image generated by photogrammetry from HiRISE imagery. These products are respectively described at URL1 ‡ and at URL2 § . The HiRISE team constructed the DEM from a pair of HiRISE images (see TABLE II) using the stereo photogrammetry software SOCET SET (BAE Systems). The initial DEM (6524x10847) present a high signal to noise but some artifacts. The best combination of meta-parameters for the retrieval is 𝜆 A =1.0e-02 and 𝜆 C =1.0e-08 based on the conjugated minimization of the cost function (1), the root mean square error between z1 and zMOLA (RMSE= 3.50 m instead of 3.54m for z0) at the available MOLA points in the scene (color-coded in Fig. 7), and the value of the GMSD indicator (GMSD=0.16) between the reconstituted image I and the original IL. There are other parameters used in the calculation i.e. the kernel values of the photometric model (k L =0.40, k G =0.03, and k V =0.13, i.e. standard dust model for Mars), the geometrical conditions (see TABLE II, left column), and the Aerosol Optical Depth (AOD=0.5). Upon convergence of the HDEM algorithm (800 iterations), the refined DEM is satisfactory as illustrated in Fig. 4 for a representative region of interest. One shaded version of the initial HiRISE DEM z0 (𝜃 " =60°, 𝜑 " =299°) and two shaded versions of the refined DEM z1 (one according to the initial geometrical conditions and the other with (𝜃 " =60°, 𝜑 " =180°) are displayed. The defects initially present in the DEM are removed while increasing the level of details as shown by the improved sharpness of the shaded image.

2) Self-consistency Check

We start from the refined HiRISE DEM of superior quality zref from which:

(i) we built, thanks to the radiative transfer model of ( 2), a synthetic reflectance image Iref according to acquisition conditions identical to the HiRISE observation ( 𝜃 " =60°, 𝜑 " =299°, 𝜃=15.8°, 𝜑=73.6°),

(ii) we filter out the highest spatial frequencies to generate a degraded smooth version ž0 based on the pyramidal decomposition of the DEM zref by the IUWT algorithm on J+1=6 scales. We keep 𝒄 " and details 𝒅 § , 𝑗 = 1,2.

Then we run the HDEM algorithm until convergence (500 iterations) in order to refine ž0 based on the intensity variations of Iref. The best combination of parameters for the retrieval is 𝜆 A =1.0e-03 and 𝜆 C =1.0e-05 based on the conjugated minimization of the cost function (1), the root mean square error between z1 and zref, and the value of the GMSD indicator between the reconstituted image I1 and Iref. The other parameters used in the calculation are identical to the previous experiment. Finally we compare the refined DEM z1 with zref for assessing the capability of the algorithm to generate the missing details by photoclinometry. First we calculate a series of indicators: the root mean square error between z1 and zref (0.08m), between the corresponding slopes (0.028) and the GMSD value (0.01) between I1 and Iref. Second we visualize portions of 3 topographical profiles respectively extracted from zref, ž0, and z1 for a region of interest (Fig. 5). We deduce from the experiment that the HDEM algorithm is indeed able to reconstruct the original topography of the scene by exploiting the intensity variation of the image, although not fully.

The corresponding slope profiles v and w are also extracted and displayed in Fig. 6. The slope vref is very well reproduced by v1 whereas the slope wref is only partially reproduced by w1. We deduce from these experiments that the retrieval of slopes by the HDEM algorithm is more sensitive in the direction aligned with the sun than in the direction orthogonal to the latter. Consequently the performance of the HDEM algorithm operating on a single image IL of the scene is expected to be satisfactory but with a limited smoothing in the direction orthogonal the sun azimuth ( [30-210°] in this HiRISE case).

3) Comparison between CTX and HiRISE

We now perform the comparison between three products: the original CTX DEM z0, the refined version derived by the HDEM algorithm z1, and, as a reference, the refined HiRISE DEM zref downscaled at 6 m.pixel-1. We start the interpretation with the topographic profiles extracted along a transect aligned with the CTX illumination direction (Fig. 7, brown line in the web version of the article). In terms of altitude (Fig. 8, top) we see that the CTX original and refined profiles share a common baseline that can be offset by up to 5 meters compared to the HiRISE baseline. This comes from the initial photogrammetric processing that implied different bundle adjustments should the CTX or the HiRISE product be considered. The refined DEM inherits the topographical properties of the initial product at large scales (several hundreds of meters or more for CTX). Nevertheless it can be noticed that at smaller scales the refined version of the CTX DEM display topographical structures that are very similar in position, shape, and amplitude compared to the reference HiRISE DEM whereas the initial DEM generally lacks such structures. The comparison of the three altimetric profiles (z0, z1, and zref) along a transect orthogonal to the CTX illumination direction (Fig. 7, orange line in the web version of the article) allows to draw a similar conclusion (Fig. 9, top) even though in this case the baseline is nearly common to all. In terms of slopes (Fig. 8, bottom), in agreement with the results of the self-consistency check, the slopes of the reference profile are much better reproduced by the refined DEM than the initial one for which high spatial frequencies are totally missing whereas they are present in the refined signal and conform to the reference. Nevertheless the amplitude of the signal is underestimated especially for the largest slopes. This is accentuated in the profile corresponding to the transect orthogonal to the CTX illumination direction (Fig. 9, bottom) likely because of the smoothing effect as noted already in the self-consistency check.

Fig. 8. First series of topographical profiles extracted from the original CTX DEM z0, the refined version of it z1, and, as a reference zref, the refined HiRISE DEM downscaled at 6m/pixel for a transect aligned with the CTX illumination direction (Fig. 7, brown line). Top: height profiles in meters. Bottom: profiles of absolute slope in degrees.

In addition we present a comparison of topographical profiles extracted from the original CTX DEM (z0), the refined version of it (z1), and, as a reference, the refined HiRISE DEM (zref) accompanied with available MOLA points along track # 19229 (Fig. 7). One can see that the baselines of the CTX profiles differ from HiRISE and MOLA by up to 10 meters (the latter two profiles being in better agreement). Thus we confirm that the refined DEM inherits the topographical properties of the initial product at large scales -not necessarily in agreement with the MOLA points -but restitutes more faithfully the topography at smaller scales. In summary, the HDEM algorithm captures satisfactorily topographical structures that can be seen in the orthorectified image but are absent or poorly represented in the photogrammetrically generated DEM. The enhancement is even more perceptible in terms of slopes. For both quantities (height and slope), the improvement is scale dependent in the sense that it rather concerns the high spatial frequencies. In the following section we study this dependency quantitatively.

4) Characteristic Spatial Scales Controlled by the Photoclinometric Process

In this subsection we carry out a multi-scale analysis of the CTX and HiRISE products (DEMs and associated orthorectified images) with the IUWT algorithm at J+1=9 levels (Section II-E). On this basis we conduct an assessment of the functioning of the HDEM algorithm at different spatial scales when refining the DEM generated by photogrammetry. For that purpose we reconstruct a series of signals from the decomposition by considering a decreasing number of spatial frequencies from the highest end of the spectrum: 𝑱′ = 𝑱, . . . , 𝟎:

𝒄 𝑱 ′ = 𝒄 𝟎 + Ä 𝒅 𝒋 𝑱 ′

𝒋³𝟏

The first signal is at full resolution whereas the following signals are increasingly smoothed. Note that ž0 is initially (J'=J) a smoothed version of the photogrammetric DEM to eliminate any spurious artifact (ž0=g*z0 , g being a Gaussian filter). First we look at the potential for improving the initial DEM ž0 using the left image IL depending on spatial scale J'. We estimate this quantity by calculating the GMSD indicator between the filtered version of the real orthorectified image and the image reconstituted by our radiative transfer model from the filtered version of ž0 (GMSD(IL,RBOA0), upper left graph of Fig. 11). Then we start quantifying what was actually modified by HDEM by calculating the GMSD indicator between the images reconstituted respectively from the filtered version of ž0 and z1 (GMSD(RBOA1,RBOA0), upper right graph of Fig. 11). We then proceed by a comparison between the initial DEM ž0 and the refined DEM z1 in terms of altitudes (lower left graph) and in terms of slopes (lower right graph) for different levels of filtering. In both cases we calculate the root mean square deviation (RMSD) between the two products: (18) where N is the total number of valid pixels in the products, s1 and š0 stand for absolute slopes corresponding to z1 and ž0 respectively. The J' index is omitted for clarity. For the CTX experiment, we see that the filtered version of images IL and RBOA0 are the most similar when the details from j=5 to j=8 are removed (Fig. 11, upper left graph). They are increasingly different when adding back progressively these components. Consequently there is a potential for increasing the level of details of the DEM at spatial scales ≈ 6 to 100 m (i.e. 1 to 17 pixels). Actually the potential for improvement was mostly exploited by the HDEM algorithm for the CTX experiment as the GMSD(RBOA1,RBOA0) curve (Fig. 11, upper right graph) is similar to GMSD(IL,RBOA0) even though all the information from IL could not be picked up: GMSD(RBOA1,RBOA0) < GMSD(IL,RBOA0) at all scales. The synthetic image RBOA1 retains some properties of RBOA0 that are not included in IL. It is confirmed that the HDEM algorithm operates at scales approximately 6 to 100m. For the CTX experiment we note that the RMSD between z1 and ž0 (Fig. 11, lower left graph) at full resolution is 2.2m. The characteristic spatial scale of the photoclinometric processcorresponding to the decrease of the quantity of interest (RMSD) by a factor e=2.72 -is approximately 160-170m (27 pixels). The relative RMSD is defined as the RMSD between the slopes of z1 and the slopes of ž0 divided by the mean of the slopes ž0. At full scale it is 77% then declining sharply to reach a threshold of 10 % at 250m (Fig. 11, lower right graph). The HDEM algorithm has thus a significant effect on slopes within a large range of spatial scales (up to 42 times the pixel size).

𝑹𝑴𝑺𝑫(𝐳 𝟏 , 𝐳 AE 𝟎 ) = Ç 𝟏 (𝑵 -𝟏) ½(𝐳 𝟏 (𝒏) -𝐳 AE 𝟎 (𝒏)) 𝟐
The same analysis is then applied to the DEM and orthorectified image generated from HiRISE imagery of the Columbia hills. Note that this experiment HiRISE_1513_1655 was performed on a region of interest 2048x2048. We see that the filtered version of images IL and RBOA0 are the most similar when the details from j=4 to j=8 are removed (Fig. 11, upper left graph). They are increasingly different when adding back progressively these components. Consequently there is a potential for increasing the level of details of the DEM at spatial scales approximately 1 to 34 m (i.e. 1 to 34 pixels). Actually the potential for improvement was mostly exploited by the HDEM algorithm in both experiments as the GMSD(RBOA1,RBOA0) curve converges at scale 32m (Fig. 11, upper right graph). For the HiRISE experiment we note that the RMSD between z1 and ž0 (Fig. 11, lower left graph) at full resolution is respectively 0.62m. The characteristic spatial scale of the photoclinometric process -corresponding to the decrease of the quantity of interest by a factor e=2.72 -is approximately 23m (23 pixels). At full scale the relative RMSD of the slopes is 66% then declining sharply to reach a threshold of 10 % at 32m. The HDEM algorithm has thus a significant effect on slopes within a large range of spatial scales (up to 32 times the pixel size).

From both the CTX and HiRISE experiments we conclude that our photoclinometric scheme HDEM exploits in a very satisfactory manner the intensity variations of the left image in order to refine the photogrammetrically generated initial DEM. We note that the HDEM algorithm operates at spatial scales approximately 1 to 30 times the original pixel size. In this range the effect on the height field and the slopes depends on the considered scale, being strong for the smallest scales and decreasing rapidly for increasingly coarser scales. If there is a potential for increasing the level of details of the DEM at spatial scales larger than the characteristic range of the HDEM method (approximately 1 to 30 pixels) then it is necessary to adopt a coarse-to-fine strategy.

B. Coarse-to-fine Strategy for Processing Challenging Scenes: the Russell Dune Imaged by HiRISE

Russell is a crater on Mars, located in the Noachis quadrangle at 54.9° south latitude and 347.6° west longitude, and gullies are found on some dunes. The initial HiRISE DEM z0 and corresponding orthorectified image IL used in the photoclinometric process are described respectively in URL3 * and in URL4 † . See The initial DEM is of poor quality with a lot of artifacts: stripes, seams, etc., challenging the standard version of our HDEM algorithm (see Fig. 12, left image). We start the experiment by calculating the GMSD indicator between the low pass filtered version of IL and the reflectance image (RBOA0) reconstituted from the low pass filtered version of z0: GMSD(IL(J'),RBOA0(J')). We see that beyond spatial scale J'=2 (65 m) the filtered left image does not convey more information on the topography than the image built on the filtered DEM z0. From experiments described in Section III-A (Characteristic spatial scales controlled by the photoclinometric process) we know that the HDEM algorithm is only able to add spatial details over 5 levels of spatial scales, at most 30 pixels. Consequently, in an attempt to improve the processing of these HiRISE data for which the initial DEM is low in quality, we propose a coarse-to-fine strategy with HDEM involving two levels of processing. In the first level, the initial DEM z0 and the left image IL are both downscaled to coarse version z0_c and IL_c (J'=6, 4 m.pixel -1 ), and the HDEM algorithm is adopted to produce a refined DEM z1_c from ž0=g*z0_c and IL_c, as shown in Fig. 14 (g* stands for a Gaussian convolution with sigma equal to 15 pixels, i.e. ž0 is at 64 m.pixel -1 ). In the second level, the z1_c and the original left image IL are put into the HDEM algorithm to produce the final refined DEM z1, as shown in Fig. 15. From the shaded version of the final DEM in Fig. 15, we see that the convergence of HDEM on these data is very satisfactory leading to a new DEM with few artifacts, sharper than z1_c and much cleaner and detailed than the initial DEM z0. The evolution of the height field from the filtered version of the initial DEM ž0 to the refined version at coarse resolution z1_c, and then from the latter to the refined version at fine resolution z1 is illustrated in Fig. 16 by the extraction of topographic profiles. The capability of the HDEM algorithm to reconstitute topography over 4 spatial scales (from 64 m.pixel -1 to 4m.pixel -1 , green to red curve) from an image is clearly demonstrated after the first level of processing, both in terms of altitude (Fig. 16 top) and absolute slope (Fig. 16 bottom). The gain is still significant (4m.pixel -1 to 1m.pixel -1 , red to blue curve) after the second level of processing especially in terms of slope.

C. Investigating Planetary Processes by Characterizing the Topography at Small Scales

In this subsection we present topographical products derived from DEMs aimed at understanding planetary processes. These products are based on the calculation of the local roughness after the definition of Wilson et al. (2007): the difference between the maximum and minimum values of the height field within a user defined n × n rectangular neighborhood surrounding the central pixel with n=3. Height fields can be decomposed into a series of spatial components d j at different scales j thanks to the IUWT algorithm (Section II-E). For each spatial component, the local roughness can be calculated leading to a multiscale topographic characterization of the terrain. By selecting scales that discriminate different features, one can build informative false color compositions. In addition to the planetary interest, the previous topographical products can further demonstrate the added value of our photoclinometric HDEM method. In the following, we consider the "Columbia hills" imaged by HiRISE (see Section III-A, "Producing the refined DEMs") and we perform the multi-scale topographical analysis of three DEM versions for experiment #1513 (initial, refined at 1m.pixel -1 , and refined at 0.25m.pixel -1 ). The latter product comes from the operation of the HDEM algorithm on an orthorectified HiRISE image available at 0.25m.pixel -1 using the refined DEM at 1m.pixel -1 as the initial guess z0. In this case no Gaussian filtering is necessary because of the high quality of z0 (ž0=z0). For this experiment we calculate the roughness of the scene at three spatial scales (1x, 2x and 4x) respectively attributed to the red, green, and blue channels of a false color composition. The latter can be superimposed on the HiRISE orthorectified image (Fig. 17). Variations of hues reveal different characteristics of the scene. We note that the level of useful topographical information is clearly higher in the refined versions (middle and lower compositions) than in the initial version (upper composition) that shows banded artefacts. The refined DEM versions clearly reveal meaningful planetary features. First reddish colors correspond to small wavy patterns affecting sand deposits inside big craters or in the hills. They turn to yellow, and then to white if their amplitude increases. Second reddish hues also form digital aureoles around the biggest craters where they reveal the presence of small blocks ejected by the impact. Third bluish colors indicate the deepest structures of the scene (e.g. highly dissected portions of the hills) or pronounced slope discontinuities (e.g. crater ramparts).

IV. CONCLUSIONS AND PERSPECTIVES

In this paper the mathematical and physical basis of a photoclinometric scheme aimed at refining DEMs generated by photogrammetry from in-orbit Martian imagery are first presented. An evolution of the HDEM algorithm for solving the underlying optimization problem is then introduced in which the sun-oriented coordinate system is used in the integrability constraint. In the sun coordinate system, we propose a weighting factor to pose more penalization on the direction normal to the sun azimuth. The concept of retaining only the largest reliable spatial components of the initial DEM is also defined and translated into a new set of numerical equations for updating the height field. The principal goal of changing the photoclinometric scheme is to reduce the sensitivity of the algorithm to artifacts and noise present in the initial DEM. Experiments are conducted to study and validate the refined DEMs generated by HDEM. First we consider CTX and HiRISE data covering the Columbia Hills region for which the photogrammetrically generated initial DEMs is already of (quite) good quality. The products generated by the improved version of the HDEM algorithm are cross-validated. We show that the algorithm is able to recover the shape and amplitude of topographical structures that are absent or poorly represented in the photogrammetrically generated DEM thanks to the intensity information contained in the associated ortho-image. Nevertheless the recovering is better in the direction parallel to the sun than in the orthogonal direction. It also occurs on a limited range of spatial frequencies. An important contribution of this paper is the careful assessment of the modifications brought by the HDEM algorithm to the height field and its derivatives compared to the initial DEM at different spatial scales. We note that a single run of the HDEM algorithm operates at horizontal spatial scales approximately 1 to 17 times the original pixel size for CTX imagery (6m/pixel) and 1 to 30 times for HiRISE imagery upon convergence of the algorithm (»500 iterations). In this range the effect on the height field at full scale expressed by the root mean squared deviation is of the order of 2m for the CTX derived DEM and 0.7 m for the HiRISE derived DEM. Relative slope deviations between the original and refined DEMs are of the order of 60-70% at full scale. Then, with an experiment conducted on HiRISE data covering the Russell dune, we consider situations for which the photogrammetrically generated initial DEM is of poor quality and only the spatial components of the height field at the largest scales are reliable. For that purpose we introduce a successful coarse-to-fine strategy based on an iterative use of HDEM at different spatial scales. Finally we open the possibility to create new noise-free topographical products aimed at characterizing small topographic scales thus contributing to the understanding of planetary processes.

Our method provides a very significant gain of precision for the quantitative topographic characterization of Mars at small horizontal scales (sub-metric to hectometric) with access to structures whose height can be measured relatively with subtle precision (down to 15 cm in practice with HiRISE). It is particularly suitable for the processing of images produced by HiRISE for which the gain is maximized. However, the partial reconstruction of the slopes in the direction perpendicular to the solar azimuth is an intrinsic, but moderate limitation of the photoclinometric approach when the inversion of the image model is performed on a single image (one geometric configuration) as it is the case in our experiments. The only way to overcome this non-critical limitation -as well as the question of a spatially variable photometric model -is an inversion of multiple images with different acquisition geometries. This will motivate new investigations.

Fig. 1 .

 1 Fig. 1. The flowchart of the proposed method.

  𝜌 A and 𝜌 C are defined as 𝜌 A (𝑠, 𝑣, 𝑤) = A Ch ∫ 𝜌(𝑠 ′ , 𝑠, 𝑣, 𝑤) Ω V 𝑑𝑠 ′ and 𝜌 C (𝑠 " , 𝑣, 𝑤) = A Ch ∫ 𝜌(𝑠 " , 𝑠 ′ , 𝑣, 𝑤) Ω j

Fig. 2 .

 2 Fig. 2. Images depicting the two slope components (v->left column, w->right column) estimated by the photoclinomeric algorithm for a small region of interest from a CTX image. First row: the algorithm is run with d=0.5. Second row: the algorithm is run with d=0.75.

  For w and v, the cost function can be written: 𝑬(𝐯, 𝐰) = 𝟏 𝟐 (𝜷 • 𝑹(𝐯, 𝐰) -𝐈) 𝟐 + 𝝀 𝟏 𝟐 [(𝟏 -𝜹)(𝐯 -𝐳 𝒗 ) 𝟐 + 𝜹(𝐰 -𝐳 𝒘 ) 𝟐 ] 𝐳 y and 𝐳 " are the partial derivative of z in the sun coordinates and we assume 𝐳 𝒗𝒘 = 𝐳 𝒘𝒗 . They derive from the derivatives 𝐳 … and 𝐳 † that are calculated in the image coordinate system by a linear transformation: 𝐳 𝐯 = 𝒔𝒊𝒏(𝝋 𝟎 ) 𝐳 𝐱 -𝒄𝒐𝒔(𝝋 𝟎 ) 𝐳 𝐲 𝐳 𝐰 = 𝒄𝒐𝒔(𝝋 𝟎 ) 𝐳 𝐱 + 𝒔𝒊𝒏(𝝋 𝟎 ) 𝐳 𝐲 Then we can use gradient descent for updating: 𝐯 𝒏OE𝟏 = 𝐯 𝒏 -𝜶 𝟏 𝒏 𝝏𝑬(𝐯,𝐰) 𝝏𝐯 , 𝐰 𝒏OE𝟏 = 𝐰 𝒏 -𝜶 𝟐 𝒏 𝝏𝑬(𝐯,𝐰) 𝝏𝐰 where 𝝏𝑬(𝐯, 𝐰) 𝝏𝐯 = 𝐯 𝒏 -𝜶 𝟏 𝒏 Ž(𝜷 𝒏 • 𝑹(𝐯 𝒏 , 𝐰 𝒏 ) -𝐈)𝜷 𝒏 𝝏𝑹(𝐯 𝒏 , 𝐰 𝒏 ) 𝝏𝐯 + 𝝀 𝟏 (𝟏 -𝜹)(𝐯 𝒏 -𝐳 𝒗 )• 𝝏𝑬(𝐯,𝐰) 𝝏𝐰 = 𝐰 𝒏 -𝜶 𝟐 𝒏 •(𝜷 𝒏 • 𝑹(𝐯 𝒏OE𝟏 , 𝐰 𝒏 ) -𝐈)𝜷 𝒏 𝝏𝑹'𝐯 𝒏V𝟏 ,𝐰 𝒏 ' 𝝏𝐰 + 𝝀 𝟏 𝜹(𝐰 𝒏 -𝐳 𝒘 )"

  , § " . For 𝛽 , by setting the derivative of the energy function E(𝛽) = A C (𝛽 • 𝑅(v, w) -I) C with respect to 𝛽 to zero, we get the closed form of the solution as: 𝜷 𝒏OE𝟏 = W𝑹'𝐯 𝒏V𝟏 ,𝐰 𝒏V𝟏 ' 𝑻 𝐈X 𝑹'𝐯 𝒏V𝟏 ,𝐰 𝒏V𝟏 ' 𝑻 𝑹'𝐯 𝒏V𝟏 ,𝐰 𝒏V𝟏 '

Fig. 3 .

 3 Fig. 3. Shaded representation of the original and refined DEMs restricted to a region of interest. Left image: initial CTX DEM with (𝜃 " =59°, 𝜑 " =311°). Middle image: refined CTX DEM with (𝜃 " =59°, 𝜑 " =311°). Right image: refined CTX DEM with (𝜃 " =59°, 𝜑 " =180°).

Fig. 4 .

 4 Fig. 4. Shaded representation of the original and refined DEMs restricted to a region of interest. Left image: initial HiRISE DEM with (𝜃 " =60°, 𝜑 " =299°). Middle image: refined HiRISE DEM with (𝜃 " =60°, 𝜑 " =299°). Right image: refined HiRISE DEM with (𝜃 " =60°, 𝜑 " =180°).

Fig. 5 .

 5 Fig. 5. Topographic profiles extracted from the reference DEM zref, the smooth version of the latter ž0, and the refined DEM z1 for different short portions of the same vertical (i.e. North-South) line through a region of interest.

Fig. 6 .

 6 Fig. 6. slope profiles extracted from the derivation of the reference DEM zref, of the smooth version ž0, and of the refined DEM z1 for the same vertical (i.e. North-South) line through a region of interest. Top: slope in the direction of the sun (vref, v0, and v1). Bottom: slope in the direction perpendicular to the sun (wref, w0, and w1).

Fig. 7 .

 7 Fig. 7. Shaded representation of the refined HiRISE DEM (bounding box) superimposed on the shaded representation of the refined CTX DEM covering a larger area of the Columbia Hills region. MOLA reference points are color-coded (in the web version of the article) depending on their orbit number (track 19229 is in indicated with bigger symbols). Two lines along which different topographical profiles are extracted are also represented: (i) a transect aligned with the CTX illumination direction (brown in the web version), (ii) a transect orthogonal to the CTX illumination direction across eroded terrains (orange in the web version).

Fig. 9 .

 9 Fig.9. Second series of topographical profiles extracted from the original CTX DEM z0, the refined version of it z1, and, as a reference zref, the refined HiRISE DEM downscaled at 6m/pixel for a transect orthogonal to the CTX illumination direction (Fig.7, orange line). Top: height profiles in meters. Bottom: profiles of absolute slope in degrees.

Figure 10 :

 10 Figure 10: Series of topographical profiles extracted from the original CTX DEM, the refined version of it, and, as a reference, the refined HiRISE DEM downscaled at 6m/pixel and the MOLA points of track 19229 (black crosses).

  𝟏 , 𝒔 AE 𝟎 ) = É 𝟏 (𝟐𝑵J𝟏) WÄ (𝐯 𝟏 (𝒏) -𝐯 𝟎 (𝒏)) 𝟐 𝑵 𝒏³𝟏 + Ä (𝐰 𝟏 (𝒏) -𝐰 𝟎 (𝒏))

Fig. 11 .

 11 Fig. 11. Series of graphs illustrating the potential and the actual functioning of the HDEM photoclinometric process at different spatial scales. Two experiments with CTX and HiRISE imagery are considered (see main text for a description). Upper left: GMSD indicator between the filtered version of the real orthorectified image and the image reconstituted by our radiative transfer model from the filtered version of ž0. Upper right: GMSD indicator between the images reconstituted from the filtered version of respectively ž0 and z1. Lower left: root mean square deviation between the initial DEM ž0 and the refined version z1 as a function of the scale of filtering. Lower right: relative root mean square deviation between the slopes of z1 and the slopes of ž0.

Fig. 12 .

 12 Fig. 12. Shaded representation of the initial HiRISE DEM (left) and corresponding orthorectified image for a region of interest 2048x2048 in dimension.

Fig. 13 .

 13 Fig. 13. Evolution of the GMSD indicator between the filtered HiRISE image and the reflectance image reconstituted from the filtered version of z0 as a function of the spatial scale.

Fig. 14 .

 14 Fig. 14. Inputs ž0=g*z0_c (left shaded image), IL_c (middle image) and result z1_c (right shaded image) obtained by running the HDEM algorithm at coarse resolution.

Fig. 15 .

 15 Fig. 15. Inputs z1_c (left shaded image), original left image IL(middle image) and result z1 (right shaded image) obtained by running the HDEM algorithm at fine resolution.

Fig. 16 :

 16 Fig.16: topographical profiles extracted from the filtered version of the original HiRISE DEM (ž0), the refined version of it at coarse resolution (z_1c), and the refined version of it at fine resolution (z1) for a transect across some gullies (Fig.15). Top: height profiles in meters. Bottom: profiles of absolute slope in degrees.

Fig. 17 .

 17 Fig. 17. False color (in the web version of the article) compositions of the Columbia hills region indicative of the roughness of the scene at three consecutive spatial scales. The compositions are superimposed with transparency on the HiRISE orthorectified image. Top: composition derived from the initial HiRISE DEM at 1 m.pixel-1. Middle: composition derived from the refined HiRISE DEM at 1 m.pixel -1 . Bottom: composition derived from a refined HiRISE DEM at 0.25m.pixel -1 .

  the center weight of filter F , 𝑐 A and 𝑐 C are calculated as:

	𝒄 𝟏 =	𝐳 𝒙 𝒏 ª	𝐅 𝒍 𝐳 𝟐𝒓 𝒏 𝟐𝒄	J	𝐅 𝒓 𝐳 𝟐𝒍 𝒏 𝟐𝒄	OE𝐅 𝒇𝒉 * 𝐳 𝒙 𝒏 -OE𝐳 𝒚 𝒏 ª (𝐳 𝒙 𝒏 ) 𝟐 OE'𝐳 𝐲 𝐅 𝒖 𝐳 𝟐𝒅 𝒏 𝟐𝒄 𝐧 ' 𝟐	J 𝐅 𝒅 𝐳 𝟐𝒖 𝒏 𝟐𝒄	OE𝐅 𝒇𝒗 * 𝐳 𝒚 𝒏 -
	𝒄 𝟐 = 𝐅 𝒇𝟎 * 𝐳 𝒏 -𝐅 * 𝐳 𝟎		(14)
	where

TABLE I IMAGE

 I INFORMATION OF THE CTX PAIR (EXPERIMENT CTX F21_043841).

			Sensor	76°	80°
			elevation
			Sun azimuth	311°	310°
			angle
			Phase angle	69°	54°
			Used	Upper-left corner, lat: -15.15 °, E lon: °;
			Region	Lower-right corner, lat: -14.17 °, E lon: 175.73 °.
		Left image	Right image
	Image ID	F21_043841_1654_XN_	F21_043907_1652_XN_
		14S184W	14S184W
	Sun	31°	29°
	elevation

TABLE II IMAGE

 II INFORMATION OF THE HIRISE COLUMBIA HILLS CASE ‡ https://www.uahirise.org/dtm/dtm.php?ID=PSP_001513_1655 (accessed 05.24.2018).

		(EXPERIMENT HIRISE #1513).
		Left image	Right image
	Image ID	PSP_001513_1655	PSP_001777_1650
	Sun elevation	30°	29°
	Sensor elevation	74.2°	86°
	Sun azimuth angle 299°	294°
	Phase angle	73.6°	57.3°
	Used Region	Upper-left corner, lat: -14.405°, E lon: 175.429°;
		Lower-right corner, lat: -14.779°, E lon: 175.568°.

§ https://www.uahirise.org/PSP_001513_1655 (accessed 05.24.2018).

  TABLE III for some basic image information.

	GMSD (RBOA_1,RBOA_0)	https://www.uahirise.org/dtm/dtm.php?ID=PSP_007018_1255 (accessed 05/24/2018) 0.25 0.05 0.00 0.20 0.15 0.10

* † https://www.uahirise.org/PSP_007018_1255 (accessed 05/24/2018)

TABLE III IMAGE

 III INFORMATION OF THE HIRISE RUSSELL DUNE CASE (EXPERIMENT HIRISE #7018).

		Left image	Right image
	Image ID	PSP_007018_1255	PSP_007229_1255
	Sun elevation	17°	13°
	Sensor elevation	86.6°	84.3°
	Sun azimuth angle	316°	314°
	Phase angle	75.6°	73.8°
	Used Region	Upper-left corner, lat: -54.204°, E lon: 12.877°;
		Lower-right corner, lat: -54.385°, E lon: 13.058 °.