FOSS4G DATE for DSMs generation from tri-stereo optical satellite images: development and first results

ABSTRACT The topic of this research is the identification of an innovative strategy for DSMs generation from optical satellite tri-stereo imagery, exploiting efficient dense matching algorithms from computer vision, without losing a rigorous photogrammetric approach. The main challenge is related to the epipolarity resampling for satellite images, for which it is not straightforward the epipolar geometry achievement, due to their multiple projection centres. The Ground quasi-Epipolar Images (GrEI) are the core of this original strategy for DSMs generation, able to return reliable and accurate solutions. The FOSS4G DATE has been modified since its original implementation, in order to obtain better results in terms of DSMs generated, and a new version has been released. The workflow provides a ground projection for the raw satellite images, exploiting RPFs model and a coarse DSM. In order to bring the ground images in the quasi-epipolar geometry, a preliminary rotation is estimated for both images to align them in the disparity prevailing direction. Then, a refinement of their relative position is performed. These orthorectified images can act as GrEI, and can undergo a dense image matching procedure. Tests have been carried out with Pléiades triplets over Trento (Northern Italy) and very good results are achieved.


Introduction
In this paper, the new procedure for digital surface models (DSMs) generation implemented in the Free and Open Source Software for Geospatial (FOSS4G) named Digital Automatic Terrain Extractor (DATE), and its application to optical tri-stereo satellite images are presented.
Since its original implementation, DATE has been modified in order to obtain better results in terms of DSMs generated. As a matter of fact, starting from the version implemented during the Google Summer of Code program, a new release has been developed, improving the achievable results. A new and more reliable procedure for epipolarity resampling has been developed and tested, along with an original methodology for pixel to meter disparity conversion. Furthermore, a merging strategy for tri-stereo acquisition exploitation has been designed and implemented: as a matter of fact, DSMs merging could be performed to exploit multi-view data from tri-stereo acquisitions, useful for minimising areas that cannot be reconstructed due to occlusions, e.g. in urban areas where a lower convergence angle is of advantage, or to eliminate outliers and leading to a more complete 3D restitution.
DATE is a totally automatic console application conceived as an Open Source Software Image Map (OSSIM) plug-in; it has been developed starting from summer 2014 in the framework of 2014 Google Summer of Code, within a project entitled "Photogrammetric image processing: DSM generation tool for OSSIM" (Di Rita, 2014). OSSIM, the software in which DATE has been designed and developed, is a FOSS4G (Moreno-Sanchez, 2012), that shows up as a set of libraries and applications for image, maps and terrain data processing, photogrammetric and remote sensing applications; it belongs to the wide OSGeo family, has been actively developed since 1996 and appears with a great amount of features and a very robust architecture (OSSIM wiki, 2015). Amongst other benefits, OSSIM is able to handle a large amount of data as it manages to process different images in a sequential and totally automatic way.
The implemented tool is based on a hybrid procedure, whereby photogrammetric and computer vision algorithms are mixed in order to take the best of both. The synergy between computer vision and photogrammetry has taken hold over a decade ago and it has already returned outstanding results (Pierrot-Deseilligny & Isabelle, 2011). The main goal was to develop a fully automatic, efficient, precise and accurate tool for DSM generation from satellite imagery, paying special attention in finding a solution for the epipolar resampling for images acquired by optical and SAR satellite sensors: this is performed in the object space thanks to the images ground projection (see section "Images ground projection"), and it is based on a rototranslation transformation model for Ground quasi-Epipolar Images (GrEI) generation (see section "GrEI generation") and on the coarse-to-fine pyramidal scheme adopted.
The goal of the paper is therefore to describe DATE new procedure and functioning and to present and discuss some results obtained through tri-stereo Pléiades satellite imagery. In section "DATE workflow", the DATE workflow is described, pointing out its main changes with respect to the previous version (Di Rita, Martina, & Nascetti, 2017); in section "DATE testing", the results related to some tests performed using Pléiades imagery are presented and discussed; finally in section "Conclusions and future prospects", some conclusions are drawn and future prospects are outlined.

