Seamless 3 D Image Mapping and Mosaicing of Valles Marineris on Mars Using Orbital HRSC Stereo and Panchromatic Images

: A seamless mosaic has been constructed including a 3 D terrain model at 50 m grid-spacing and a corresponding terrain-corrected orthoimage at 12.5 m using a novel approach applied to ESA Mars Express High Resolution Stereo Camera orbital (HRSC) images of Mars. This method consists of blending and harmonising 3 D models and normalising reﬂectance to a global albedo map. Eleven HRSC image sets were processed to Digital Terrain Models (DTM) based on an opensource stereo photogrammetric package called CASP-GO and merged with 71 published DTMs from the HRSC team. In order to achieve high quality and complete DTM coverage, a new method was developed to combine data derived from different stereo matching approaches to achieve a uniform outcome. This new approach was developed for high-accuracy data fusion of different DTMs at dissimilar grid-spacing and provenance which employs joint 3 D and image co-registration, and B-spline ﬁtting against the global Mars Orbiter Laser Altimeter (MOLA) standard reference. Each HRSC strip is normalised against a global albedo map to ensure that the very different lighting conditions could be corrected and resulting in a tiled set of seamless mosaics. The ﬁnal 3 D terrain model is compared against the MOLA height reference and the results shown of this intercomparison both in altitude and planum. Visualisation and access mechanisms to the ﬁnal open access products are described.


