Shading aware DSM generation from high resolution multi-view satellite images

ABSTRACT

In many cases, the Digital Surface Models (DSMs) and Digital Elevation Models (DEMs) are obtained with Light Detection and Ranging (LiDAR) or stereo matching. As an active method, LiDAR is very accurate but expensive, thus often limiting its use in small-scale acquisition. Stereo matching is suitable for large-scale acquisition of terrain information as the increase of satellite stereo sensors. However, underperformance of stereo matching easily occurs in textureless areas. Accordingly, this study proposed a Shading Aware DSM GEneration Method (SADGE) with high resolution multi-view satellite images. Considering the complementarity of stereo matching and Shape from Shading (SfS), SADGE combines the advantage of stereo matching and SfS technique. First, an improved Semi-Global Matching (SGM) technique is used to generate an initial surface expressed by a DSM; then, it is refined by optimizing the objective function which modeled the imaging process with the illumination, surface albedo, and normal object surface. Different from the existing shading-based DEM refinement or generation method, no information about the illumination or the viewing angle is needed while concave/convex ambiguity can be avoided as multi-view images are utilized. Experiments with ZiYuan-3 and GaoFen-7 images show that the proposed method can generate higher accuracy DSM (12.5–56.3% improvement) with sound overall shape and temporarily detailed surface compared with a software solution (SURE) for multi-view stereo.

KEYWORDS:

• Shape from Shading (SfS)
• multi-view stereo
• Digital Surface Model (DSM)
• high resolution multi-view satellite images

1. Introduction

The Digital Surface Models (DSMs) record the elevation values of the regularly spaced object surface. They are the base data for many applications, such as hydrologic modeling, urban planning, route selection, and earth volume calculation, to name a few. DSMs with higher spatial resolution and accuracy are needed considering the applications above. Generally, DSM can be generated with the Light Detection and Ranging (LiDAR) technique or photogrammetry method using multi-view optical images (Gabet, Giraudon, and Renouard 1997; Toutin 2004). The latter is becoming popular in large-scale DSM generation because of its low cost, the increase in high-resolution satellite images, and their extensive applications. Typically, the photogrammetry method generates the elevation values through image matching which finds the corresponding pixels between images according to the similarities of image texture (Shan et al. 2020). When the image texture is rich, image matching performs well. However, image matching may underperform in textureless areas. Consequently, the detailed shape in the textureless area cannot be generally reconstructed effectively with the image matching technique (Furukawa and Hernández 2015).

By contrast, Shape from Shading (SfS) (Frankot and Chellappa 1988; Horn 1970; Horn and Brooks 1986; Lions, Rouy, and Tourin 1993) can recover the detailed shape with one or more images even in a textureless area. Thus, it is complementary to the image matching technique. Some studies focus on the refinement of existing low-resolution DEM/DSM using SfS technique with single or multiple images (Heipke, Piechullek, and Ebner 2001; Lohse, Heipke, and Kirk 2006; Peng, Zhang, and Shan 2015). However, the concave/convex ambiguities can be obvious when using SfS with a single image which may cause wrong shape recovery, and no shading-based DSM generation method exists directly from multi-view satellite images. In this paper, we proposed a Shading Aware DSM GEneration Method (SADGE) for high-resolution multi-view satellite images. First, an improved Semi-Global Matching (SGM) (Hirschmüller 2007) technique is used to generate an initial DSM. Then, the imaging process is modeled as an interaction of the illumination, surface albedo, and surface normal while introducing an objective function that considers the multi-view images and all the parameters. At last, the final DSM is generated by optimizing the objective function. The experiments with ZiYuan-3 (Tang et al. 2013; Li, Wang, and Jiang 2021) and GaoFen-7 (Xie et al. 2020) stereo images show that our method can generate high accuracy DSM with fine details.

The paper is organized as follows: Section 2 briefly introduces the related work. Section 3 describes the proposed method in detail, while section 4 presents the quantitative and qualitative experiments. Section 5 concludes our work.

2. Related work