DATE workflow
As extensively described in Di Rita et al. (2017), a tool for automatic DSM generation has been developed and various algorithms have been adapted and implemented in order to generate a DSM from an optical satellite stereopair in a totally automatic way. In particular, here an improvement of the developed strategy, in order to obtain a robust solution in terms of quasi-epipolar images achieved from the original stereopairs, is defined and described.
The research founding ideas are detailed in Di Rita et al. (2017), hereafter only two well-known assumptions on which the strategy bases on are recalled. First of all, assuming to have two or more satellite images acquired over the same area, if we know exactly the orientation for each image and we know perfectly the morphology of the area (DSM), all these images orthorectified using this DSM produce equal orthoimages in the same reference frame, apart from local occlusions and radiometric differences. Secondly, it is well known that satellite pushbroom sensors acquire different images on the same orbit (including stereopairs and triplets) with smooth changes of attitude over time; these changes can be modelled with low-degree polynomial functions (usually up to third degree), so that the images can be reprojected in a quasi-epipolar geometry (Toutin, 2004;Poli, 2005;Crespi, Fratarcangeli, Giannone, & Pieralice, 2012). Therefore, the orthoimages generated using a coarse a priori DSM are affected by local pixel disparities lined up in a prevailing direction. The idea is to bring the ground images in the quasi-epipolar geometry basing on two steps: firstly estimating a preliminary rotation for both images, in order to align them in the disparity prevailing direction, then performing a refinement of the image relative position thanks to a set of well-identified TPs, exploited in order to correct possible residual roto-translations between master and slave images. These orthorectified images are referred to as GrEI.
Following the good results of previous works (such as Zhang & Fraser, 2008;Zhang & Gruen, 2006;Rong, Dazhao, Song, & Huiqin, 2008), a classical coarse-to-fine pyramidal scheme (Figure 1) is adopted to take advantage of an iterative solution based on multi-resolution imagery approach: intermediate DSMs at multiple resolutions are obtained and matches at low-resolution serve as approximations to restrict the search space. In propagating through the pyramid, matching continues until the highest desired resolution is reached.
Through the new procedure for GrEI generation, described in the following sections, a more accurate result in terms of GrEI is achieved and for the user there is no longer need to distinguish between acrosstrack or along-track images since the strategy has been generalised for being able to process images with whatever scanning direction.
The DSM generation processing chain identified consists of the following main steps hereafter illustrated: • images ground projection • GrEI generation • dense image matching and disparity map generation • DSM generation Figure 1. Image pyramidal approach: low resolutions serve as approximations to restrict the search space for the following levels.

Images ground projection
The first step implemented is the images ground projection: raw satellite images are projected in a ground geometry, thanks to rational polynomial coefficients (RPCs) provided by the vendors along with the images, using a coarse DSM (e.g. SRTM, ASTER GDEM, AW3D30) or other freely available DSMs for the area of interest.
Compared to classical approaches, here the quasiepipolarity is achieved in the object space, producing an easier procedure for DSM generation. In some previous works, such as (Haala, Stallmann, & Statter, 1998), a constant height plane is used for images projection, but in general it is better to know the elevation range of the area the scenes cover in order to improve the accuracy of the successively generated epipolar lines. Gyer (1981) suggested the use of a coarse DSM, and in more recent times also other authors proposed a projection in the object space, exploiting standard area-based matching techniques (e.g. Radhadevi et al., 2010).
Thanks to its modular structure, DATE is easily upgradable to upcoming sensors since it is sufficient to format the images metadata files in a standardised way in order to be able to manage images from new satellites. But not always satellite imagery are supplied with RPCs; when this does not happen, it is possible to generate the RPCs. SISAR, a software developed at the Geodesy and Geomatics Division, University of Rome "La Sapienza" (Capaldo, Crespi, Fratarcangeli, Nascetti, & Pieralice, 2011), among the other things, is able to generate RPCs starting from the rigorous sensor model, exploiting a terrain-independent approach.
If the orientations of all the images were fully consistent, so that also all the images relative orientation would not display global transversal parallax errors at the adopted resolution, the orthorectified images generated in this way could be considered as GrEI at each pyramidal level and would satisfy the quasi-epipolar condition. This is a general assumption, but usually an additional check, with possible transformation estimation, must be done for the actual GrEI generation. Furthermore, the scan direction has to be taken into account in order to correctly rotate the images to achieve the quasi-epipolar geometry.

