## 3 DEM Extraction  
In practice, a 3D reconstruction from a set of images can be performed in either the image space or the object space. The image space approach is mostly used with stereo pairs or triplets of images. The object space approach is commonly used with multiple overlapping images, such as images acquired with drone and close range photogrammetry.  
In this study, we adopt the object space approach, as it is well adapted to process images from push frame systems like those delivered by SkySat and PlanetScope. We also use the image-space method to a subset of data, using tri-stereo SkySat images, for comparison.  
The main steps for 3D extraction in image space are (Figure 4):  
1. Epipolar rectification: it consists in resampling stereo pairs based on the adjusted RFM, so that the two images have a common orientation and the matching features between the images appear along a common axis [27].
2. Stereo-matching: it consists in computing the correspondences between pixels of the image pairs. These correspondences are computed using a correlation technique (e.g., NCC, FFT) or using a Semi-global matching scheme [28]. Results are displayed as disparity maps.
3. Disparity maps fusion: intermediate results generated from each possible stereo are merged to produce a final DSM map. The fusion is performed using local approaches (e.g., mean, median) or global optimization (e.g., total variation, gradient decent) [26]. In the object space (Figures 4 and 5), the main steps are:  
1. Multi-image matching: an object-based matching algorithm, e.g. OSGM [29], is applied directly in the object space, hence the epipolar rectification is no longer necessary; the transformation between object space and image space relies on the refined RFMs.
2. Spatial forward intersection: this leads directly to dense 3D point cloud.
3. Meshing: it consists in deriving 3D surfaces by interpolating the dense point cloud.
4. Mesh-based DEM: gridded terrain model (i.e., 2.5D raster map) is derived from the 3D mesh. In addition, for rendering, a 3D-textured model is generated using (see Supplementary Materials):  
- The 3D mesh model, which has rich geometric information;  
- The image collection (RGB spectral bands), which provides high photorealistic details about the texture of the objects;  
- The free-bias RFMs.  
SkySat and PlanetScope images have a high-dynamic range. This makes it difficult to visualize and compress the full range for display, which can result in loss of details. Therefore, before generating the texture map, a contrast enhancement is performed. We apply Contrast Limited Adaptive Histogram Equalization (CLAHE) in order to amplify the image contrast while preserving details [30]. The operator is also implemented in the COSI-Corr RFM package. Each image is divided into small tiles; then the histogram is equalized by tile. The radiometric noise gets amplified during this process. Therefore, contrast limiting is applied to limit the contrast below a specific limit threshold (\(th_{climit}\)). Then, if any histogram bin is above the specified contrast limit (\(th_{climit}\)), those pixels are clipped and distributed uniformly to other bins before applying histogram equalization. Ultimately, to remove artifacts in tile borders, bilinear interpolation is applied. The 3D-textured models are used to generate the flyover animations provided in Supplementary Materials.  
Figure 4: Schematic workflow for SkySat image geometry refinement and 3D extraction.