3D Gaussian Splatting for Large-scale Surface Reconstruction from Aerial Images

Recent advances in Neural Radiance Fields (NeRF) [7] have significantly influenced the Novel View Synthesis (NVS) domain by employing neural networks to learn and render high-quality 3D representations of continuous volumetric scenes from images taken from multiple viewpoints. NeRF achieves high-fidelity scene representation by predicting density and color. However, the huge computational demands and the time required for both training and rendering present significant challenges for real-time rendering in large-scale scenarios. To address NeRF’s computational demands, Plenoxels [29] uses 3D sparse grids to represent scene points, reducing computational complexity and storage requirements while enhancing speed compared to vanilla NeRF. However, the use of voxel grids in Plenoxels leads to a degradation of fine details. Efforts like Mip-NeRF [30] and Tri-MipRF [31] enhance rendering quality with multi-scale representations and anti-aliasing, achieving better visual quality without sacrificing efficiency.

More recently, 3DGS [12] has gained attention for its ability to achieve high-fidelity and real-time rendering. Subsequent work, such as Mip-Splatting [13], employs 3D smoothing filters to improve rendering quality. Scaffold-GS [14] leverages anchor points to regulate local 3D Gaussian distributions, allowing for real-time adjustments to both the distribution and density of Gaussians.

Image-based reconstruction has advanced considerably over the decades. Semi-global matching (SGM) [32] and patch-based methods [33] have been widely used for dense image matching. Complete surface reconstruction solutions, such as Colmap [15][16] and OpenMVS [17], utilize dense matching techniques to generate dense point clouds or triangulated meshes. With the rapid development of deep learning, learning-based MVS methods [18] [19] have been developed to predict depth maps from multi-view images, and some of these methods are integrated into comprehensive surface reconstruction frameworks [24].

More recently, several works have attempted to apply NeRF or 3DGS-based methods to surface reconstruction, mainly for small-scale or foreground objects. NeuS [10] combines the strengths of both volumetric and surface rendering by optimizing Signed Distance Functions (SDF) and a color field to reconstruct fine surface details, but it requires substantial computational resources and inference time. Neuralangelo [11] improves surface reconstruction fidelity by combining multi-resolution 3D hash grids with neural surface rendering, but it also demands significant computational power.

The emergence of 3DGS [12] has introduced new approaches for surface reconstruction. SuGaR [34] compresses 3D Gaussian spheres into approximate 2D ellipses during training and utilizes Poisson reconstruction to extract continuous mesh from 3D point clouds, which are sampled based on the Gaussian density field. However, the absence of geometric constraints leads to dispersed point clouds and numerous holes in the meshes. 2D Gaussian Splatting [23] replaces 3D Gaussian ellipsoids with 2D disks, addressing the limitations of vanilla 3DGS in surface reconstruction. Nonetheless, 2DGS still struggles to capture the intricate details of large-scale scenes.

At present, most GS-based surface reconstruction methods predominantly focus on small-scale foreground targets, and large-scale surface reconstruction remains largely underexplored.

In rendering research, the applications of the radiance field technology have extended from rendering close-range foreground objects to large-scale scenes. Block-NeRF [35] and Mega-NeRF [27] use data chunking strategies to handle large-scale scenes, while UE4-NeRF [36] integrates NeRF with Unreal Engine 4 for scalable rendering and scene editing. Switch-NeRF [37] uses a mixture of expert models for scene decomposition, enhancing large-scale scene reconstruction.

A few studies use 3DGS for large-scale scene rendering. CityGaussian [38] explores large-scale scene rendering using chunking and Level of Detail (LoD) methods to address challenges related to rendering efficiency and scalability. Similarly, VastGaussian [20] addresses lighting effects on rendering, further enhancing the visual quality of large-scale scenes.

Most of the aforementioned NeRF and 3DGS-based studies primarily focus on image rendering, with few addressing large-scale scene surface reconstruction. In contrast, our work seeks to extend 3DGS-based methods to surface reconstruction from large-scale aerial MVS images. We borrow the idea of block chunking in the large-scale rendering method VastGaussian [20] and develop a viewpoint selection and culling strategy to bridge the huge computational resource demand and limited GPU capacity. Additionally, we introduce the ray-gaussian intersection method [21][22] to determine the intersection point between the Gaussian and the ray, thereby enabling the acquisition of accurate depth and normal vector information. Furthermore, we apply multi-view geometric consistency constraints to ensure geometric consistency across different views, enhancing high-fidelity and high-precision large-scale scene reconstruction.