GrEI generation
In order to carry out a dense matching, it is necessary to perform some transformations to bring the ground images in the quasi-epipolar geometry. In order to achieve a more accurate result in terms of GrEI, with respect to what is described in Di Rita et al. (2017), a twostep transformation is applied. At first, a preliminary approximate rotation is estimated for both images in order to align them in the disparity prevailing direction. This is performed according to forward and inverse rational polynomial functions (RPFs), on a square image grid 10 × 10. Consider a point on the grid P 1 in the master image, and project it by inverse RPF (imageto-ground transformation) at two heights h max and h min giving respectively P 2 0 and P 3 . Then P 3 is forward projected onto slave image giving P 2 ; P 2 is then inverse projected at h max giving P 0 2 ( Figure 2). The points P 0 1 and P 0 2 , both at height h max , are considered for the computation of the rotation angle that takes into account the orientation of the satellite scan direction. As a matter of fact the vector P 1 0 À P 0 2 define the quasi-epipolar direction in object space (being checked its being practically constant throughout the image) and thus the α angle (azimuth with respect to the north) to be applied to the images to act as GrEI ( Figure 3).
The applied azimuth is the mean of those obtained over the grid points. This approach works both with along and across track data, and no additional information is needed apart from RPCs. Then, as a second step, a refinement of the images relative position is performed through the model already defined in Di Rita et al. (2017). Using a set of well-identified TPs and according to Equation (1), a (small) rotation and a translation to be applied to the slave image, in order to correct possible residual roto-translations between master and slave images due to RPCs inaccuracy, are estimated: where x 0 , y 0 are TPs master coordinates; x, y are TPs slave coordinates; D is the disparity value in each point; θ is the mean rotation angle and Ty 0 is the y mean transversal parallax. The model is clearly non-linear and it should be linearised in order to solve it with a least squares approach (see Di Rita et al., 2017 for the analytics details).
Orthoimages are aligned in the disparity prevailing direction and, applying the model in Equation (1), a potential shift (mean rotation and transversal parallax), caused by systematic error in the image relative orientation based on metadata RPCs, is removed. As a matter of fact, the model takes into account that each TP can have a whatever disparity value on the x coordinates in order to return null y mean transversal parallax.
In particular, as illustrated in more detail in Di Rita et al. (2017), TPs are homogeneously identified on the orthoimages and filtered in order to keep only the most reliable points: 2% of the TPs is selected based on the higher feature descriptor values computed with the Brief Descriptor Extractor OpenCV algorithm.
This new approach, based on two subsequent transformations for GrEI generation, works both with along and across track data, and no additional information is needed apart from RPCs.