Image matching means determining the corresponding pixels of the same point in two or more images. Usually, image similarity measurements, such as the normalized correlation coefficients, Mutual Information (MI), Census, among others, are used as matching costs to judge photo consistency. According to the output scene representations, image matching methods can be roughly divided into four types (Furukawa and Ponce 2010): depth map based (Goesele, Curless, and Seitz 2006; Kolmogorov and Zabih 2001; Kwon, Tai, and Lin 2015), voxel based (Faugeras and Keriven 1998; Heise et al. 2015), patch based (Gallup, Frahm, and Pollefeys 2010; Habbecke and Kobbelt 2006; Hou et al. 2018), and mesh based (Esteban and Schmitt 2003; Furukawa et al. 2009). It can reconstruct highly detailed surface in a textured area (Hu et al. 2021; Liu and Ji 2020; Ma et al. 2021). However, it often performs poorly in a textureless area where accurately measuring photo consistency is difficult.

Different from image matching technique, SfS technique extracts the surface normal to recover the detailed surface by modeling the imaging process with the illumination, surface albedo, and surface normal. In the early days, SfS technique is mostly developed in the reconstruction of synthetic objects and planetary surfaces. Over few decades’ improvement, SfS technique has been used in many other scenes (e.g. human faces (Kemelmacher-Shlizerman and Basri 2010), close-range objects (Panagopoulos, Hadap, and Samaras 2012) and Earth surface (Peng, Zhang, and Shan 2015)). Most SfS-related algorithms are based on one single image. However, SfS is an ill-posed problem given that the illumination, surface albedo, and surface normal are all unknown while the image is the only observation, especially when only one image is available. Therefore, assumptions about the illumination and surface albedo are often imposed to make the problem solvable. For example, Ikeuchi and Horn (1981) assumed uniform surface albedo and known illumination. Introduction of the spherical harmonics (Basri and Jacobs 2003) also eases the limits of illumination, but it assumes far and point light source. Another method called Photometric Stereo (PS) (Woodham 1980) aims to increase the number of illuminations, i.e. capturing images from a fixed position under different illuminations. It can produce compelling results, but capturing images is difficult from a fixed position with different illuminations in natural scenes which prevents its application range.

As mentioned above, image matching technique can effectively reconstruct the overall shape in a well-textured area, while SfS technique can recover the detailed shape even in a textureless area. Thus, it is of great interest to combine these two techniques and generate detailed surface in high accuracy. One method (Liu, Wu, and Wöhler 2018; Kim, Torii, and Okutomi 2016; Peng, Zhang, and Shan 2015; Wu et al. 2018, 2011; Xu et al. 2018) used the SfS technique to refine the existing surface derived from global open terrain data, such as the Shuttle Radar Topography Mission (SRTM) and the ALOS Global Digital Surface Model “ALOS World 3D − 30 m” (AW3D30) or generated through image matching or other methods, such as InSAR. The other method (Langguth et al. 2016; Maurer et al. 2018; Quéau et al. 2017) attempted to fuse image matching and SfS into one model and then solved it together with close-range images. The former actually adopted a two-stage strategy, first initial surface reconstruction with image matching, InSAR or LiDAR, and then SfS-based surface refinement. The latter combined image matching and SfS into one stage. However, most of the algorithms considered the planetary geometry (Lohse, Heipke, and Kirk 2006) or the close-range surface, such as faces and statues (Langguth et al. 2016). For DEM on earth, the work most related to our method is that of Peng, Zhang, and Shan (2015) who refined the SRTM GL3 of 90 m resolution to 30 m with one single Landsat image and subsequently achieved fine results, but the concave/convex ambiguities persisted which might cause wrong surface detail recovery. In addition, the illumination direction is assumed to be known, and uniform surface albedo is assumed in a local area in their algorithm.

In this paper, we proposed a method that aims to accurately generate a detailed DEM with high resolution from a high-resolution multi-view satellite stereo image. No need arises for illumination direction in our algorithm because spherical harmonics is used to represent the illuminations, and the surface albedo can be different in a different point.

3. Proposed method

This section rigorously introduces the proposed method. As shown in Figure 1, an initial surface is firstly reconstructed with the multi-measurement Semi-Global Matching (mSGM). Then, an objective function considering the illumination, surface albedo, and surface normal is built. Finally, the initial surface reconstructed earlier is refined by optimizing the objective function. Different from the existing DEM refinement algorithms, the refinement procedure in the proposed method used multiple images; moreover, eliminating the convex/concave ambiguity is helpful (Quéau et al. 2017).

Figure 1. An overview of the proposed SADGE: (a) the input multiple satellite images, (b) the initial DSM reconstructed with mSGM, (c) the final DSM after shading based refinement.