Introduction
Valles Marineris on Mars are the largest system of canyons in the Solar System. The system is more than 4,000 km long, 200 km wide, and up to 10 km deep with respect to the surrounding plateaus. Because of their sheer size Valles Marineris have been a subject of growing number of studies since the Viking era and we now know that they have been subject to a wide variety of processes throughout their history. Estimated to be approximately 3.5 Ga old [1], the superstructure is related to tectonic rifting of the crust related in part to the loading induced by the Tharsis Volcanic province [1][2][3][4]. This overall tectonic structure has been modified by outflows [5][6][7][8], fluvial processes [7,[9][10][11] (delta, melas), volcanic processes [12][13][14], aeolian processes [15][16][17] and it also thought to contain evidence of glacial modification [18,19]. Its extreme topographic relief influences the atmospheric circulation and creates weather systems unique to the canyon [20,21]. Vast sediment accumulations have formed on its floor, whose geochemical signature attests

Datasets
In this work, we use the MOLA areoid (https://pds-geosciences.wustl.edu/missions/ mgs/megdr.html (accessed on 21 March 2021)) data as our baseline reference system (both spatially and vertically). MOLA (from 1997 to 2001) was one of the five instruments onboard MGS that measured the time-of-flight of laser pulses at 1064 nm, over the distance from the spacecraft to the Martian surface. Individual MOLA laser tracks were interpolated and extrapolated, and corrected, by the MOLA team to produce a global DTM of Mars at 463 m grid spacing [37], which is available through the NASA Planetary Database System (PDS; https://pds-geosciences.wustl.edu/missions/mgs/megdr.html (accessed on 21 March 2021)).
The HRSC stereo imaging instrument onboard the ESA's Mars Express (MeX) spacecraft is still orbiting Mars having acquired its first images in January 2004. HRSC was designed for photogrammetric mapping of the topography of the Martian surface via continuous acquisition of stereo push-broom images from multiple angles on a single orbital pass. HRSC captures multi-angular and multi-colour images mostly at 12.5-100 m pixelscale over various swath-widths, providing near global (~98%) coverage of the Martian surface at 12.5-100 m/pixel. The DLR processed HRSC level 4 DTMs cover as of the time of submission~50% of the Martian surface. In this work, we use the HRSC level 2 stereo images for UCL level 4 DTM processing and the DLR HRSC level 4 DTMs (version 50 and higher) [42] for co-registration, co-alignment with MOLA, and mosaicing. It should be noted that only the nadir image is captured at 12.5 m with the off-nadir stereo channels usually captured 2 times coarser at 25 m.

Previous Work
In [39], we previously compared the quality and performance of CTX DTMs produced using different existing planetary 3 D reconstruction software and introduced the CASP-GO system, which is based on the open-source NASA Ames Stereo Pipeline (ASP) [43], tie-point based multi-resolution image co-registration [44], the Gotcha [45] sub-pixel refinement method, and other optimisation algorithms. The CASP-GO system produces high quality DTMs in terms of improved completeness and accuracy, as well as minimising matching and interpolation artefacts.
For HRSC data, DLR generate HRSC single-strip DTMs [42] in a sinusoidal projection system at 50-150 m grid-spacing for level-4 products. More recently, DLR have developed a level-5 processing chain to produce a large (>100 strip) HRSC  DLR processing of the raw HRSC level-2 data includes radiometric de-calibration, noise reduction, image matching, geo-referencing, photogrammetric processing and alongtrack Bundle Adjustment (BA). The final products are labelled as "Level-4 v50+" when they reach a satisfactory level of quality [42]. The co-registration of HRSC DTM and MOLA DTM is checked using average height differences and usually reveal residual horizontal offsets smaller than the grid spacing of the HRSC DTM and residual vertical offsets below 10 m [42]. In addition, bundle block adjustment and joint interpolation of multi-scale 3 D point data sets are used for DTM mosaicing. For the ORI mosaic, [46] describes the standard DLR method using single strip scaling (to calibrated radiances), brightness normalisation according to Lambert reflection, and local contrast adjustment and brightness adjustment to the MGS's TES albedo [39,41].
Another related work on wide area HRSC DTM/ORI production and mosaicing is [47], in which the authors demonstrated the production of a HRSC DTM and ORI mosaic for the south polar region using the DLR/NASA VICAR (https://www-mipl.jpl.nasa.gov/vicar_ os/v1.0/vicar-docs/VICAR_guide_1.0.pdf (accessed on 21 March 2021)) based pipeline and brightness correction through a high altitude HRSC reference image as well as the use of the 128 m/pixel MOLA DTM to provide height normalisation. In this previous work, individual HRSC orbital strips were separately bundle-adjusted by FUB due to little coverage of DLR HRSC products [47].

Data Preparation
For HRSC, 71 DLR processed HRSC level 4 DTMs and ORIs were available until April 2020. These include 57 with DTMs at 50 m/pixel and ORIs at 12.5 m/pixel; 6 with DTMs at 100 m/pixel and ORIs at 25 m/pixel and 8 with DTMs at 150 m/pixel and ORIs at 50 m/pixel. A list of the image IDs with download links are provided as Supplementary Material. The 71 DLR HRSC level 4 DTMs cover most of the Valles Marineris region but leave several significant gaps. Therefore, we checked all available HRSC level 2 images and selected a set of 11 HRSC level 2 stereo pairs at 12.5 m/pixel to be processed to DTMs and mosaiced at 50 m/pixel. The corresponding level 2 HRSC orbit IDs are h2112_0000, h2123_0000, h2134_0000, h2211_0000, h2145_0000, h2178_0000, h1896_0000, h0438_0000, h3195_0000, h2182_0000, h2402_0001. The image footprints of the HRSC DLR products and UCL products are shown in Figure 1.

The ASP and CASP-GO Auto-DTM Generation System
The Ames Stereo Pipeline (hereafter referred to as ASP) software was developed by the Intelligent Robotics Group at the NASA Ames Research Center, Mountain View, CA [43]. ASP is a suite of open-source automated geodesy and stereo-photogrammetry tools designed for processing planetary imagery captured either from planetary (robotic) orbiters or landers or from Earth observing satellites. The CASP-GO pipeline is based on the ASP framework with specific enhancements, dealing with matching artefacts, gaps, and co-alignment issues, from in-house software previously developed at the Imaging Group at UCL. CASP-GO uses the ASP I/O interface and takes USGS Integrated Software for Imagers and Spectrometers (USGS-ISIS) formatted stereo images (see Figure 2; see explanation of the USGS-ISIS Cube format in https://isis.astrogeology.usgs.gov/documents/LogicalCubeFormatGuide/LogicalCubeFormatGuide.html (accessed on 21 March 2021)) and a reference base map as inputs. Converting from raw HRSC PDS data into the USGS-ISIS formatted Cube includes the process of format conversion, image pre-denoising, SPICE (see https://naif.jpl.nasa.gov/naif/spiceconcept.html (accessed on 21 March 2021)) kernel initialisation, and map projection processes (see https://isis.astrogeology.usgs.gov/Application/index.html#Mars_Express). The core matching algorithm uses a combination of cross-correlation and a 5 th generation of adaptive least squares correlation (ALSC) and region growing matcher called Gotcha [45], which provides accurate and robust sub-pixel conjugate points.
The complete CASP-GO workflow (see Figure 2) for DTM production includes the following 10 steps which are shown in the aforementioned figure and described as follows: (a) ASP image pre-processing including image enhancement (e.g., Laplacian of Gaussian filtering and image normalisation) and stereo rectification; (b) ASP disparity map initialisation using coarse cross-correlation; (c) UCL fast Maximum Likelihood (f-ML) matching and the construction of a "float" initial disparity map; (d) ASP Bayes Expectation Maximisation (EM) weighted affine adaptive sub-pixel cross-correlation; (e) UCL re-defined outlier rejection and gap erosion scheme to detect, remove and eliminate

The ASP and CASP-GO Auto-DTM Generation System
The Ames Stereo Pipeline (hereafter referred to as ASP) software was developed by the Intelligent Robotics Group at the NASA Ames Research Center, Mountain View, CA [43]. ASP is a suite of open-source automated geodesy and stereo-photogrammetry tools designed for processing planetary imagery captured either from planetary (robotic) orbiters or landers or from Earth observing satellites. The CASP-GO pipeline is based on the ASP framework with specific enhancements, dealing with matching artefacts, gaps, and co-alignment issues, from in-house software previously developed at the Imaging Group at UCL. CASP-GO uses the ASP I/O interface and takes USGS Integrated Software for Imagers and Spectrometers (USGS-ISIS) formatted stereo images (see Figure 2; see explanation of the USGS-ISIS Cube format in https://isis.astrogeology.usgs.gov/documents/ LogicalCubeFormatGuide/LogicalCubeFormatGuide.html (accessed on 21 March 2021)) and a reference base map as inputs. Converting from raw HRSC PDS data into the USGS-ISIS formatted Cube includes the process of format conversion, image pre-denoising, SPICE (see https://naif.jpl.nasa.gov/naif/spiceconcept.html (accessed on 21 March 2021)) kernel initialisation, and map projection processes (see https://isis.astrogeology.usgs.gov/ Application/index.html#Mars_Express (accessed on 21 March 2021)). The core matching algorithm uses a combination of cross-correlation and a 5 th generation of adaptive least squares correlation (ALSC) and region growing matcher called Gotcha [45], which provides accurate and robust sub-pixel conjugate points.
The complete CASP-GO workflow (see Figure 2) for DTM production includes the following 10 steps which are shown in the aforementioned figure and described as follows: (a) ASP image pre-processing including image enhancement (e.g., Laplacian of Gaussian filtering and image normalisation) and stereo rectification; (b) ASP disparity map initialisation using coarse cross-correlation; (c) UCL fast Maximum Likelihood (f-ML) matching and the construction of a "float" initial disparity map; (d) ASP Bayes Expectation Maximisation (EM) weighted affine adaptive sub-pixel cross-correlation; (e) UCL re-defined outlier rejection and gap erosion scheme to detect, remove and eliminate mis-matched and unreliable disparity values; (f) UCL ALSC sub-pixel refinement for previous matches; (g) UCL Gotcha based densification method for filling-in disparity gaps; (h) ASP camera triangulation and DTM creation; (i) UCL co-kriging grid-point interpolation and calculation of height un-Remote Sens. 2021, 13, 1385 6 of 26 certainties for interpolations; (j) UCL ORI co-registration/geocoding with reference to the input base map and DTM adjustment. A detailed description of each of these processing steps can be found in [39]. This processing chain operates without any user interaction and was shown in [39] to be sufficiently robust that it worked without any manual intervention for over 5000 CTX stereo-pairs. Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 28 mis-matched and unreliable disparity values; (f) UCL ALSC sub-pixel refinement for previous matches; (g) UCL Gotcha based densification method for filling-in disparity gaps; (h) ASP camera triangulation and DTM creation; (i) UCL co-kriging grid-point interpolation and calculation of height uncertainties for interpolations; (j) UCL ORI co-registration/geocoding with reference to the input base map and DTM adjustment. A detailed description of each of these processing steps can be found in [39]. This processing chain operates without any user interaction and was shown in [39] to be sufficiently robust that it worked without any manual intervention for over 5,000 CTX stereo-pairs.

Figure 2.
Flow diagram of the UCL CASP-GO auto-DTM processing system. Grey indicates the ASP system, green for UCL processing, and pink for USGS-ISIS processing.

The Hybrid Stereo Processing Chain
As previously mentioned in [39], different DTMs from different stereo pipelines always display different kinds of matching artefacts or gaps. Although these can be partially reduced with down-sampling and interpolation, minimising matching artefacts and better matching coverage are the keys to obtain a high quality DTM with the best possible resolution (finest details).
The quality of the final disparity map used for camera triangulation dominates the quality of the final DTM and is mostly determined by the matching blocks in a DTM pipeline. As the stereo matching algorithms have become more mature over the past decade, practical matching issues mostly come from one or more natural issues of the stereo inputs, for example, robustness to illumination differences, noise level differences, changes or degradations in surface textures, a good stereo convergence angle, and the differences in their native resolutions (although raw images are generally resampled to a fixed reso- Figure 2. Flow diagram of the UCL CASP-GO auto-DTM processing system. Grey indicates the ASP system, green for UCL processing, and pink for USGS-ISIS processing.

The Hybrid Stereo Processing Chain
As previously mentioned in [39], different DTMs from different stereo pipelines always display different kinds of matching artefacts or gaps. Although these can be partially reduced with down-sampling and interpolation, minimising matching artefacts and better matching coverage are the keys to obtain a high quality DTM with the best possible resolution (finest details).
The quality of the final disparity map used for camera triangulation dominates the quality of the final DTM and is mostly determined by the matching blocks in a DTM pipeline. As the stereo matching algorithms have become more mature over the past decade, practical matching issues mostly come from one or more natural issues of the stereo inputs, for example, robustness to illumination differences, noise level differences, changes or degradations in surface textures, a good stereo convergence angle, and the differences in their native resolutions (although raw images are generally resampled to a fixed resolution in planetary imaging, the individual image clarity still varies within the same instrument over the same region). How well a stereo matching algorithm can handle these different types of input issues ultimately determines the quality of the resulting DTM. Although some matching algorithms are more robust to bad quality inputs or differences Remote Sens. 2021, 13, 1385 7 of 26 in image conditions, there is not a single algorithm that can handle all different types of input issues.
In particular, the ASP Bayes EM optimised cross-correlation matching method is a collection of algorithms that are based on image cross-correlation but has been augmented and tuned in efforts to produce more accurate results and obtain robustness over many different input issues. These optimisations, on top of the general "cross-correlation" method, include a coarse-to-fine (pyramidal) approach, filtering, segmentation (tiling), searchingrange optimisation that are based on pre-processing feature matching, and finally the Bayes EM subpixel refinement to adapt to affine and noise variations. The ASP Bayes EM optimised cross-correlation matching method has shown robustness to image noise and good performance on semi-low to highly textured areas, but it is more likely to be affected by the differences in local lighting conditions and contrast, and has poor performance on very-low to non-textured regions (e.g., strongly shaded regions). Also, the ASP Bayes EM optimised cross-correlation matching method is generally robust to many different planetary observation instruments but has some inherent issues from the local (window-based) matching methods, such as block (staircase) artefacts.
On the other hand, the ASP More Global Matching (MGM) method is based on the semi-global matching (SGM) algorithm [48] and the original MGM algorithm that is described in [49]. MGM has been optimised during its implementation in ASP, in terms of 2 D disparity search (instead of 1 D search), multi-resolution hierarchical search, and smoothing of high frequency artefacts for textureless regions. The ASP MGM based matching method has shown less robustness to different instruments (compared to Bayes EM optimised cross-correlation) but has very high matching accuracy and good performance on textureless regions as well as steep slopes. Also, the ASP MGM based matching method has shown higher robustness to handling differences in lighting conditions and contrasts, however, it may produce seam artefacts and diagonal artefacts for large (continuous) low textured regions.
Finally, the UCL Gotcha matcher [45], which was introduced in [39] as a complementary matching algorithm to the ASP Bayes EM optimised cross-correlation matching method, starts from a group of MSA-SIFT "corner" features and gradually grows to their neighbouring pixels with local constraints. The Gotcha matcher is not suitable for global matching due to its high computational cost but is able to produce high accuracy matching results for small images/regions. The Gotcha matcher often produces over-smoothed results on low textured regions, but on the other hand, has shown robustness to local illumination differences, perspective degradations, and a sub-optimal stereo convergence angle. Also, there is no known systematic artefact from Gotcha [45] but it requires a clean (artefact-free) initial matching result as input (from tie-points or other stereo matching disparity maps).
In this work, we experimented with a hybrid method based on CASP-GO and the ASP's new MGM matcher in order to try to exploit the strengths of different matching algorithms whilst minimising the artefacts to produce the best possible quality DTM. The three matching algorithms in the proposed hybrid method contain two local methods and one global method. They provide solutions for different conditions and difficult regions of HRSC stereo images. Although it may not always be feasible to have three parallel algorithms in a general 3 D mapping task, the proposed hybrid method has demonstrated better overall performance compared to the three methods alone.
The proposed HRSC hybrid DTM processing chain follows the standard USGS-ISIS/ASP data I/O methods and takes ISIS formatted HRSC level 2 stereo images as inputs. The same ASP pre-processing methods that were described in [43] and more specifically in [39] are also used.
The complete HRSC hybrid DTM processing chain, is shown in Figure 3, and has the following 6 steps, replacing the original steps (a) to (g) in Figure 2: (a) ASP stereo preprocessing (image normalisation, Laplacian of Gaussian filtering, initial feature matching to decide a searching range); (b) application of the fast-Maximum Likelihood (f-ML) matching to build a "float" initial disparity map (see [39] for more details) and then application of the ASP Bayes EM weighted affine adaptive sub-pixel cross-correlation matching to refine the disparity map; (c) in parallel to (b), Gotcha matcher is applied based on an initial disparity map created from the ASP Normalised Cross-Correlation (NCC) method; (d) in parallel to (b) and (c), the ASP MGM matcher is directly applied to the pre-processed stereo inputs; (e) application of a slightly stricter outlier filtering and border erosion scheme to remove potential bad-matches, as we have three disparity maps from three matching methods in parallel; (f) finally, the three disparity maps are blended together using a "closest-to-median" scheme (described in the next paragraph). Note that the closest-to-median blending is able to discard a bad-match by considering the three matching results and neighbouring disparities but is highly dependent on (e) to discard the majority of bad-matches.
The complete HRSC hybrid DTM processing chain, is shown in Figure 3, and has the following 6 steps, replacing the original steps (a) to (g) in Figure 2: (a) ASP stereo preprocessing (image normalisation, Laplacian of Gaussian filtering, initial feature matching to decide a searching range); (b) application of the fast-Maximum Likelihood (f-ML) matching to build a "float" initial disparity map (see [39] for more details) and then application of the ASP Bayes EM weighted affine adaptive sub-pixel cross-correlation matching to refine the disparity map; (c) in parallel to (b), Gotcha matcher is applied based on an initial disparity map created from the ASP Normalised Cross-Correlation (NCC) method; (d) in parallel to (b) and (c), the ASP MGM matcher is directly applied to the pre-processed stereo inputs; (e) application of a slightly stricter outlier filtering and border erosion scheme to remove potential bad-matches, as we have three disparity maps from three matching methods in parallel; (f) finally, the three disparity maps are blended together using a "closest-to-median" scheme (described in the next paragraph). Note that the closest-to-median blending is able to discard a bad-match by considering the three matching results and neighbouring disparities but is highly dependent on (e) to discard the majority of bad-matches. The closest-to-median blending scheme can be described in the following steps: (a) for any given pixel at (x,y), we have three X-direction disparity values, denoted as dx_1(x,y), dx_2(x,y), dx_3(x,y), and three Y-direction disparity values, denoted as dy_1(x,y), dy_2(x,y), dy_3(x,y), from the three disparity maps; (b) Sort the three x disparity values into each of their 8 adjacent (neighbouring) disparity values. Note that any of the disparity values and their 8 nearest neighbours could be "Nodata" if the matcher fails at that location or at a neighbouring location. "Nodata" is not counted in later steps; (c) Forcibly remove 4 "outliers" that are present at the two ends of the sorted vector from (b) and calculate a median; (d) Any of the three X-direction disparity values, i.e. dx_1(x,y), dx_2(x,y), dx_3(x,y), that are closest to the median from (c) is set as the final blended disparity value. If any of them is "Nodata" or was removed at (c), they will not be counted in this step. If all three of dx_1(x,y), dx_2(x,y), dx_3(x,y) are "Nodata" or were removed at (c), the median value calculated at (c) will be set as the final disparity value for this pixel; (e) Repeat (b)-(d) for Y-direction disparities. The closest-to-median blending scheme can be described in the following steps: (a) for any given pixel at (x,y), we have three X-direction disparity values, denoted as dx_1(x,y), dx_2(x,y), dx_3(x,y), and three Y-direction disparity values, denoted as dy_1(x,y), dy_2(x,y), dy_3(x,y), from the three disparity maps; (b) Sort the three x disparity values into each of their 8 adjacent (neighbouring) disparity values. Note that any of the disparity values and their 8 nearest neighbours could be "Nodata" if the matcher fails at that location or at a neighbouring location. "Nodata" is not counted in later steps; (c) Forcibly remove 4 "outliers" that are present at the two ends of the sorted vector from (b) and calculate a median; (d) Any of the three X-direction disparity values, i.e. dx_1(x,y), dx_2(x,y), dx_3(x,y), that are closest to the median from (c) is set as the final blended disparity value. If any of them is "Nodata" or was removed at (c), they will not be counted in this step. If all three of dx_1(x,y), dx_2(x,y), dx_3(x,y) are "Nodata" or were removed at (c), the median value calculated at (c) will be set as the final disparity value for this pixel; (e) Repeat (b)-(d) for Y-direction disparities.
The advantages of using the hybrid matching system are demonstrated in Figure 4 for 5 different areas that are cropped from one of the HRSC single-strip (h2145_0000) stereo images displaying the X-direction disparity maps (dx) as examples. These examples show different characteristics of the stereo inputs, including (Area-1) hills and dunes with Remote Sens. 2021, 13, 1385 9 of 26 illumination differences and high noise-levels; (Area-2) low-textured flat surfaces, cliffs with good contrast, and small and medium sized craters; (Area-3) heavily shaded regions and steep slopes; (Area-4) sand dunes with contrast differences; (Area-5) bright cliffs and flat rocks with little contrast differences. The results are shown of the 5 cropped areas from the input HRSC level 2 stereo image, the X-disparity maps from the 3 different (aforedescribed) matching algorithms alone, and the final merged X-disparity map using the proposed hybrid system.
The advantages of using the hybrid matching system are demonstrated in Figure 4 for 5 different areas that are cropped from one of the HRSC single-strip (h2145_0000) stereo images displaying the X-direction disparity maps (dx) as examples. These examples show different characteristics of the stereo inputs, including (Area-1) hills and dunes with illumination differences and high noise-levels; (Area-2) low-textured flat surfaces, cliffs with good contrast, and small and medium sized craters; (Area-3) heavily shaded regions and steep slopes; (Area-4) sand dunes with contrast differences; (Area-5) bright cliffs and flat rocks with little contrast differences. The results are shown of the 5 cropped areas from the input HRSC level 2 stereo image, the X-disparity maps from the 3 different (aforedescribed) matching algorithms alone, and the final merged X-disparity map using the proposed hybrid system.   Figure 3) standalone (Disparity-1, Disparity-2, Disparity-3) and merged using the Hybrid DTM processing scheme (Disparity-Final), for 5 different cropped areas, demonstrated using h2145_0000 (location can be found from Figure 8).
We can observe (images at original resolution are provided in supplementary data) that, in general, the f-ML + Bayes EM cross-correlation method (Disparity-1) is able to capture finer texture details but has more gaps (failed matching), the NCC + Gotcha method (Disparity-2) shows higher completeness (fewer gaps) but slightly less detail (smoother), while the MGM method (Disparity-3) shows good overall quality and completeness, but with less detail and some artefacts resulting from inherent (un-smoothed) SGM [48] issues (diagonal lines/8-path streaking artefact [49]). Nevertheless, each of the three methods separately are subject to their own advantages and disadvantages given different types of terrains, lighting and contrast conditions/differences, noise, low-textured and hightextured surfaces. The final merged disparity (Disparity-final), using the afore-described "closest-to-median" scheme shows the highest quality, i.e., best overall completeness and level of detail, and the least artefacts.

DTM Height Correction and Mosaicing
Following the successful creation of 11 single-strip HRSC DTMs, the next stage focuses on spatial and vertical alignment of the new HRSC DTMs to existing DLR HRSC DTMs and ultimately to the only global Mars co-ordinate reference, the MOLA DTM. Within the DLR streamlined system, co-alignment of HRSC DTMs and MOLA is achieved using Bundle Adjustment (BA) and Sequential Photogrammetric Adjustment (SPA) [46]. More specifically, BA uses tie-points and least squares adjustment to model systematic and random errors in a single mathematical model to correct the 3 D geometry, whereas SPA takes a comprehensive set of adjustment procedures and iteratively refines the 3 D geometry. The SPA procedures [46] include (a) minimising 3 D intersection error for individual points; (b) calculation of mean lateral and vertical shifts to MOLA to correct orbit position; (c) calculation of along-track vertical deviations to MOLA to correct pitch angle; and (d) calculation of systematic across-track tilts to MOLA to correct roll angles.
In this work, we propose an alternative method to perform HRSC-to-HRSC and HRSC-to-MOLA co-alignment using a joint 3 D and image co-registration procedure. The joint 3 D and image co-registration process uses existing image and 3 D data as references, and corrects a target image and 3 D data simultaneously, according to image-to-image co-registration for spatial transformation and 3 D-to-3 D co-alignment for vertical transformation. This method uses a combination of image co-registration and 3 D co-alignment techniques and, we believe, is more intuitive and straightforward than traditional photogrammetric modelling approaches. The joint 3 D and image co-registration method also guarantees high co-alignment accuracy both globally and locally. However, it is only applicable when there is partial overlap between the target data and reference data, e.g., adjustment of new HRSC single-strip products to fill in gaps of already corrected HRSC single-strip products. The method can also be used for CTX-to-HRSC, CTX-to-CTX, HiRISE-to-CTX, and HiRISE-to-HiRISE co-alignment and mosaicing. However, for a "blind" HRSC-to-MOLA alignment task (without supporting reference HRSC products), we still suggest following the DLR BA and SPA procedures.
The overall processing chain for the joint 3 D and image co-registration method is illustrated in Figure 5. The joint 3 D and image co-registration method takes our processed HRSC ORIs and DTMs, MOLA DTM, and the existing DLR HRSC level 4 v50+ ORIs and DTMs as inputs, corrects both the UCL HRSC level 4 single-strip ORIs and DTMs as well as the existing DLR HRSC single-strip level-4 DTMs, and creates a DTM mosaic for the whole area. This is achieved through the 5 steps shown in Figure 5 with annotated labels relating to the following: (a) Firstly, a cropped version of the MOLA DTM is upsampled to 50 m/pixel and mosaiced with the existing DLR HRSC level 4 DTMs for ease of computation. In this step, we have prepared a reference DTM mosaic, for follow-up steps, covering the target DTMs, consisting of both high-resolution regions that come from the original DLR HRSC level 4 DTMs and lower resolution (smoother) regions that were filled with up-sampled MOLA DTMs, together covering the target UCL HRSC level 4 DTMs.
(b) Co-registration of the UCL HRSC level 4 ORIs with existing DLR HRSC level 4 ORIs for overlapping regions and recording of the spatial correspondences for all 2 D points (P), here denoted as T P (dX 0 , dY 0 ), for the next step. The ORI co-registrations are performed using Mutual Shape Adapted Scale Invariant Feature Transform (MSA-SIFT), which was introduced in [44] and was also implemented as part of the CASP-GO pipeline [39]. Note that the reason for using image co-registration prior to 3 D co-alignment is that images always have higher spatial resolution with more detailed textures comparing to any derived stereo-derived DTM, and thus result in more accurate spatial adjustments than simply performing a point-cloud alignment using Iterative Closest Point (ICP) based methods. Figure 6 shows that the 3 D co-alignment results for using point-cloud alignment only and using the proposed joint 3 D and image co-registration method. We can see that in both cases, the two HRSC DTMs are "co-aligned" at the DTM scale and could be "seamlessly" blended together. However, when looking at the ORIs, which are at a finer scale of resolution, the differences become much more obvious. (a) Firstly, a cropped version of the MOLA DTM is upsampled to 50 m/pixel and mosaiced with the existing DLR HRSC level 4 DTMs for ease of computation. In this step, we have prepared a reference DTM mosaic, for follow-up steps, covering the target DTMs, consisting of both high-resolution regions that come from the original DLR HRSC level 4 DTMs and lower resolution (smoother) regions that were filled with up-sampled MOLA DTMs, together covering the target UCL HRSC level 4 DTMs.
(b) Co-registration of the UCL HRSC level 4 ORIs with existing DLR HRSC level 4 ORIs for overlapping regions and recording of the spatial correspondences for all 2 D points (P), here denoted as TP (dX0, dY0), for the next step. The ORI co-registrations are performed using Mutual Shape Adapted Scale Invariant Feature Transform (MSA-SIFT), which was introduced in [44] and was also implemented as part of the CASP-GO pipeline [39]. Note that the reason for using image co-registration prior to 3 D co-alignment is that images always have higher spatial resolution with more detailed textures comparing to any derived stereo-derived DTM, and thus result in more accurate spatial adjustments than simply performing a point-cloud alignment using Iterative Closest Point (ICP) based methods. Figure 6 shows that the 3 D co-alignment results for using point-cloud alignment only and using the proposed joint 3 D and image co-registration method. We can see that in both cases, the two HRSC DTMs are "co-aligned" at the DTM scale and could be "seamlessly" blended together. However, when looking at the ORIs, which are at a finer scale of resolution, the differences become much more obvious. (c) Take the spatial correspondences from (b), assuming an initial vertical transformation equals to the absolute difference of height values for the corresponding 2 D points, P, here denoted as dZ 0 , as the initial transformation, i.e., T P (dX 0 , dY 0 , dZ 0 ), adding equally distributed initial transformations as T O (0, 0, 0) for 3 D points (O), and iteratively perform a Non-rigid Iterative Closest Point (NICP) fitting process. The final correspondences between the target UCL HRSC level 4 DTMs (before co-alignment) and the reference DTM mosaic created at step (a) is recorded as T Q (dX, dY, dZ) for all 3 D points (Q; Q=P+O). Note that the final set of 3 D correspondences T Q (dX, dY, dZ) is denser than the initial set of transformations T P (dX 0 , dY 0 , dZ 0 ), as it has added in equally distributed initial "correspondences", i.e., T O , between HRSC and MOLA DTM, where there is no overlapping HRSC ORIs (see Figure 7 for an intuitive illustration), but it is still only a subset (we use between 1/50 and 1/100) of the total HRSC DTM 3 D points due to computational limitations.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 28 (c) Take the spatial correspondences from (b), assuming an initial vertical transformation equals to the absolute difference of height values for the corresponding 2 D points, P, here denoted as dZ0, as the initial transformation, i.e., TP (dX0, dY0, dZ0), adding equally distributed initial transformations as TO (0, 0, 0) for 3 D points (O), and iteratively perform a Non-rigid Iterative Closest Point (NICP) fitting process. The final correspondences between the target UCL HRSC level 4 DTMs (before co-alignment) and the reference DTM mosaic created at step (a) is recorded as TQ (dX, dY, dZ) for all 3 D points (Q; Q=P+O). Note that the final set of 3 D correspondences TQ (dX, dY, dZ) is denser than the initial set of transformations TP (dX0, dY0, dZ0), as it has added in equally distributed initial "correspondences", i.e., TO, between HRSC and MOLA DTM, where there is no overlapping HRSC ORIs (see Figure 7 for an intuitive illustration), but it is still only a subset (we use between 1/50 and 1/100) of the total HRSC DTM 3 D points due to computational limitations. (d) Non-rigid point cloud transformation is performed for the target UCL HRSC level 4 DTMs according to TQ (dX, dY, dZ) and a corresponding spatial transformation is performed for the target UCL HRSC level 4 ORIs according to TQ (dX, dY). Note that in most of the cases for this work, the target HRSC DTM/ORI have overlapping reference HRSC DTMs/ORIs orbital strips on either side. However, where more than one new orbital strip is needed to span a gap, we co-align these one at a time, first selecting at each step, the strip that has the greatest overlap with the data already in place. Then the corrected HRSC DTM/ORI is used as the reference for subsequent co-alignment of the target HRSC DTM/ORI that previously has smaller or no overlapping area. See Figure 7 for a schematic representation, i.e., in case 2, we start from Target-1 (T-1) and Target-3 (T3) HRSC DTM/ORI, then correct T-2 using the corrected T-1, T-3 and MOLA DTM as references.
(e) Finally, the differences between MOLA DTM and each corrected UCL HRSC level 4 DTMs as well as the existing DLR HRSC level 4 DTMs are checked for any remaining systematic issues, e.g., jitter issues or "warped edges". A B-Spline fitting algorithm is used to correct any remaining systematic issues. The B-Spline fitting algorithm uses the computed planar B-Spline surface to represent the HRSC and MOLA DTMs, and then uses the same NICP method that is used in step (c) and (d) to fit the two B-Spline surfaces, and then applies the calculated non-rigid transformation to do a final correction of all HRSC DTMs (both DLR & UCL). More specifically, the NICP algorithm used in step (c) and (d) is based on the Optimal Step NICP (OS-NICP) method [50] which is briefly described as follows. Firstly, assign one initial affine 3 × 4 transformation matrix and a weight term (w) to each 3 D correspondence in Q, to model the distance between the target and reference 3 D point clouds. Note that the weights are set to 0 for the target 3 D points where no corresponding reference 3 D points can be found. The weighted distance term and two additional regularisation terms together form the cost function of OS-NICP (sum of the three terms). The two (d) Non-rigid point cloud transformation is performed for the target UCL HRSC level 4 DTMs according to T Q (dX, dY, dZ) and a corresponding spatial transformation is performed for the target UCL HRSC level 4 ORIs according to T Q (dX, dY). Note that in most of the cases for this work, the target HRSC DTM/ORI have overlapping reference HRSC DTMs/ORIs orbital strips on either side. However, where more than one new orbital strip is needed to span a gap, we co-align these one at a time, first selecting at each step, the strip that has the greatest overlap with the data already in place. Then the corrected HRSC DTM/ORI is used as the reference for subsequent co-alignment of the target HRSC DTM/ORI that previously has smaller or no overlapping area. See Figure 7 for a schematic representation, i.e., in case 2, we start from Target-1 (T-1) and Target-3 (T3) HRSC DTM/ORI, then correct T-2 using the corrected T-1, T-3 and MOLA DTM as references.
(e) Finally, the differences between MOLA DTM and each corrected UCL HRSC level 4 DTMs as well as the existing DLR HRSC level 4 DTMs are checked for any remaining systematic issues, e.g., jitter issues or "warped edges". A B-Spline fitting algorithm is used to correct any remaining systematic issues. The B-Spline fitting algorithm uses the computed planar B-Spline surface to represent the HRSC and MOLA DTMs, and then uses the same NICP method that is used in step (c) and (d) to fit the two B-Spline surfaces, and then applies the calculated non-rigid transformation to do a final correction of all HRSC DTMs (both DLR & UCL).
More specifically, the NICP algorithm used in step (c) and (d) is based on the Optimal Step NICP (OS-NICP) method [50] which is briefly described as follows. Firstly, assign one initial affine 3 × 4 transformation matrix and a weight term (w) to each 3 D correspondence in Q, to model the distance between the target and reference 3 D point clouds. Note that the weights are set to 0 for the target 3 D points where no corresponding reference 3 D points can be found. The weighted distance term and two additional regularisation terms together form the cost function of OS-NICP (sum of the three terms). The two regularisation terms are the stiffness term to penalise transformations of neighbouring 3 D points, and the landmark term to initialise and penalise the transformations of landmarks. The landmarks, in our case, are the 3 D correspondences calculated via HRSC ORI-to-ORI co-registration, i.e., T P (dX 0 , dY 0 , dZ 0 ). OS-NICP iteratively finds the optimal 3 D correspondences T Q , while lowering the stiffness weights (to allow more localised transformations), until the cost function is minimised. Note that the correct alignment can also be found without landmarks in OS-NICP with a wider range of initial conditions.
In the last step (e), a B-Spline fitting algorithm is used to remove any remaining systematic errors for both the UCL and DLR HRSC level 4 DTMs. The motivation for this step is based on our observation that some systematic errors still exist for both the DLR HRSC level 4 DTMs (after BA/SPA correction) and the UCL HRSC level 4 DTMs (after the joint 3 D and image co-registration process). Considering two DTMs with the same resolution, size, and zero noise, if they are spatially co-registered, then the absolute height difference between the two DTMs would be the theoretical transformation to co-align each point of the two DTMs. However, we cannot use an absolute difference mapping to perform point-to-point corrections for HRSC DTM against MOLA DTM in our case, since that will remove all the details in the higher resolution HRSC DTM and make it identical to the MOLA DTM. Therefore, in this work, we propose using a B-Spline fitting algorithm to correct remaining systematic errors in both the DLR and UCL HRSC level 4 DTMs.
The B-Spline fitting algorithm firstly computes two planar B-Spline surfaces, to represent the large-scale topography of HRSC and MOLA DTMs at a density of 2 km/pixel (~5 times the MOLA scale), then fit the two planar curves with NICP, and finally apply the calculated non-rigid transformation to the HRSC DTMs. In terms of computing the planar B-Spline curve representations, we used the asymmetric distance minimisation method described in [51] (implementation available through the Point Cloud Library's Surface module at: http://pointclouds.org/documentation/tutorials/#surface (accessed on 21 March 2021)).
A high-level workflow for fitting a B-Spline surface to the point-clouds has the following 4 steps: (1) initialise the B-Spline surface by calculating a bounding circle of the target point cloud and set the initial control points; (2) increase the number of control-points and subsequently fit the refined B-spline surface; (3) project the target point-cloud into the parametric domain using the closest points to the B-spline surface and cut off the overlapping regions of the surface; (4) repeat (2) and (3) until the desired density of control points is reached.
Finally, the updated UCL and DLR HRSC level 4 DTMs are used to create a mosaic for the whole Valles Marineris area. We use the ASP "priority-blending" method to create the mosaic, in which case, DTMs are blended in the order of highest-to-lowest resolutions.

Automated Co-Registration and Orthorectification
Initially, the DLR HRSC level 4 ORIs (ND4 products) and the UCL HRSC level 4 ORIs, which were orthorectified and co-registered using HRSC level 2 nadir panchromatic images (ND2 products), were used to create the ORI mosaic of Valles Marineris. However, we found some of the UCL level 4 ORIs, that were produced from the ND2 products, have missing data in different areas of the strip (the missing data comes from the original raw HRSC data records). Also, some of the DLR level 4 ORIs (ND4 products) are narrower than the level 2 stereo images (S12 and S22 products), leaving a seam gap between two adjacent ORIs (see Figure 8). Therefore, additional HRSC level 2 stereo images (S12 and S22 products) are used to cover these gap regions, and further orthorectification and co-registrations are performed. Also, due to missing data, gaps appear at different locations in the different HRSC level 2 images (the level 2 nadir panchromatic and stereo images, i.e., ND2, S12, and S22 products) for a same strip. Multiple HRSC level 2 images are therefore sometimes required in order to create the final Valles Marineris ORI mosaic without any gaps. The additional HRSC level 2 images (ND2, S12, and S22) that were needed to produce overlapping new level 4 ORIs in order to fill in the gaps of the original 11 UCL level 4 ORIs are h2178_0000, h2145_0000, h2134_0000, and h2112_0000. The additional HRSC level 2 images (S12 and S22) that were needed to produce new "wider" level 4 ORIs in order to fill the thin seam gaps of the original DLR level 4 ORIs are h0931_0000 and h0920_0000. The issues are demonstrated in Figure 8 for the original HRSC level 4 ORIs from UCL and DLR using the nadir panchromatic images only (ND4 and ND2 images) showing missing data and gaps, in comparison to adding in the additional overlapped UCL HRSC level 4 ORIs (coregistered and orthorectified using ND2, S12, and S22 images).
In this work, we use an automated image co-registration and orthorectification workflow to produce multiple overlapping HRSC level 4 ORIs from both level 2 nadir panchromatic and stereo images (ND2, S12, S22 products). The implementation uses our inhouse MSA-SIFT tie-point based image co-registration pipeline [39,44] and the ASP's image orthorectification workflow (the "mapproject" function) [43]. The processing follows the showing completeness before (a,c) and after (b,d) adding more additional overlapping UCL level 4 ORIs, produced from both the nadir and stereo level 2 images (ND2 and S12/S22 products), to fill in the missing data and seam gaps. Therefore, additional HRSC level 2 stereo images (S12 and S22 products) are used to cover these gap regions, and further orthorectification and co-registrations are performed. Also, due to missing data, gaps appear at different locations in the different HRSC level 2 images (the level 2 nadir panchromatic and stereo images, i.e., ND2, S12, and S22 products) for a same strip. Multiple HRSC level 2 images are therefore sometimes required in order to create the final Valles Marineris ORI mosaic without any gaps. The additional HRSC level 2 images (ND2, S12, and S22) that were needed to produce overlapping new level 4 ORIs in order to fill in the gaps of the original 11 UCL level 4 ORIs are h2178_0000, h2145_0000, h2134_0000, and h2112_0000. The additional HRSC level 2 images (S12 and S22) that were needed to produce new "wider" level 4 ORIs in order to fill the thin seam gaps of the original DLR level 4 ORIs are h0931_0000 and h0920_0000. The issues are demonstrated in Figure 8 for the original HRSC level 4 ORIs from UCL and DLR using the nadir panchromatic images only (ND4 and ND2 images) showing missing data and gaps, in comparison to adding in the additional overlapped UCL HRSC level 4 ORIs (co-registered and orthorectified using ND2, S12, and S22 images).
In this work, we use an automated image co-registration and orthorectification workflow to produce multiple overlapping HRSC level 4 ORIs from both level 2 nadir panchromatic and stereo images (ND2, S12, S22 products). The implementation uses our inhouse MSA-SIFT tie-point based image co-registration pipeline [39,44] and the ASP's image orthorectification workflow (the "mapproject" function) [43]. The processing follows the steps of map projection with MOLA, image co-registration with already processed and co-registered HRSC level 4 ORIs, and orthorectification with the HRSC DTM Valles Marineris mosaic created at the previous stage. A flow diagram of the automated image co-registration and orthorectification processing chain is shown in Figure 9.

Brightness Correction and Mosaicing of ORIs
A particular difficulty in compiling image mosaics for the surface of Mars arises from the varying optical properties of the atmosphere. The quantity of dust carried by the atmosphere changes as a consequence of seasonal weather phenomena, affecting camera observations to a greater or smaller degree over different parts of the planet. This is seen in individual images as variations in clarity and contrast, resulting from how much light is scattered both onto the ground -from directions other than the source of illumination -as well as into and out of the camera's line of sight. The image strips shown in Figure 8 illustrates the problem. It is difficult to make a physical correction for the effect, because the images are the only source of information on the instantaneous atmospheric dust content. A model-based measurement of dust content from an image requires some prior knowledge of how the surface should appear: once more, the image is the only source for the information, making the problem underdetermined.
We are able to approximate a solution by making use of the fact that large scale albedo variations on Mars are already known. The MGS's TES was used to construct a global albedo map at about 7 km/pixel [39]. By introducing the constraint that an HRSC mosaic should conform to this albedo at coarse resolution, we are able to bring all the images into mutual brightness consistency. The main steps of the procedure are as follows [41]: 1) Lambert correction. The main effect of this step is to compensate the variation of illumination across an image caused by the curvature of the planet.

Brightness Correction and Mosaicing of ORIs
A particular difficulty in compiling image mosaics for the surface of Mars arises from the varying optical properties of the atmosphere. The quantity of dust carried by the atmosphere changes as a consequence of seasonal weather phenomena, affecting camera observations to a greater or smaller degree over different parts of the planet. This is seen in individual images as variations in clarity and contrast, resulting from how much light is scattered both onto the ground-from directions other than the source of illumination-as well as into and out of the camera's line of sight. The image strips shown in Figure 8 illustrates the problem. It is difficult to make a physical correction for the effect, because the images are the only source of information on the instantaneous atmospheric dust content. A model-based measurement of dust content from an image requires some prior knowledge of how the surface should appear: once more, the image is the only source for the information, making the problem underdetermined.
We are able to approximate a solution by making use of the fact that large scale albedo variations on Mars are already known. The MGS's TES was used to construct a global albedo map at about 7 km/pixel [39]. By introducing the constraint that an HRSC mosaic should conform to this albedo at coarse resolution, we are able to bring all the images into mutual brightness consistency. The main steps of the procedure are as follows [41]: (1) Lambert correction. The main effect of this step is to compensate the variation of illumination across an image caused by the curvature of the planet.   (2), but this time using a smaller cell size (9 × n).
Overlapping images are feathered together over a 40-pixel range. (4) Image sequence. The overlapping sequence for the mosaic is optimised through visual inspection. The starting sequence places the shortest ground sampling distance at the top of the mosaic. Lower quality images are moved down in the sequence by estimating a longer effective ground sampling distance by comparison with neighbouring images. (5) Contrast adjustment. Images with higher atmospheric scattering appear in the mosaic as areas of low contrast. The contrast is increased either uniformly or with a multipoint interpolated along-track change. (6) Iteration. The mosaic is regenerated and steps (4) and (5) reconsidered until the result is as close to visually homogeneous as possible.

HRSC DTM Mosaic for Valles Marineris
In this work, we created 11 high-quality HRSC level 4 DTMs, using a hybrid DTM processing chain that combine matching results from three different matchers, using the HRSC level 2 stereo images, to cover the gap regions from existing DLR HRSC level 4 DTMs for the Valles Marineris area. The resulting UCL level 4 DTMs are then simultaneously coaligned with the DLR products and MOLA DTM, via a joint 3 D and image co-registration process, to allow seamless mosaicing for the whole of Valles Marineris. The mosaic which just includes the DLR DTM products is shown in Figure 10. The final resulting HRSC DTM mosaic (50 m/pixel), using the 71 existing DLR DTMs and 11 newly created DTMs from UCL, is shown in Figure 11. We note that the new HRSC DTM mosaic is much more complete covering the whole of Valles Marineris area. Some cropped examples are provided in Figure 12 for zoom-in views over different terrain types, including hill/rock chaos, craters, flat surface, and valley/cliffs, within the Valles Marineris area.
Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 28 5) Contrast adjustment. Images with higher atmospheric scattering appear in the mosaic as areas of low contrast. The contrast is increased either uniformly or with a multipoint interpolated along-track change. 6) Iteration. The mosaic is regenerated and steps (4) and (5) reconsidered until the result is as close to visually homogeneous as possible.

HRSC DTM Mosaic for Valles Marineris
In this work, we created 11 high-quality HRSC level 4 DTMs, using a hybrid DTM processing chain that combine matching results from three different matchers, using the HRSC level 2 stereo images, to cover the gap regions from existing DLR HRSC level 4 DTMs for the Valles Marineris area. The resulting UCL level 4 DTMs are then simultaneously co-aligned with the DLR products and MOLA DTM, via a joint 3 D and image coregistration process, to allow seamless mosaicing for the whole of Valles Marineris. The mosaic which just includes the DLR DTM products is shown in Figure 10. The final resulting HRSC DTM mosaic (50 m/pixel), using the 71 existing DLR DTMs and 11 newly created DTMs from UCL, is shown in Figure 11. We note that the new HRSC DTM mosaic is much more complete covering the whole of Valles Marineris area. Some cropped examples are provided in Figure 12 for zoom-in views over different terrain types, including hill/rock chaos, craters, flat surface, and valley/cliffs, within the Valles Marineris area.
Also, the final resulting HRSC DTM mosaic has shown improved co-alignment accuracy against the MOLA DTM, which can be observed in Figure 13 and Figure 14, in comparison to the HRSC-to-MOLA difference maps (at the scale of 200 m/pixel) created before and after the proposed DTM co-alignment method (before mosaicing).      Also, the final resulting HRSC DTM mosaic has shown improved co-alignment accuracy against the MOLA DTM, which can be observed in Figures 13 and 14, in comparison to the HRSC-to-MOLA difference maps (at the scale of 200 m/pixel) created before and after the proposed DTM co-alignment method (before mosaicing).
alignments have been significantly reduced. Even though some local disagreement still exists in a few areas, most of the area, i.e., ~80.4% of all DTMs, has less than 30 m difference comparing to MOLA and about half, i.e., ~55.7% of all DTMs, has less than 15 m difference comparing to MOLA. On the other hand, without 3 D co-alignment, only ~60.3% of the HRSC DTM pixels has less than 30 m difference to MOLA, and ~49.1% of the HRSC DTM pixels has less than 15 m difference to MOLA. This suggests the vertical accuracy, after the joint 3 D and image co-registration, has had a significant positive improvement for large misalignments (>30 m), and some fair improvement for small misalignments (<30 m). Figure 13. HRSC-to-MOLA difference map before 3 D co-alignment. Figure 13. HRSC-to-MOLA difference map before 3 D co-alignment. As for the mosaicing process, we used the higher resolution HRSC level 4 DTMs (57 from DLR and 11 from UCL at 50 m/pixel) in front and the lower resolution ones on the back (6 from DLR at 100 m/pixel, then 8 from DLR at 150 m/pixel). This means that although the Valles Marineris DTM mosaic is gridded at 50 m/pixel, some regions are mosaiced with upsampled 100 m/pixel or 150 m/pixel DTMs. Therefore, for those areas, they will show less detail than the areas that were mosaiced with native 50 m/pixel DTMs. Here we provide a colourised footprint figure, see Figure 15, to distinguish the native mapping resolutions of the Valles Marineris DTM mosaic. The lower DTM resolution is due to the As shown in the difference maps comparing the single-strip HRSC level 4 DTMs, before (refer to Figure 13) and after (refer to Figure 14) our proposed joint 3 D and image co-registration, against MOLA DTM, not only removes the systematic errors of the singlestrip HRSC level 4 DTMs, e.g. jitter [38] and "tilted ends" (refer to Figure 13 for an example of these issues with jitter shown in the blue-to-red-to-blue differences and "tilted ends" where the end of a strip shows blue differences), have been corrected, but also local misalignments have been significantly reduced. Even though some local disagreement still exists in a few areas, most of the area, i.e.,~80.4% of all DTMs, has less than 30 m difference comparing to MOLA and about half, i.e.,~55.7% of all DTMs, has less than 15 m difference comparing to MOLA. On the other hand, without 3 D co-alignment, only~60.3% of the HRSC DTM pixels has less than 30 m difference to MOLA, and~49.1% of the HRSC DTM pixels has less than 15 m difference to MOLA. This suggests the vertical accuracy, after the joint 3 D and image co-registration, has had a significant positive improvement for large misalignments (>30 m), and some fair improvement for small misalignments (<30 m).
As for the mosaicing process, we used the higher resolution HRSC level 4 DTMs (57 from DLR and 11 from UCL at 50 m/pixel) in front and the lower resolution ones on the back (6 from DLR at 100 m/pixel, then 8 from DLR at 150 m/pixel). This means that although the Valles Marineris DTM mosaic is gridded at 50 m/pixel, some regions are mosaiced with upsampled 100 m/pixel or 150 m/pixel DTMs. Therefore, for those areas, they will show less detail than the areas that were mosaiced with native 50 m/pixel DTMs. Here we provide a colourised footprint figure, see Figure 15, to distinguish the native mapping resolutions of the Valles Marineris DTM mosaic. The lower DTM resolution is due to the lower resolution input stereo images, i.e., 12.  As for the mosaicing process, we used the higher resolution HRSC level 4 DTMs (57 from DLR and 11 from UCL at 50 m/pixel) in front and the lower resolution ones on the back (6 from DLR at 100 m/pixel, then 8 from DLR at 150 m/pixel). This means that although the Valles Marineris DTM mosaic is gridded at 50 m/pixel, some regions are mosaiced with upsampled 100 m/pixel or 150 m/pixel DTMs. Therefore, for those areas, they will show less detail than the areas that were mosaiced with native 50 m/pixel DTMs. Here we provide a colourised footprint figure, see Figure 15, to distinguish the native mapping resolutions of the Valles Marineris DTM mosaic. The lower DTM resolution is due to the lower resolution input stereo images, i.e., 12.5 m/pixel images produce 50 m/pixel DTMs, 25 m/pixel images produce 100 m/pixel DTMs, 50 m/pixel images produce 150 m/pixel DTMs, so Figure 15 also indicates the high-resolution and low-resolution distributions of the Valles Marineris ORI mosaic (see Section 3.2). Examples of HRSC-to-MOLA profile locations are shown in Figure 16 and the results are shown in Figure 17 for Profile A, Figure 18 for Profile B, Figure 19 for Profile C, Figure Examples of HRSC-to-MOLA profile locations are shown in Figure 16 and the results are shown in Figure 17 for Profile A, Figure 18 for Profile B, Figure 19 for Profile C, Figure 20 for Profile D, and Figure 21 for Profile E. Profile A covers a broad extent, crossing multiple DLR HRSC level 4 DTMs and UCL HRSC level 4 DTMs, and has shown smooth transitions between UCL and DLR products for the whole mosaic product. Profile B is a shorter line that contains only one UCL products in between two DLR products, a smooth transition from UCL to DLR products can be clearly observed at both sides. Also, Profile B shows no systematic cross-track offset at both ends of the UCL product. Profile C is a long vertical line crossing only one UCL product and shows no systematic along-track offset over a larger extent. On the other hand, Profile D is a shorter vertical line crossing the edge of the same UCL product shown in Profile C. Profile D shows the details of the HRSC DTM compared to MOLA DTM as well as good alignment between the two. Also, Profile D has shown that, after our 3 D co-alignment process, there is no along-track warping issue near the DTM edge, which has been a common issue for stereo-derived DTMs. Finally, Profile E shows a short horizontal line, crossing both sides of the border of the UCL product. Profile E also shows the details of the good alignment of the HRSC DTM against the MOLA DTM. Also, Profile E shows that there are no across-track warping errors at either side of the DTM edges.
HRSC DTM compared to MOLA DTM as well as good alignment between the two. Also, Profile D has shown that, after our 3 D co-alignment process, there is no along-track warping issue near the DTM edge, which has been a common issue for stereo-derived DTMs. Finally, Profile E shows a short horizontal line, crossing both sides of the border of the UCL product. Profile E also shows the details of the good alignment of the HRSC DTM against the MOLA DTM. Also, Profile E shows that there are no across-track warping errors at either side of the DTM edges.   long vertical line crossing only one UCL product and shows no systematic along-track offset over a larger extent. On the other hand, Profile D is a shorter vertical line crossing the edge of the same UCL product shown in Profile C. Profile D shows the details of the HRSC DTM compared to MOLA DTM as well as good alignment between the two. Also, Profile D has shown that, after our 3 D co-alignment process, there is no along-track warping issue near the DTM edge, which has been a common issue for stereo-derived DTMs. Finally, Profile E shows a short horizontal line, crossing both sides of the border of the UCL product. Profile E also shows the details of the good alignment of the HRSC DTM against the MOLA DTM. Also, Profile E shows that there are no across-track warping errors at either side of the DTM edges.

HRSC ORI Mosaic for Valles Marineris
The mosaic from the orthoimages after Lambertian correction, brightness referencing, manual contrast adjustment and image sequencing is shown in Figure 22. In general, the transitions from one image to another are seamless, although in a few cases it is possible to see some change in image quality. It should be noted that HRSC observations, due to its elliptical orbit has equator crossing overpasses at different times of day, so that the solar illumination is not usually consistent between one strip and the next. There are also

HRSC ORI Mosaic for Valles Marineris
The mosaic from the orthoimages after Lambertian correction, brightness referencing, manual contrast adjustment and image sequencing is shown in Figure 22. In general, the transitions from one image to another are seamless, although in a few cases it is possible to see some change in image quality. It should be noted that HRSC observations, due to its elliptical orbit has equator crossing overpasses at different times of day, so that the solar illumination is not usually consistent between one strip and the next. There are also

HRSC ORI Mosaic for Valles Marineris
The mosaic from the orthoimages after Lambertian correction, brightness referencing, manual contrast adjustment and image sequencing is shown in Figure 22. In general, the transitions from one image to another are seamless, although in a few cases it is possible to see some change in image quality. It should be noted that HRSC observations, due to its elliptical orbit has equator crossing overpasses at different times of day, so that the solar illumination is not usually consistent between one strip and the next. There are also a couple of images which show morning fog on the floor of the valley for which there is still no alternative coverage: the transition from fog to no fog remains visible. The product was constructed at 12.5 m/pixel in 16-bit greyscale, and totals 120 GB.

Access to the Products
The resultant HRSC DTM and ORI mosaic for the Valles Marineris area will be viewable through an interactive webGIS system (http://www.i-mars.eu/web-gis) developed at the Free University Berlin [52]. The mosaic products are also being made publicly availa-

Access to the Products
The resultant HRSC DTM and ORI mosaic for the Valles Marineris area will be viewable through an interactive webGIS system (http://www.i-mars.eu/web-gis (accessed on 21 March 2021)) developed at the Free University Berlin [52]. The mosaic products are also being made publicly available through the ESA Guest Storage Facility (GSF) via their DOI [53] (see list of DOIs at https://www.cosmos.esa.int/web/psa/ucl-mssl_meta-gsf (accessed on 21 March 2021)). A link will also be established between the webGIS and the source-file and vice versa to ease visualisation and access/downloading.

Discussion
This paper presented our 3 D mapping and mosaicing work of Valles Marineris using HRSC level 2 stereo images and pre-existing level 4 products. We highlighted a collection of methods, in order to produce the resultant HRSC DTM and ORI mosaic, which include a new hybrid HRSC DTM processing chain for single-strip DTM production, a joint 3 D and image co-registration system to co-align all single-strip products with each other and simultaneously co-align with MOLA DTM, and finally an ORI co-registration and mosaicing method to create the ORI mosaic.
The proposed hybrid stereo matching method may not be practical for general stereo mapping tasks due to its high computational load. However, it is capable of achieving better quality for small-to-medium sized 3 D reconstruction jobs comparing to the experimental single-matcher driven methods. The hybrid method can also be used to combine other unique matchers as far as they provide unique advantages over different circumstances of stereo matching. Similarly, if a stereo matching algorithm is highly dependent on the set parameters, then this hybrid method can also be applied to improve the stereo result by merging different disparity maps produced with different parameters of the same matcher. In future experiments, we plan to explore the advantages and disadvantages of more global matchers and, as well as deep-learning based methods, that have been proven successful in the field of stereo vision (https://vision.middlebury.edu/stereo (accessed on 21 March 2021)) and try different combinations for more robust and more accurate matching of planetary images. The disparity blending scheme proposed at the disparity merging step, can also be improved in future studies, as currently it only uses a small number (nearest 8 neighbours) of spatial information of the disparities. For example, if the disparity blending scheme considers a larger group of relevant neighbouring disparity values using the idea of co-kriging [54], the blended disparity might have better robustness towards matching artefacts and gaps.
Note that the HRSC stereo images do not normally have significant brightness/contrast /clarity/shadow differences between different off-nadir views as the HRSC instrument acquires continuous along-track acquisitions in a single orbit for stereo images (minimal time interval of up to~30 seconds). However, for across-track repeat images such as CTX and HiRISE stereo, there can be substantial differences in image brightness/contrast/shadow, and the aforementioned clarity differences caused by different native resolution and noise level is a common challenge for HiRISE stereo. The proposed HRSC hybrid DTM processing chain can be further applied to produce robust results over CTX and HiRISE stereo images.
As for the proposed joint 3 D and image co-registration method, one of the key advantages is its high spatial and vertical alignment accuracy to the reference data, as we have demonstrated in the difference map shown in Figure 14 and follow-up profile analysis. The joint 3 D and image co-registration method makes use of the finer-scale image-to-image correspondences as initial and "landmark" correspondences in a coarserscale 3 D-to-3 D matching scenario. However, the main limitation for this method is that it always requires a reference data. Where there is more than one sequential co-alignment steps needed, as described in Section 2.4 (d), there may be an accumulated HRSC-to-HRSC error in the centre, even though the HRSC-to-MOLA error is minimised. Also, the BA/SPA procedure updates the orientation data to correct image and DTM distortions. Therefore, it can be advantageous to perform such relative adjustment as a first step before any further improvements on 3 D co-alignment, and even essential if there is a large distortion in any of the input HRSC level 2 stereo images. No such incidences were observed here.
In this work, we assume the existing DLR HRSC level 4 v50+ products (initially) and MOLA DTM are trustable reference data. Note that our proposed joint 3 D and image co-registration method is not meant as a replacement of the BA/SPA procedure for initial (and batch) HRSC corrections. Instead, it is meant to act as an alternative method for new target HRSC corrections, where there are sufficient overlapping reference products, and also a complementary method to improve the HRSC-to-MOLA co-alignment after BA/SPA. All HRSC single-strip DTMs (both UCL and DLR) are finally corrected against MOLA DTM, assuming there are no significant artefacts or topographical changes at the MOLA scale.
Regarding the final HRSC DTM mosaic, overall quality has been visually checked and co-alignment accuracy has been quantified. However, the dataset is very large, and the stereo images are limited, so some local artefacts cannot be fully removed. There are known artefacts, including: tiling artefacts (seam lines), small interpolated-regions, smoothed area that were caused by using "foggy" HRSC stereo images (no other options were available as we have checked), occasional artefacts inside small craters, quilting artefacts at heavily shaded cliff areas, and the junction of lower resolution with higher resolution single-strip DTMs. These artefacts are illustrated in Figure 23. Methods have been developed for their removal, where feasible but were beyond the scope of this paper.

Conclusions
In this paper, we have introduced a complete processing system to produce HRSC level 4 DTMs over an extensive area, co-alignment with existing HRSC level 4 and MOLA DTMs, and the creation of a DTM and ORI mosaic that covers the large canyon system of Valles Marineris. We used a hybrid DTM processing chain that combines the matching result from three different matchers to produce high-quality HRSC single-strip DTMs. Then we introduced our joint 3 D and image co-registration processing chain for highaccuracy HRSC-to-HRSC and HRSC-to-MOLA co-alignment and DTM mosaicing. Finally, automated image co-registration and orthorectification, and brightness correction methods are introduced for the production of Valles Marineris ORI mosaic. (c) large "smoothed out" region inside the valley caused by using foggy input images that are not replaceable; (d) some minor artefact in very small craters; (e) quilting artefact at heavily shaded cliff areas; (f) the effect of joining higher resolution (50 m/pixel; from the left side) and lower resolution (150 m/pixel; from the right side) single-strip DTMs.

Conclusions
In this paper, we have introduced a complete processing system to produce HRSC level 4 DTMs over an extensive area, co-alignment with existing HRSC level 4 and MOLA DTMs, and the creation of a DTM and ORI mosaic that covers the large canyon system of Valles Marineris. We used a hybrid DTM processing chain that combines the matching result from three different matchers to produce high-quality HRSC single-strip DTMs.
Then we introduced our joint 3 D and image co-registration processing chain for highaccuracy HRSC-to-HRSC and HRSC-to-MOLA co-alignment and DTM mosaicing. Finally, automated image co-registration and orthorectification, and brightness correction methods are introduced for the production of Valles Marineris ORI mosaic.
As the result, we generated 11 new HRSC level 4 single-strip DTMs and ORIs that have filled the gap regions from existing DLR HRSC level 4 products. We then co-aligned the new HRSC DTM/ORIs with existing DLR HRSC DTM/ORIs, and for both of the UCL and DLR products, co-aligned with MOLA. The Valles Marineris DTM mosaic is therefore created using 82 co-aligned single-strip DTMs and the Valles Marineris ORI mosaic is subsequently created using 94 co-registered, brightness and contrast corrected single-strip ORIs. HRSC-to-MOLA difference map and profile analysis are demonstrated, examples are shown for evaluation of the mosaic product, and known artefacts are summarised and demonstrated.
The Valles Marineris mosaic product will be publicly available through the iMars webGIS and ESA GSF system. In the near future, we will look at applying the shape-fromshading method such as [55] to further improve quality and reduce artefact in the DTMs. In the future, we also plan to also build a HRSC colour ORI mosaic, and as well as the CTX DTM/ORI mosaic by co-aligning the CTX products we have processed for Valles Marineris, some of which have already been released through (https://www.cosmos.esa.int/web/ psa/UCL-MSSL_iMars_CTX_v1.0 (accessed on 21 March 2021)).

Funding:
The research leading to these results is receiving funding from the UKSA Aurora programme (2018-2021) under grant no. ST/S001891/1. Partial funding has been received from the STFC "MSSL Consolidated Grant" ST/K000977/1. SJC is thankful to the CNES for supporting her HRSC-related work.

Data Availability Statement:
The data products are all available via the DOI being added to the summary page for UCL products at https://www.cosmos.esa.int/web/psa/ucl-mssl_meta-gsf (accessed on 21 March 2021).