Dense image matching and disparity map generation
In the last years, methods first developed in computer vision, based on epipolar imagery and dense stereo matching (Scharstein & Szeliski, 2002;Hirschmüller, 2008;Krauss, Reinartz, Lehner, Schroeder, & Stilla, 2005;d'Angelo, Lehner, Krauss, Hoja, & Reinartz, 2008), have been more and more exploited (Reinartz, D'Angelo, Krauss, & Chaabouni-Chouayakh, 2010) to generate high-resolution DSMs. For a long time, local correlation-based matching methods dominated photogrammetry. Due to the size of the correlation window, resulting surface models have a lower spatial resolution than the input images, fine structures are lost and object boundaries are smoothed. In contrast, global methods permit pixel-wise matching, which is guided by a global smoothness constraint. The resulting DSMs have the same resolution as the input images, but global methods typically show a high processing time (Wohlfeil, Hirschmüller, Piltz, Börner, & Suppa, 2012). Since dense matching is typically slow, the search range should be limited: this can be done thanks to the initial projection using the a priori DSM, which allows to shrink the dense matching search range. The Semi-Global Matching (SGM) algorithm by Hirschmüller (2008), performing matching of corresponding pixels along the epipolar lines, provides a very good trade-off between accuracy and processing speed.
In order to compute a disparity map for the stereopair, at each pyramidal level the SGM algorithm, as implemented in the OpenCV library (OpenCV library overview, 2015), is used. This is actually the Semi-Global Block Matching (SGBM) algorithm, which differs from the original one for considering only five search directions for cost aggregation instead of eight (even if it is possible to run the full variant of the algorithm, yet consuming a lot of memory), for matching not individual pixels but blocks of pixels (even if changing the default parameter it is possible to reduce the blocks to single pixels) and for having implemented a simpler Birchfield-Tomasi sub-pixel metric (Birchfield & Tomasi, 1998) instead of a mutual information (MI) cost function (Viola & Wells, 1997). It must be recalled that, for its implementation, the OpenCV SGBM code only allows the use of 8-bit images instead of the 11-, 12-or 16-bit satellite imagery. This could lead to a compression of the range of radiometric values. Nevertheless, the majority of satellite images indeed present a non-uniform and mostly compressed radiometric distribution. The idea is to use the OpenCV SGBM having preliminarily reduced the bit depth considering a 95th percentile and successively scaling the image radiometry between 0 and 255. For the majority of the satellite images, this leads to a negligible loss of radiometric information.
Note that using GrEI the disparity maps obtained at each pyramidal level are in the same reference frame of the a priori DSM, and the local pixel disparity values are proportional to the height correction that will be applied to the a priori DSM.
In Figure 4 an example of a disparity map computed from Pléiades imagery on the Trento test site is reported. Some unmatched areas (black pixels) are visible due to occlusions and shadows: as a matter of fact, in general, DSMs, which are automatically generated from stereo satellite image matching procedures, tend to have many null locations so-called voids or holes since occlusions, especially in urban scenes, result in some deficiencies in the stereo matching phase (Bafghi, Tian, d'Angelo, & Reinartz, 2016). At the moment, in the no-matched areas, the SRTM value is left unchanged without summing up any height correction.

From pixel disparity to DSM generation
Generally, to obtain a DSM from stereo image pairs, conjugate points are firstly automatically matched on the images; then, these points are mapped by forward intersection into the object space and the DSM is retrieved from these points by triangulation and interpolation. Since here disparities are computed in the object space, the projection step can be avoided: in order to get the final DSM, it is sufficient to sum the found disparities to the coarse DSM exploited as a reference for the initial projection. But first, the obtained local pixel disparity values have to be converted into height corrections with respect to the used a priori DSM, through an efficient procedure based on the computation of a disparity-to-height conversion factor (C factor ) exploiting RPCs-based orientation model ( Figure 5).
In particular, a constant C factor is computed starting from the metadata RPCs, remarkably improving the procedure explained in Di Rita et al. (2017) and exploiting the procedure for the preliminary rotation angle computation, aimed at the images alignment in the disparity prevailing direction for epipolar geometry achievement (see section "GrEI generation"). As a matter of fact, the C factor is obtained dividing the magnitude of the vector P 0 1 À P 0 2 (see Figure 3) to the height variation (h maxh min ): the final C factor used to transform disparities to height corrections is just the mean of the values obtained over the grid where C factor is the disparity-to-height conversion factor; ΔE k is the difference in the object space (at h max ) in east coordinates between two given points; ΔN k is the difference in the object space (at h max ) in north coordinates between two given points and Δh k is the height difference considered (h maxh min ). It has been checked that the height variation set (h maxh min ) has no appreciable influence on the final C factor value nor on the preliminary rotation angle computation. As a matter of fact, both using a small height variation (e.g. 50 m) or the height variation of the entire tile (here about 2000 m), deductible from the coarse DSM, the final C factor value does not change significantly.
The obtained height corrections are then summed up to the height of the a priori DSM at each pyramidal-level iteration. As last step of the deployed workflow, the geocoded DSM is saved in a standard geotiff format, exploiting OSSIM file handling functionalities.

DATE testing
The results obtained with Pléiades-HR (Pléaides infopage 2015) satellite optical imagery are reported and analysed. To show and evaluate DATE behaviour and performance, DSMs in areas with different land use/ land cover and morphology and in the whole image footprint have been generated and evaluated.
In order to assess the obtained results, a 2.5D comparison with respect to a reference DSM was performed. The reference data set used to assess the generated DSMs in this work is a LiDAR DSM with grid posting 1.0 × 1.0 m and a mean elevation accuracy of 0.25 m (Figure 6), freely available on the website of "Provincia Autonoma di Trento" (http:// www.territorio.provincia.tn.it).

Analysed satellite imagery dataset
The analysed data set consists of three along-track Pléiades-HR optical images referring to Trento test site, Northern Italy.
Located in the northeast part of Italy, the test site lies in a valley at an average altitude around 200 m, surrounded by very steep slopes and Alps peaks up to 2000 m. It includes urban areas with residential, industrial and commercial buildings with different sizes and heights, agricultural and forested areas, offering therefore a very heterogeneous landscape in terms of morphological complexity, land use and land cover (Agugiaro, Poli, & Remondino, 2012).
The images composing the triplet, hereafter called image 01, image 02 and image 03, were acquired in August 2012 in north-south scan direction. The average incidence angles of the three images are, respectively, 13 , À 13 and À 19 in along-track direction with respect to the nadir and close to zero in across-track direction. Thanks to the fast scanning, the illumination conditions are almost constant during the acquisitions. The data include metadata files together with their original RPCs. The three images cover an area of about 390 km 2 on the Trento test site, which already includes images from other very high resolution (VHR) satellite sensors. The native mean ground sample distance (GSD) varies between 0.72 and 0.78 m depending on the viewing direction, even if the images are oversampled and supplied to the user with a 0.50 m GSD.
The three images will be analysed by pairs, testing all possible combinations, this is: • pair 1, composed by images 01 and 03 • pair 2, composed by images 01 and 02  • pair 3, composed by images 02 and 03 In Table 1, the main features of the optical data set are reported. Furthermore, in order to assess the potential of the tri-stereo acquisition, a merging procedure has been adopted and its results are presented in this section.

Results analysis and DSMs validation
The height differences between the DSM generated with DATE and the LiDAR DSM were computed together with their statistics [root mean square error (RMSE), standard deviation, bias and LE95]. These statistics used for the DSMs assessment are those usually adopted for the assessment of the global DSMs (as SRTM, ASTER, TANDEM-X DSMs) (Rodriguez, Morris, & Belz, 2006;, Gesch, Oimoen, & Evans, 2014).
The RPCs provided by the vendors along with the images have been corrected in order to remove any potential orientation error: a set of ground control points (GCPs) has been used and the RPCs offset corrected considering the image coordinates differences obtained in sample and line directions (see Table 2). In Table 3 are illustrated the results obtained with the corrected RPCs on the GCPs used. Furthermore, in order to check the reliability of the RPCs bias correction, a set of well-distributed check points (CPs) has been used. The image (here the results for image 02 are reported) is firstly projected on the LiDAR DSM, using both the original and the corrected RPCs, generating two orthoimages; then, visible features have been detected and their ground coordinates on the two orthoimages compared to the corresponding coordinates on the LiDAR DSM. A planimetric assessment has been carried out (Table 4), checking, for each CPs, the differences between the LiDAR coordinates and the orthoimage coordinates (generated both with original and corrected RPCs) in east and north direction: with the corrected RPCs, for the east coordinates the standard deviation is reduced and for the north coordinates a sharp improvement in the mean is evidenced.
Different land covers and morphologies have been investigated. In particular, the entire tile, an urban, a mainly flat (even if in this scene it is difficult to find a wide flat area) and a mountainous areas are analysed (see Figure 6). The urban tile is characterised by small adjacent buildings and narrow streets, whereas the mountainous one by very steep slopes.
In Table 5 the results for each Pléiades DSM for the different areas analysed are reported, together with the statistics for the merged DSM. As a matter of fact, a merging procedure has been applied in order to exploit the redundant information from the tri-stereo acquisition. The three generated DSMs have been merged, giving them a different weight percentage on the basis of the images intersection angles: a smaller intersection angle implies a weaker acquisition geometry (and thus a lower reliability in terms of reconstruction accuracy), whereas a larger one entails a stronger geometry (and a higher reliability in terms of reconstruction accuracy). On the contrary, pairs with small intersection angles can give a significant contribution to the DSM completeness since occlusions are reduced, especially in the urban area.
In general, a better accuracy is achieved in the flat area, whereas the mountain and the urban areas behave   pretty much in the same way. Moreover, it is evident that DSM from pair 3 shows always a very high bias with respect to the other two pairs. However, it is quite clear that with the merging procedure an improvement both in the statistical parameters and in the completeness level (see Table 6) of the final DSM is achieved. With respect to the three DSMs, the Pléiades pair 3 shows the higher degree of completeness (>99%), while for the other two DSMs it is substantially the same, varying between 91% and 94%: this is due to the fact that pair 3 offers the lowest intersection angle and hence the dense matching performs much better. The merged DSM shows a completeness percentage slightly greater than DSM from pair 3. It can be seen that, in this acquisition configuration with very different angles of intersection between images, the triplet exploitation gives very good results: as a matter of fact, the merged DSM takes advantage both of the DSM from pair 3, with a small intersection angle, thanks to which the dense matching gives few occluded areas, and of the DSMs with a larger intersection angle, that guarantee an acquisition geometry favourable for a good photogrammetric reconstruction, and hence a good final accuracy.
Also by visual inspection, an improvement of the reconstruction thanks to the merging procedure can be detected. An interesting example of this can be seen in the urban area, where a slightly higher reconstruction completeness and a better building reconstruction is observed (see Figures 7 and 8).
Thanks to an overview, it is possible to check that the terrain profile appears well modelled, even in very steep areas, for each of the generated DSMs, with different completeness level (see Figures 9-Figure 12).
The height discrepancies histograms of the generated DSMs are reported in Figures 13 and 14. It is evident that the histograms of the entire area and of the mountain tile present an asymmetry towards the negative zone: probably the entire tile behaviour is affected by the mountain area, which is the pre-eminent in the scene analysed. Presumably this negative trend is due to the different trees height in the forest between LiDAR and Pléiades imagery acquisitions. Furthermore, as far as concerns the urban tile, the DSM from pair 3 shows a more Gaussian shape than the other two. As regards the flat tile, it is evident that the DSM from pair 3 displays a slightly bimodal behaviour, probably due to the small images intersection angle, which undermines a proper 3D reconstruction.
Some investigations to evaluate VHR optical data sets on Trento test site have been carried out in recent years (Agugiaro et al., 2012;Poli, Remondino, Angiuli, & Agugiaro, 2013): the results here obtained outperform by about 10-20% those presented in these works, and this confirms the quality of the developed workflow. However, it is necessary to underline that DATE shows high performances in terms of efficiency: as a matter of fact, as illustrated in Di Rita et al. (2017), it is about three times faster than well-known commercial software.

Conclusions and future prospects
A new procedure for GrEI generation from satellite optical imagery, within DATE FOSS4G, is presented. The main purpose was to overcome the issues related to the epipolarity resampling for satellite images, for which it is not straightforward the epipolar geometry achievement, due their multiple projection centres.
The strategy for the epipolar resampling for satellite imagery relies on a preliminary ground projection using an a priori DSM (such as SRTM, a coarse and freely available global DSM), thanks to which the search space for the successive dense matching is significantly restricted. Then, an approximate rotation is computed through RPCs in order to align the images in the disparity prevailing direction. Afterwards, a roto-translation is estimated, thanks to a set of well-identified TPs, and applied, to refine the images relative orientation and to achieve the epipolar geometry: these orthorectified images can act as GrEI and can undergo a dense image matching procedure. Furthermore, in order to be able to keep the transversal parallax errors well below the imagery resolution, a coarse-to-fine pyramidal scheme is adopted: this allows to iteratively refine the a priori DSM until the desired final resolution is achieved. The final DSM at each pyramidal level is the input for the next pyramidal level, acting as the a priori DSM. In this way, at lower resolution it is possible to detect Note that the bias is positive when the reference LiDAR height is above the extracted DSM. larger structures whereas at higher resolution small details are progressively added to the already obtained DSM. Apart from this original procedure to obtain the quasi-epipolar imagery, there are other two DATE key features which have not changed with respect to the previous version: the use of computer vision algorithms in order to improve the processing efficiency and make the DSMs generation process fully automatic; the free and open-source aspect of the developed code (freely available at the URL of the GitHub repository, https://github.com/martidi/ opencv dsm/tree/imageStack).
The potential of the developed strategy has been demonstrated through application to DSMs generation from Pléiades-HR tri-stereo imagery in the Trento test site. Accuracies achieved are in the order of 4 m, with a completeness level greater than 99%, for the entire tile.
The urban and the mountain tiles present similar standard deviation (a little more than 4 m) and a bias respectively of 0.3 and −0.5 m, whereas the flat tile shows a lower standard deviation (1.5 m) and a bias of about 0.3 m. These results, in terms of RMSE, are better by about 10-20% compared to those obtained in previous works (Agugiaro et al., 2012;Poli et al., 2013).
The good results obtained with the merging procedure encourage, as a future work, to develop a merging functionality able to preserve edges. This implies that it could be worthwhile to improve DATE capability for tri-stereo processing in order to obtain more accurate and complete DSM starting from the pair-derived DSMs.
Moreover, the MI could be implemented as a similarity criterion for the dense image matching step (see section "Dense image matching and disparity map generation") since, especially for images with Urban area DSM pair 1 (a) Urban area DSM pair 2 (b) Urban area DSM pair 3 (c) Urban area merged DSM (d) Figure 7. DSMs generated in the urban area with pair 1 imagery (a), pair 2 imagery (b), pair 3 imagery (c) and with the merging procedure (d).
Error map for the urban area pair 1 (a) Error map for the urban area pair 2 (b) Error map for the urban area pair 3 (c) Error map for the urban area merged DSM (d) Figure 8. Error maps, with respect to the LiDAR DSM, generated for the urban tile with pair 1 imagery (a), pair 2 imagery (b), pair 3 imagery (c) and with the merging procedure (d). Figure 9. DSM obtained from Pléiades-HR imagerypair 1. Figure 10. DSM obtained from Pléiades-HR imagerypair 2.  a large stereo baseline, which show stronger radiometric differences, the Birchfield-Tomasi metric might not be the most appropriate similarity. Besides, as a further development, a solution to fill gaps for the unmatched areas both through interpolation and fusion with other information will be investigated.