3.1 Image matching: mSGM

For the image matching technique, we followed Tao’s work (Tao 2016) and applied mSGM to obtain the initial DEM with high-resolution satellite images. Similar to the traditional SGM, mSGM applied pyramid matching in different image layers. However, mSGM combined two kinds of matching cost, namely, MI and Census, to take advantage of them. Census, as a window-based matching cost, is very robust. It can measure the similarities of two images without needing a prior probability lookup table required by MI. MI can calculate the similarities of pixel pairs instead of two image windows. Thus, it has a huge advantage in capturing the structure details of objects. In the top layer, the matching measure of SGM has been changed from MI to Census to increase the robustness of SGM. Except for the top layer, MI is still used as matching measures in the other layers. In addition, mSGM adopts the coarse-to-fine strategy, and matching results are passed from top to bottom. The matched disparity map of the current image layer is used to interpolate the initial disparity map of the next image layer, and the range of the disparity candidates for each pixel is estimated. The matching efficiency is largely improved, and the possibility of mismatch is reduced.

3.2 Refinement of DSM

As for the refinement of DEM, considering the reflection model as the first step is necessary because the SfS technique is essentially an inverse of imaging process. The Lambertian reflection model that reflects the same in all directions is used similar to many other DEM generations or refinement algorithms. However, different from those algorithms on Earth or other planets, the spherical harmonics is utilized instead of the classic “direction model” that needs the incident angle from the light source. The spherical harmonics is first introduced by Basri and Jacobs (Basri and Jacobs 2001) and is widely used in close-range surface refinement. They proved that the images of the convex Lambertian object can be accurately approximated with a low-dimensional linear subspace. However, a prerequisite is set, i.e. the light source should be distant which means the complex illumination or close light source might influence the stability of the spherical harmonics. Fortunately, the light source of the surface on Earth in satellite images is mainly from the sun and the sky, both of which are sufficiently far. Therefore, second-order spherical harmonics basis functions are temporarily used to represent the illumination in the proposed method considering the complexity. Based on the spherical harmonics mentioned above, the imaging process can be modeled as

(1) $I_{i,j}^m={\rho }_{\mathit{i},\mathit{j}}\cdot \left[\begin{array}{ccc}{L}_0^m& \dots & L_8^m\end{array}\right]\cdot {\left[\begin{array}{ccccccccc}1& n_x& n_y& n_z& n_x{n}_y& n_x{n}_z& n_y{n}_z& n_x^2-n_y^2& 3n_z^2-1\end{array}\right]}^{\mathrm{T}}$(1)

where $I_{i,j}^m$ is the corresponding pixel value in the $\mathit{m}$th image, ${\rho }_{\mathit{i},\mathit{j}}$ is the albedo, $\left(n_x,n_y,n_z\right)$ is the unit normal vector of the surface point $\left(\mathit{i},\mathit{j}\right),$ and $L_0^m$ to $L_8^m$ are the coefficients of the spherical harmonics basis of the $\mathit{m}$th image. For the same surface point $\left(\mathit{i},\mathit{j}\right)$, the surface albedo ${\rho }_{\mathit{i},\mathit{j}}$ and normal $\left(n_x,n_y,n_z\right)$ are the same for different images, whereas the pixel value $I_{i,j}^m$ and spherical harmonics basis $L_0^m$ to $L_8^m$ are different in different images. In addition, different surface points can have different surface albedo which is closer to the real situation.

Considering the imaging process above, an objective function is built:

(2) $\underset{d_z,\mathit{\rho },\mathit{L}}{min}{E}_{\mathit{data}}\left(d_z,\mathit{\rho },\mathit{L}\right)+\mathit{\alpha }{E}_{\mathit{geometry}}\left(d_z\right)+\mathit{\beta }{E}_{\mathit{albedo}}\left(\mathit{\rho }\right)$(2)

where $d_z$ is the change in the Z-axis coordinates of the surface point, $\mathit{\rho }$ is the albedo of the surface point, and $\mathit{L}$ is the coefficient of the spherical harmonics basis. $\mathit{\alpha }$ and $\mathit{\beta }$ are the weights to balance the data term $E_{\mathit{data}}$, geometry smooth term $E_{\mathit{geometry}}$, and reflectance smooth term $E_{\mathit{albedo}}$. Next, the three terms will be introduced in detail.