In the field of aerial photogrammetry, certain basic principles guide the partitioning of a scene into blocks [24] to alleviate the computational burden. In this paper, our method is inspired by VastGaussian [20] and involves three main steps. The first two steps are derived from VastGaussian. The process, as shown in Fig. 2, begins by dividing the scene based on camera positions and extending the block boundaries. Specifically, the scene, containing a total of n viewpoints, is divided into M×N blocks. First, the viewpoints are horizontally divided into M blocks, with each block containing n/M viewpoints. Then, these M blocks are further divided vertically, resulting in each block containing n/(M×N) viewpoints. After this initial partitioning, the point clouds within each region are aggregated into distinct point cloud blocks. To minimize artifacts across the scene, each region is extended by a certain proportion as in VastGaussian, ensuring that each block is adequately optimized.

However, a coarse data block partitioning method based solely on camera positions may result in insufficient optimization within each block due to suboptimal viewpoints. To address this, we develop a viewpoint selection and culling strategy, as shown in Fig. 2(c) and (d). This strategy removes erroneous viewpoints from the data block and supplements it with additional effective viewpoints. This process involves projecting all point clouds within a data chunk onto all images and calculating projection scores to determine the suitability of each viewpoint. If a point falls within the central scope (here 70%) of an image, the projection score for that image is incremented by one. For each region, the top N images with the highest scores are selected as viewpoints to optimize the data block. Finally, to ensure sufficient points for initialization and to mitigate artifacts, the sparse point cloud generated from the SfM of all new viewpoints is incorporated into the data block.

In vanilla 3DGS, the depth of each Gaussian primitive is assigned based on its distance from the screen, and accurate normal vector information is not provided. However, for surface reconstruction tasks, reliable depth information and normal vector data are essential for geometric constraints. To address these issues, we introduce the ray-gaussian intersection technique [21] [22], as shown in Fig. 3, to improve the accuracy of surface reconstruction. The intersection point $t_{max}$ between a Gaussian primitive and a ray, corresponding to the maximum Gaussian value along the ray, can be computed as follows [22]:

$o_{g}$, $r_{g}$ are $o$ (the camera center) and $r$ (the ray direction) converted into the Gaussian local coordinate system. An arbitrary point along the ray is defined as $x=o+tr$, where $t$ is the depth of the ray.

Once the depth value is computed, the Gaussian’s normal is derived as the normal of the intersection plane relative to the given ray direction. For image rendering, rather than projecting 3D Gaussians onto 2D screen space as done in the original 3DGS, we utilize the ray-gaussian intersection method to determine the intersection point between the Gaussian and the ray, allowing us to compute the contribution of a Gaussian to a given ray in 3D space.

After obtaining the depth and normal vector information, we apply 2DGS’s depth normal consistency constraints[23]. Specifically, this involves calculating the error between the normal map and the gradient values derived from the depth map, which is used as the loss function.

Here, $i$ represents the Gaussian index, $\omega_{i}$ denotes the blending weight, $n_{i}$ represents the normal vector of the Gaussian, and $N$ is the normal vector calculated from the depth map. The normal vector in the depth map at a given point is computed as Eq.(5), in which $p$ represents the pixel coordinates and $\nabla$ represents the gradient calculation.

As 3DGS-based surface construction is still in its early stages, certain beneficial empirical approaches, such as multi-view geometric consistency constraints, have yet to be fully developed. In this study, we introduce multi-view geometric consistency constraints to ensure geometric coherence across multiple views. As shown in Fig. 4, we render the depth maps for two adjacent viewpoints $V_{r}$ and $V_{n}$, denoted as $D_{r}$ and $D_{n}$. Firstly, a pixel $P$ in the reference view $V_{r}$ is projected onto the adjacent view $V_{n}$ through its depth value $Dr(P)$ and the intrinsic and extrinsic parameters, yielding the projected point $P^{\prime}$ in $V_{n}$. Subsequently, the projected point $P^{\prime}$ is reprojected onto the reference view based on its rendered depth value $Dn(P^{\prime})$, resulting in the reprojected pixel $P_{r}^{\prime}$. The distance between the coordinates of $P$ and $P_{r}^{\prime}$ is calculated as the geometric consistency constraints:

When calculating the loss, only the non-zero values are averaged, as shown in Eq.(6), where $V$ represents the valid pixels. To reduce the impact of occlusion, a distance threshold T (It is usually set to 1.) is applied to identify valid pixels, with distances $\|P-P_{r}^{\prime}\|$ exceeding T set to zero.

After parameters in each data block are optimized separately, all blocks are merged to form a coherent scene. This is achieved by removing the expanded regions of each block prior to merging.