The data term $E_{\mathit{data}}$ is used to measure the difference between the observation, i.e. the image value and the estimator, i.e. the inverse rendered pixel value.

(3) $E_{\mathit{data}}\left(d_z,\mathit{\rho },\mathit{L}\right)=\sum _{\mathit{m}=0}^M\sum _{\left(\mathit{i},\mathit{j}\right)\in S}\frac{I_{i,j}^m\left(\mathit{Z}+d_z\right)-R_{i,j}^m\left(\mathit{Z}+d_z,\mathit{\rho },\mathit{L}\right)}{\mathit{M}}$(3)

where $\mathit{M}$ is the number of the visible satellite images, $\mathit{S}$ is the initial surface, i.e. the DEM reconstructed by SGM, $\left(\mathit{i},\mathit{j}\right)$ denotes the horizontal coordinates of the surface point, $\mathit{Z}$ is the vertical coordinate of the surface point, $I_{i,j}^m$ is the corresponding image pixel value of the surface point $\left(\mathit{i},\mathit{j}\right)$, and $R_{i,j}^m$ is the inverse rendered pixel value of the surface point $\left(\mathit{i},\mathit{j}\right)$. Notably, the satellite images used in this paper are captured from distant places, and the view angles are small; thus, the surface point will be considered visible as long as it is in the range of the image. In addition, $I_{i,j}^m$ will be re-computed with the change in $\mathit{Z}+d_z$ as well as $R_{i,j}^m$ which is defined in Equation (1)(1) $I_{i,j}^m={\rho }_{\mathit{i},\mathit{j}}\cdot \left[\begin{array}{ccc}{L}_0^m& \dots & L_8^m\end{array}\right]\cdot {\left[\begin{array}{ccccccccc}1& n_x& n_y& n_z& n_x{n}_y& n_x{n}_z& n_y{n}_z& n_x^2-n_y^2& 3n_z^2-1\end{array}\right]}^{\mathrm{T}}$(1) with the illumination, surface albedo, and surface normal. The data term $E_{\mathit{data}}$ is set following the principle that the inverse rendered pixel value should be equal to the image pixel value.

To keep the surface smooth and eliminate the obvious error from the SGM, the geometry smooth term $E_{\mathit{geometry}}$ is set.

(4) $E_{\mathit{geometry}}\left(d_z\right)=\sum \underset{\left(\mathit{i},\mathit{j}\right)\in S}{\mathit{\hspace{0.17em}}}\mathit{\hspace{0.17em}}\left(\sum \underset{\left(i^{\prime },j^{\prime }\right)\in A_{\mathit{ij}}}{\mathit{\hspace{0.17em}}}{w}_{i^{\prime }{j}^{\prime }}\frac{∥{\left(Z+d_z\right)}_{\mathit{ij}}-{\left(Z+d_z\right)}_{i^{\prime }{j}^{\prime }}∥}{\mathit{gsd}\cdot N_{{\mathit{A}}_{\mathit{ij}}}}+\mathit{\gamma }∥{\left(Z+d_z\right)}_{\mathit{ij}}-G_{\mathit{ij}}∥\right)$(4)

where $A_{\mathit{ij}}$ is the neighbor of the surface point $\left(\mathit{i},\mathit{j}\right)$; $\left(i^{\prime },j^{\prime }\right)$ is the neighbor point belonging to $A_{\mathit{ij}}$; $w_{i^{\prime }{j}^{\prime }}$ is the weight of the neighbor point $\left(i^{\prime },j^{\prime }\right)$; $\mathit{gsd}$ is the resolution of the DEM; $N_{{\mathit{A}}_{\mathit{ij}}}$ is the number of the neighbor points in $A_{\mathit{ij}}$ which is 8, considering that a 3 × 3 window is chosen as the neighbor in our experiment; ${\left(Z+d_z\right)}_{\mathit{ij}}$ and ${\left(Z+d_z\right)}_{\mathit{ij}}$ are the Z-axis coordinate of the surface point $\left(\mathit{i},\mathit{j}\right)$ and its neighbor point $\left(i^{\prime },j^{\prime }\right),$ respectively; $G_{\mathit{ij}}$ is the filtered vertical coordinate of the surface point $\left(\mathit{i},\mathit{j}\right)$ with guided filter; and $\mathit{\gamma }$ is the weight. The first term is set to smooth the surface while keeping the sharp edge, and the second term is set to eliminate the possible obvious error while considering the image edge similar to the first term. $w_{i^{\prime }{j}^{\prime }}$ and $G_{\mathit{ij}}$ are calculated using a guided filter with the image as guide to filter the DEM. As mentioned above, a 3 × 3 window is used as neighbor for the surface point. However, the window size is not big enough when the obvious error is not a point but within an area, and the increase in the window size means memory and complexity increase. Thus, the second term is set and calculated with large window size, i.e. 10 × 10 in our experiment. Accordingly, the second term in the geometry smooth term can eliminate the obvious error in the initial DEM, whereas the first term can keep the sharp edges in the image. The calculation of $w_{i^{\prime }{j}^{\prime }}$ is shown below:

(5) $w_{\mathit{i}{\mathit{ }}^{\mathrm{\prime }}\mathit{j}{ }^{\mathrm{\prime }}}=\frac{1}{{{\mathit{N}}_{{\mathit{A}}_{\mathit{ij}}}}^2}\sum _{\left(\mathit{i}{\mathit{ }}^{\mathrm{\prime }},\mathit{j}{ }^{\mathrm{\prime }}\right),\left(\mathit{i},\mathit{j}\right)\in A_{\mathit{ij}}}\left[1+\frac{\left({\mathit{I}}_{\mathit{ij}}-\mathit{\mu }\right)\left(I_{\mathit{i}{\mathit{ }}^{\mathrm{\prime }}\mathit{j}{ }^{\mathrm{\prime }}}-\mathit{\mu }\right)}{{\sigma }^2+\mathit{\epsilon }}\right]$(5)

where $\mathit{\epsilon }$ is a regularization parameter, $\mathit{\mu }$ and ${\sigma }^2$ are the mean and variance of the image value in $A_{\mathit{ij}}$, respectively.

As for the surface albedo $\mathit{\rho }$ in Equation (2)(2) $\underset{d_z,\mathit{\rho },\mathit{L}}{min}{E}_{\mathit{data}}\left(d_z,\mathit{\rho },\mathit{L}\right)+\mathit{\alpha }{E}_{\mathit{geometry}}\left(d_z\right)+\mathit{\beta }{E}_{\mathit{albedo}}\left(\mathit{\rho }\right)$(2) , it is allowed to be different for a different surface point. However, $\mathit{\rho }$ is multiplied with the spherical harmonics as shown in Equation (1)(1) $I_{i,j}^m={\rho }_{\mathit{i},\mathit{j}}\cdot \left[\begin{array}{ccc}{L}_0^m& \dots & L_8^m\end{array}\right]\cdot {\left[\begin{array}{ccccccccc}1& n_x& n_y& n_z& n_x{n}_y& n_x{n}_z& n_y{n}_z& n_x^2-n_y^2& 3n_z^2-1\end{array}\right]}^{\mathrm{T}}$(1) . To better separate the surface albedo $\mathit{\rho }$ with the spherical harmonics, an albedo smooth term is set:

(6) $E_{\mathit{albedo}}\left(\mathit{\rho }\right)=\sum _{\left(\mathit{i},\mathit{j}\right)\in S}\sum _{\left(i^{\prime },j^{\prime }\right)\in A_{\mathit{ij}}}{e}^{-k{\left(I_{ij}^0-I_{i^{\prime }{j}^{\prime }}^0\right)}^2}\frac{∥{\mathit{\rho }}_{\mathit{ij}}-{\rho }_{i^{\prime }{j}^{\prime }}∥}{N_{{\mathit{A}}_{\mathit{ij}}}}$(6)

where $\mathit{k}$ is a constant value, $I_{ij}^0$ is the mean pixel value of all the visible images for surface point $\left(\mathit{i},\mathit{j}\right)$, $I_{i^{\prime }{j}^{\prime }}^0$ is the mean pixel value of all the visible images for surface point $\left(i^{\prime },j^{\prime }\right)$, ${\rho }_{\mathit{ij}}$ and ${\rho }_{i^{\prime }{j}^{\prime }}$ are the albedo of the surface points $\left(\mathit{i},\mathit{j}\right)$ and $\left(i^{\prime },j^{\prime }\right)$, respectively.

After the objective function is built, the Levenberg-Marquardt implemented in the Ceres Solver (Agarwal and Mierle 2018) is used to optimize it, i.e. the best coefficients of the spherical harmonic basis, surface albedo, and normal are determined under the constraints of geometry and albedo smoothness.

4. Experiment and discussion

Two sets of experiments are designed to verify the effectiveness of the proposed method. The first experiment is conducted with the ZiYuan-3 stereo images in a desert area. However, no ground control points or high accuracy DEM exists for quantitative evaluation. Accordingly, a qualitative evaluation is made. As shown in Figure 1(b), image matching method like mSGM can effectively reconstruct the overall shape. However, when examining the details in the DSM generated by mSGM, it is noisy with several obvious errors. On the contrast, the proposed method which combines image matching and SfS technique can effectively reconstruct the overall shape while recovering the detailed shape and eliminating some obvious errors.

The second experiment is conducted with the GaoFen-7 and the ZiYuan-3 stereo images. For the experiment area of the GaoFen-7 images, the DSM of 2.5 m resolution obtained with high accuracy LiDAR is used as the reference. For the ZiYuan-3 images, the DSM in Sainte-Maxime of France based on French DSM from very high resolution aerial images is used to evaluate the DSM generated with the proposed method quantitatively. Its spatial resolution has been reduced from 0.4 to 10 m. Therefore, the surface details might be missing, but the accuracy should remain sufficiently high to serve as reference. As baseline, we choose a software solution for multi-view stereo SURE (Rothermel and Wenzel 2012). The libtsgm with the core functionality for dense matching provided by Rothermel and Wenzel (2012) is used to match the disparity map, while the generations of point clouds from disparity maps and the generation of grid DSM are the same as the proposed method.

For the GaoFen-7 and ZiYuan-3 stereo images, mSGM is first used to generate the initial DSM. Then, the initial DSM is refined by optimizing the objective function mentioned above, and the final spatial resolutions of the reconstructed DSM are 1 and 5 m, respectively. The proposed method is implemented using C++ with external dependencies: Ceres Solver (Agarwal and Mierle 2018) and OpenCV (Bradski 2000). We experiment on a standard Windows 10 computer with Intel Xeon CPU and 64 GB of memory without GPU optimization. Throughout the experiments, we use the same values for all optimization parameters in SADGE: $\mathit{\alpha }$ = 0.1, $\mathit{\beta }$ = 0.3, and $\mathit{k}$ = 10 with the basic rule that no substantial difference should exist among the terms in the objective function.

To quantitatively evaluate the proposed method, three areas with different land covers are chosen as the experiment area. The three land covers are field, urban, and hill. In addition, SURE (Rothermel and Wenzel 2012) is used as a baseline method to verify the accuracy of the proposed SADGE.

For the GaoFen-7 images, a 300 × 300 m field area is chosen for quantitative evaluation. As shown in Figure 2, the DSM reconstructed by the proposed SADGE is more similar to the reference DSM than the DSM reconstructed by SURE in the field area. Apparently, many noises are found in the DSM reconstructed by SURE, whereas the DSM generated by the proposed SADGE is considerably more flat. As for the detailed shape, the proposed SADGE can successfully recover the detailed shape such as the road in the field area. This result demonstrates that the proposed method not only can recover fine details but also reduce noises in the field area.

Figure 2. Results of the field area with the GaoFen-7 images. (a) The reference image, (b) the ground truth obtained with LiDAR, (c) the DSM generated with SURE, (d) the DSM generated with the proposed SADGE.

For the ZiYuan-3 images, two 2.5 × 2.5 km areas with urban and hill land covers are chosen for quantitative evaluation. As shown in Figure 3, the DSM reconstructed by the proposed SADGE is “sharper” than that reconstructed by SURE in the urban area. Compared with SURE, the proposed SADGE evidently reconstructs more detailed shapes, especially in the edge areas of the buildings. In addition, the DSM by the proposed SADGE is smoother than SURE in the flat area shown in the second and third columns in Figure 3.

Figure 3. Results in Sainte-Maxime of France of the urban area. From top to bottom are the nadir ZiYuan-3 images (a), the rendered reference DSM based on French DSM (b), the rendered DSM generated by the proposed method (c), and the rendered DSM generated by SURE (d). The second and third columns are for the up and down small squares in the scene, respectively.