WHU-OMVS: This dataset covers an area of Guizhou, China, with a ground resolution of 10 cm. The images are captured using a camera rig with one nadir and four oblique viewpoints, totaling 268 images, each with a resolution of 3712×5504 pixels. The flight height is 550 meters, covering an area of 850×700 $m^{2}$. Due to GPU memory limitations, we apply a 4x downsampling to the images during the training of all methods, as done in previous studies [38] [20][37][27] for rendering and reconstruction. The depth map is used as ground truth, and we evaluate the geometric accuracy of the proposed method based on the rendered depth map.

Tianjin Dataset: This dataset is captured by a camera rig with one nadir and four oblique viewpoints over Tianjin city, China, with an image size of 3840×2560 pixels and a ground resolution of 20 cm captured at a height of 200 meters. The dataset consists of 342 images and covers an area of 400×350 $m^{2}$. The image overlap is 80% along the heading direction and 60% in the side direction. The ground truth data is derived from LiDAR point cloud scans. We apply a 4x downsampling operation to all images.

Mill-19 and UrbanScene3D: We apply the proposed method to three open-source large-scale scenes: Rubble and Building from the Mill-19 dataset [27] and Residence from the UrbanScene3D dataset [28], containing 1,678, 1,940, and 2,582 images, respectively. Following previous rendering methods [38] [20][37][27], we perform a 4x downsampling operation to the input image during training. Due to the absence of ground truth, we conduct a qualitative analysis of the surface reconstruction results from these datasets.

When training the proposed method, we first perform Manhattan alignment on the target scene, aligning the y-axis to be perpendicular to the ground to facilitate the chunking process. Each block is expanded by 20%. During training, each data block is independently optimized for 50,000 iterations. The densification process begins after 500 iterations and ends at 30,000 iterations. The multi-view geometric consistency constraints and depth normal consistency constraints are introduced at 7,000 iterations. Experiments are performed on the RTX4090. Other settings remain consistent with those used in the original 3DGS method [12]. For the sparse point cloud generation, we use the SfM module in Colmap [16] [15]. For surface reconstruction, we follow the 2D Gaussian Splatting approach[23] and utilize the Truncated Signed Distance Function (TSDF) [39]. When evaluating 3DGS [12], we set the densification interval to 250 instead of the original 100 to avoid the out-of-memory problem.

We conduct surface reconstruction experiments on the WHU-OMVS [24], Tianjin, Mill-19[27] and UrbanScene3D[28] dataset. Due to the limited number of 3DGS-based surface reconstruction methods capable of handling scenes with significant depth variation, we select three for comparison: 2D Gaussian Splatting [23], 3D Gaussian Splatting [12], and Gaussian Opacity Fields [22]. Additionally, we include comparisons with the widely recognized open-source MVS software Colmap [16][15] and OpenMVS [17].

Beyond surface reconstruction, we also evaluate the rendering quality of the proposed method in Mill-19, UrbanScene3D, and WHU-OMVS. For rendering comparisons, we select four state-of-the-art methods for large-scale rendering: Mega-NeRF [27], Switch-NeRF [37], GP-NeRF [40], and CityGaussian [38].

We conduct ablation experiments mainly on the WHU-OMVS dataset.

Viewpoint Selection and Culling: As shown in Table V, disabling the viewpoint selection and culling strategy results in a significant reduction in reconstruction accuracy for the proposed method. This decline primarily stems from the insufficient number of views. As depicted in Fig. 12, the lack of proper viewpoints leads to severe errors in localized areas of the depth error map, and certain parts of buildings in the normal map appear transparent and under-optimized. Thus, viewpoint selection and culling are crucial components of the chunking strategy, ensuring better coverage and optimization.

Ray-Gaussian Interaction: As shown in Table VI, the introduction of ray-gaussian interaction significantly improves reconstruction quality, evidenced by a 19.60% improvement in $PAG_{0.6m}$. In the original 3DGS method [12], depth maps are generated using the initialized depth from Gaussians, and normal vector information is unavailable. In contrast, our approach accurately captures both depth and normal vectors, enabling the application of geometric constraints. As a result,the ray-gaussian intersection effectively mitigates the challenges associated with surface reconstruction, particularly those caused by the irregular distribution of Gaussian primitives.

Multi-View Geometric Consistency Constraints: The original 3DGS [12] applies loss within a single view, which can lead to overfitting and fails to ensure consistency across multiple viewpoints. To address this issue, we introduce multi-view geometric consistency constraints. Table VI demonstrates that the multi-view geometric consistency constraints effectively enhance reconstruction accuracy by ensuring coherence across multiple views.