As shown in Figure 4, the DSM reconstructed by the proposed SADGE is “clearer” than that reconstructed by SURE in the hill area. Compared with SURE, the proposed SADGE obviously reconstructs the shapes well and is more similar to the reference, especially in the second and third columns. As mentioned above, the reference DSM has been reduced from 0.4 to 10 m in the spatial resolution. The lack of details is anticipated, whereas the proposed SADGE can reconstruct more detailed shapes consistent to the nadir image.

Figure 4. Results in Sainte-Maxime of France of the hill area. From top to bottom are the nadir ZiYuan-3 images (a), the rendered reference DSM based on French DSM (b), the rendered DSM generated by the proposed method (c), and the rendered DSM generated by SURE (d). The second and third columns are for the left and right small squares in the scene, respectively.

Figure 5 further shows the accumulated distribution of DSM errors of the DSM reconstructed by SADGE and SURE. Apparently, SADGE has more points with smaller DSM error than SURE in all the experiment areas.

Figure 5. Accumulative distribution of the DSM errors of the DSM generated by SADGE and SURE in the field (a), urban (b), and hill areas (c).

To evaluate the accuracies of the DSMs reconstructed by the proposed SADGE and SURE, Table 1 calculates and lists the mean, standard deviation and root mean square of the DSM errors in the field, urban and hill experiment areas. As shown in Table 1, the proposed SADGE achieves the highest accuracy in all the evaluated metrics and has an improvement of about 56.3% in the field area, 15.7% in urban area, and 12.5% in hill area compared with SURE in terms of Root Mean Square Error (RMSE).

Table 1. Quantitative evaluation of SADGE and SURE.

	Method	Mean±Std (m)	RMSE (m)	Improvement (%)
Area 1 (field)	SADGE	0.42±0.47	0.63	56.3%
	SURE	1.30±0.60appe	1.44
Area 2 (urban)	SADGE	1.91±1.58	2.48	15.7%
	SURE	2.16±1.88	2.87
Area 3 (hill)	SADGE	2.44±1.82	3.05	12.5%
	SURE	2.75±2.06	3.43

5. Conclusions

The current study proposes a novel Shading Aware DSM GEneration (SADGE) method based on high-resolution satellite images. The DSM is reconstructed with a combination of image matching and SfS, i.e. the initial surface is reconstructed with image matching, and the final DSM is obtained through the refinement using SfS. For image matching, an improved semi-global matching method is used. As for DSM refinement, it is achieved through solving an object function considering the illumination, surface albedo and surface normal.

Two kinds of satellite images are utilized to evaluate the proposed SADGE: ZiYuan-3 and GaoFen-7 images. To better evaluate the reconstruction accuracy, three areas with different types of landcover are used: field, urban and hill. All the experiments demonstrate that this method can recover more detailed shapes and improve the accuracy of the DSM up to 56.3% in the field area, 15.7% in the urban area, and 12.5% in hill area than a software solution for multi-view stereo (SURE).

Owing to the compelling performance of combining image matching and shading in DSM reconstruction, this study recommends considering shading as an essential part of DSM generation, such as SADGE. However, the assumptions of SADGE are still not ideal, e.g. the Lambertian reflectance model and the lack of handling shadows. Thus, exploring more accurate reflectance models and better handling shadows can be interesting future work.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The Sainte-Maxime dataset used in the study is available from https://www.isprs.org/data/zy-3/Default-HongKong-StMaxime.aspx, other datasets are available from China Center For Resources Satellite Data and Application (CRESDA), http://www.cresda.com/CN/index.shtml, but restrictions apply to the availability of these data, which were used under license for this study.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [grant number 41801390] and the National Key R&D Program of China [grant number 2018YFD1100405].

Notes on contributors

Zhihua Hu

Zhihua Hu received the PhD degree in photogrammetry and remote sensing from Wuhan University in 2020. His current research interests include mesh refinement, and multi-view images 3D reconstruction.

Pengjie Tao

Pengjie Tao is currently an associate research fellow. His research interests include photogrammetry, registration of optical images and LiDAR points, and multi-view images 3D reconstruction.

Xiaoxiang Long

Xiaoxiang Long is currently a research fellow. His research interests include satellite photogrammetry, and multi-view images 3D reconstruction.

Haiyan Wang

Haiyan Wang is currently an engineer in the field of satellite mapping. His research interests include satellite photogrammetry, and multi-view images 3D reconstruction.