City-Mesh3R: Simulation-Ready City-Scale 3D Mesh Reconstruction
from Multi-View Images
May 28, 2026
City-scale 3D surface reconstruction from multiview images for downstream 3D simulation, poses highly challenging problems due to the scale and complexity of urban scenes. Existing city-scale 3D reconstruction methods based on NeRF, Gaussian Splatting etc. often fail to recover 3D meshes ready for simulation due to incomplete/missing geometry and irregular, noisy surfaces. Scaling existing small-scale 3D reconstruction methods to arbitrarily large urban scenes is highly infeasible due to their computational complexity. We present City-Mesh3R, a scalable framework for reconstructing watertight surface meshes directly from large unordered image collections. Unlike recent methods which use global sparse SfM point-cloud initialization followed by a distributed 3D dense reconstruction of large-scale scenes, our method follows an end-to-end images-to-mesh 3D reconstruction approach using a divide-and-conquer strategy. The sparse city map is reconstructed via topological image clustering, cluster-wise independent sparse SfM and map merging, without need for exhaustive image feature matching. Then this map is partitioned spatially to perform geometry-aware camera selection, followed by dense surface reconstruction and surface refinement using curvature-aware adaptive vertex density remeshing. These partition meshes are then stitched together to produce the global mesh of the city. The proposed end-to-end framework is evaluated on city-scale reconstruction datasets. As demonstrated by our qualitative and quantitative results, our proposed method yields high-fidelity watertight 3D meshes with regular geometry, capturing fine surface details, and is suitable for scaling to arbitrarily large scenes owing to the end-to-end processing in a distributed setting.
3D reconstruction from multiview images of city-scale scenes City10892091? is a challenging research problem which directly impacts applications such as urban digital twin [1], 15-minute city planning[2], infrastructure monitoring and disaster management. These applications often require simulation of high-fidelity city-scale 3D meshes that can capture the structural details of buildings, public infrastructure etc. in real-world urban scenes. However, scaling to arbitrarily large city scenes is challenging owing to a trade-off between geometric fidelity and computational complexity.
Recent research use NeRF [3]–[5] or Gaussian Splatting representations [6]–[8] for city-scale 3D reconstruction. 3D Surface reconstruction using Gaussian Splat representations can result in ill-defined surface boundaries, making it challenging to extract consistent normals, preserve sharp features, or perform structural edits such as adding or reshaping buildings. For downstream tasks like for e.g. urban planning, infrastructure modification, and hydrodynamic flood simulation; explicit mesh-based representations remain the practical and simulation-ready solution. Extracting 3D mesh surfaces from NeRF or Gaussian Splats often results in noisy, incomplete surfaces (holes or tears) that need manual postprocessing for simulation. Moreover, NeRF-based methods are computationally slow, while Gaussian-Splat–based approaches incur additional overhead for mesh extraction. Inaccuracy in 3D geometry directly impacts applications such as urban infrastructure monitoring, disaster management, etc. Moreover, methods based on gaussian-splatting rely on a good SfM point cloud initialization [9] typically computed using exhaustive image feature matching with huge computational complexity. Also, existing methods for 3D mesh reconstruction cannot scale to arbitrarily large city scenes captured by thousands or millions of images.
In this work, we present a novel pipeline for scalable 3D mesh reconstruction that produces simulation-ready watertight 3D meshes of city-scale scenes. Following large-scale 3D reconstruction methods, we apply a distributed approach to partition the 3D urban scene, reconstruct individual partitions and merge the individual maps into a global 3D city map. Our proposed pipeline applies the distributed reconstruction in a sparse-to-dense manner. It consists of a distributed SfM stage to generate a city-scale sparse point cloud followed by a distributed approach for 3D mesh reconstruction and surface optimization. Unlike the recent city-scale mesh reconstruction techniques that use a given initial sparse global point cloud from a non-scalable standard SfM pipeline [9], our proposed end-to-end pipeline reconstructs both sparse point cloud and 3D mesh in a distributed manner, making it easy to scale to arbitrarily large scenes. Unlike city-scale reconstruction methods using NeRF/GS which perform global 3D mesh surface extraction and are not scalable, our proposed method follows a distributed approach for 3D mesh reconstruction and optimization, followed by mesh stitching.
Instead of relying on exhaustive pairwise feature matching in a standard SfM pipeline [9], [10], we apply community-based image clustering to obtain coarse topological clusters of the scene for distributed sparse SfM. The sparse point clouds reconstructed for overlapping clusters are merged to reconstruct a
globally aligned sparse map of the entire city. We then apply geometry-aware partitioning of the globally aligned sparse map, leveraging camera poses and co-visibility statistics to produce data partitions suitable for downstream dense reconstruction. This
hierarchical partitioning strategy yields well-formed map partitions, mitigates redundancy and ambiguity caused by irregular topological cluster boundaries, and significantly improves the robustness of subsequent dense map merging. We adapt a fast mesh
optimization and remeshing procedure based on differentiable rendering for each partition. Our formulation introduces adaptive vertex density remeshing based on local scene complexity instead of uniform distribution, to produce watertight meshes with fine
structural details. This adaptive formulation facilitates accurate geometric recovery, better optimization stability and memory efficiency, making it particularly well suited for large-scale urban environments.
Our main contributions are as follows:
We present a distributed hierarchical sparse-to-dense pipeline for reconstructing city-scale, simulation-ready, high-fidelity watertight surface meshes from large unordered image collections, enabling scalability to arbitrarily large scenes.
We propose a two-stage divide-and-conquer strategy for large-scale reconstruction that decouples sparse and dense partitioning to match their different data requirements: topological image clustering for distributed SfM, followed by spatial partitioning with geometry-aware camera selection for dense reconstruction.
We introduce a curvature-aware adaptive remeshing strategy for large urban scenes that allocates more vertices to geometrically complex regions and fewer to flatter areas, improving reconstruction quality, stability, and memory efficiency.
We establish the efficacy of our proposed pipeline in reduced computational complexity and improved mesh fidelity on urban-scale 3D benchmark datasets.
Large-scale 3D reconstruction has been explored through diverse representations, including point clouds, volumetric NeRFs, Gaussian splats, and meshes. Early systems such as [11] demonstrated city-scale sparse point cloud reconstruction from 150,000 images via distributed matching and parallel SfM, while large-scale Structure-from-Motion [9], [12] and MVS methods [13] methods do not focus on reconstruction mesh surfaces. [14] introduced the first million-image global SfM using a divide-and-conquer camera partitioning framework.
Recent research on large-scale 3D reconstruction increasingly adopts neural scene representations. NeRF-based approaches [3]–[5] perform block-wise training for large scenes, and SDF-based models [15] leverage dense fields for geometry extraction. Recent methods apply Gaussian-splatting [6]–[8], [16]–[19] to enable real-time rendering of large-scale urban scenes, with CityGaussian [7] and CityGaussianV2 [8] scaling 3DGS training to multi-square-kilometer scenes; CityGS-X [20] accelerates this using multitask rendering; and DOGS [17] applies ADMM-based block optimization. Horizon-GS [21] fuses aerial and street views using 3DGS for rendering and 2DGS for geometry.
Surface extraction from Gaussian or neural fields has been explored through Poisson surface reconstruction in SuGaR [22], SDF-based depth fusion in 2DGS [23], opacity-field iso-surfacing in GOF [24], planar-disk splatting in PGSR [25], and hybrid SDF–GS methods such as 3DGSR [26] and SurfaceSplat [27]. Additional formulations include sparse-voxel rasterization (GeoSVR [28], [29]) and differentiable mesh extraction in MILo [30]. MeshSplatting [31] offers an efficient representation but doesn’t guarantee watertight mesh or manifold-ness, while RadianceMeshes [32] produces meshes at the cost of significantly longer reconstruction time. However, both methods remain unsuitable for distributed city-scale reconstruction due to their high per-block computational overhead. Large-scale mesh methods [18], [20], [21] typically rely on global TSDF fusion, which is computationally expensive, fails to capture thin structures, and often produces irregular surfaces with holes [8].
Moreover, Gaussian-based surface extraction often suffers from ill-defined boundaries, noisy geometry, and difficulty preserving sharp edges, limiting its applicability to downstream tasks such as urban planning, infrastructure modeling, or flood simulation. Similarly, mesh extraction from NeRF-based methods often produces irregular or incomplete surfaces and incurs substantial computation time, making these approaches impractical for large-scale reconstruction. In contrast, our method is explicitly designed for practical scalability and produces simulation-ready mesh outputs, making it directly suitable for real-world city-scale digital-twin applications.
Given a sequence of unordered images \(\mathcal{I} = \{I_1, I_2, \dots, I_N\}\) capturing a large-scale scene, our proposed method (shown in Fig. 1) reconstructs the global mesh representing the scene. Our proposed pipeline first reconstructs a sparse 3D point cloud of the entire scene using a distributed SfM approach (Images to Sparse SFM stage). Then the sparse 3D point cloud is converted to a watertight mesh via dense surface reconstruction followed by mesh surface refinement in a partitioned approach (Sparse SFM to Surface Reconstruction stage). The next sections describe the different stages of our proposed pipeline.
Global Feature Extraction and Similarity Graph Construction: Traditional large-scale SfM systems (e.g., [33]) rely on local feature matching (e.g., SIFT) to construct the image view graph, requiring exhaustive pairwise descriptor matching, which becomes computationally prohibitive at scale. To address this, modern pipelines employ global image descriptors such as AP-GeM [34] and distilled NetVLAD [35] used in [36] and [37] for efficient image retrieval. These approaches adopt a two-stage strategy: global descriptors first identify candidate image pairs, followed by local feature matching, significantly reducing computational cost and improving scalability. We construct the view graph using global image descriptors. Specifically, we employ the DINOv2 [38] vision transformer model to extract global features from images due to its strong generalization capability across diverse datasets.
For each image \(I_i\), we compute a global feature vector \(\mathbf{f} = \Phi(I) \in \mathbb{R}^d\), where \(\Phi(\cdot)\) denotes the DINOv2 feature extractor. Using these feature vectors, we construct an image similarity graph \(G = (V, E)\), where each vertex \(v_i \in V\) corresponds to an image \(I_i\). An edge \((i,j) \in E\) is created if the similarity score satisfies \(d_S(\mathbf{f}_i,\mathbf{f}_j) \geq \tau\) ,where \(d_S(\cdot,\cdot)\) denotes the cosine similarity value and \(\tau\) is a predefined similarity threshold.
Image clustering with overlaps: Commonly used clustering approaches such as Normalized Cuts (Ncut) in [33] or Agglomerative Clustering [39] typically produce disjoint clusters and do not guarantee overlaps between them. However, for our distributed sparse SfM we need some overlap between clusters for a reliable merging into a global 3D map. To achieve this, we apply the Speaker–Listener Label Propagation Algorithm (SLPA) [40], which naturally detects overlapping communities in graphs, thereby producing clusters that share common images and facilitate subsequent merging of individually reconstructed point clouds for each cluster.
In SLPA each nodes of the similarity graph maintains a memory of observed labels over time. At iteration \(t\), a node \(i\) (speaker) selects a label from its memory according to a
speaking rule (typically the most frequent label) and sends it to a randomly chosen neighbor \(j\) (listener). Formally, let \(M_i^{(t)}\) denote the label memory of node \(i\). The speaker samples a label \(\ell \sim P(\ell \mid M_i^{(t)})\) and transmits it to node \(j\). The listener then updates its memory as \(M_j^{(t+1)} = M_j^{(t)} \cup \{\ell\}\). After convergence, the probability of label \(\ell\) at node \(i\) is \(P_i(\ell) =
\frac{\text{count}_i(\ell)}{|M_i|}\). A threshold is applied to remove weak memberships on the above calculated probability. This produces a set of overlapping communities \(\mathcal{C} = \{C_1, C_2, \dots, C_m\}\)
which naturally capture shared viewpoints across image groups. To ensure connectivity between clusters, we enforce a minimum overlap constraint.
Cluster Graph Construction, Cluster-wise sparse SfM, and Hierarchical Cluster Map Merging: We build a cluster graph \(G_c = (V_c, E_c)\) where nodes denote clusters and edge weights \(w_{ij} = |C_i \cap C_j|\) capture shared images. A Minimum Spanning Tree (MinST) (via Kruskal) defines the merging order. The base cluster is selected by solving a Minimum Height Tree (MHT) on the MinST through iterative leaf
removal, choosing the remaining node (or the larger one if two remain). Clusters are then merged hierarchically from leaves to the base.
Within each cluster, feature matching is performed using MASt3R [41] restricted to edges in the original similarity graph, and sparse reconstruction is obtained using COLMAP [9]. Inter-cluster alignment is achieved by estimating a \(\text{Sim}(3)\) transformation (rotation, translation, scale) via RANSAC on shared camera poses. The source reconstruction is transformed into the target frame, followed by track merging based on shared observations and removal of high reprojection error outliers to ensure consistency.
Area partitioning of large-scale sparse SfM: To scale reconstruction to large scenes, we first spatially partition the sparse COLMAP point cloud and then select, for each partition, a compact set of informative cameras. We parameterize the scene using a dominant support plane, assuming only a single dominant supporting surface rather than Manhattan-world structure. Projecting the sparse 3D points onto this plane yields a 2D footprint that is better suited to area-based decomposition than the original SfM frame. We then rotate the in-plane coordinates to compact the footprint, divide the resulting domain into a regular grid, and slightly enlarge each cell to define the effective partition support. This introduces moderate overlap between neighboring partitions and improves robustness near partition boundaries. Full details of support-plane fitting, planar parameterization, orientation refinement, and inflated grid construction are provided in the supplementary (Sec. 6.1).
For partition \((r,c)\) with sparse support \(\mathcal{P}_{r,c}\), we define the candidate camera pool as \(\mathcal{I}^{\mathrm{cand}}_{r,c}=\bigcup_{p\in\mathcal{P}_{r,c}}\mathcal{I}(p)\), namely, all cameras that observe at least one point in the partition. Rather than ranking cameras individually, we rank co-visible camera pairs, since local reconstruction quality depends on both visibility and geometric complementarity. We retain only camera pairs with sufficient local co-visibility within the partition and assign each admissible pair a prior that favors moderate relative viewpoint change and strong local overlap. For each partition point \(p\), each admissible observing pair is further scored using two point-dependent cues: triangulation quality and image-plane centrality. Thus, preferred pairs are well supported within the partition, provide a useful baseline for \(p\), and observe \(p\) away from image boundaries. To avoid over-counting many near-redundant pairs around the same point, we retain only the top-\(K\) pairs for each point and normalize their scores into a unit vote. Each selected pair then splits its vote equally between its two cameras, and votes are accumulated over all partition points to obtain a camera score \(S_i\) for each candidate view. Finally, the top \(M\) cameras form the partition-specific subset \(\mathcal{I}^{\mathrm{top}}_{r,c}\).
Compared with simple visibility counting, this ranking favors cameras that are not only present in a partition, but also geometrically informative for local reconstruction. Detailed definitions, scoring functions, and normalization steps are provided in
supplementary material Sec. 6.2.
Dense Reconstruction and Surface Initialization:-
For each partition, we retain the top-\(M\) cameras from the ranking stage and fix their SfM intrinsics and poses. We then densify the partition using MASt3R depth predictions and pixel correspondences [41], reusing cached outputs from the original similarity graph whenever possible. As a result, only pairs newly introduced by the final partitioning, typically near former inter-cluster boundaries, require additional inference. To address residual inter-view scale mismatch and local geometric noise in the per-view dense predictions, we perform a two-stage geometry alignment procedure inspired by MASt3R-SfM [42], but tailored to fixed SfM cameras. We first estimate one global scale per view to align matched dense 3D points across overlapping images, then refine the dense geometry by enforcing multi-view reprojection consistency under the frozen cameras. This removes the dominant scale inconsistency, suppresses local depth errors, and improves cross-view agreement. The full formulation is deferred to the supplementary (7).
Finally, we fuse the aligned depth maps with TSDF integration, extract oriented points from the fused volume, and reconstruct a watertight mesh using Screened Poisson Surface Reconstruction. The resulting mesh provides the initialization for the
subsequent differentiable-rendering-based refinement stage. Overall, the module transforms sparse SfM cameras and largely cached MASt3R predictions into a dense, globally aligned, and denoised surface initialization without re-estimating camera
poses.
Mesh Refinement: Starting from the watertight initial mesh \(M_0=(V_0,F_0)\), we refine the surface by directly optimizing
the mesh using losses based on differentiable rendering. For view \(j\), let \(S_j(u)\in[0,1]\) be the target silhouette and \(\hat{S}_j(u;M^k)\) the
rendered silhouette. The silhouette loss is \[\mathcal{L}_{\mathrm{sil}}^{(j)}(M^k)=
\frac{1}{|\Omega_{\mathrm{img}}|}
\sum_{u\in\Omega_{\mathrm{img}}}
\big(\hat{S}_j(u;M^k)-S_j(u)\big)^2.\] Let \(\mathcal{N}_j(u)\in\mathbb{S}^2\) be the predicted unit normal map in camera coordinates from a pretrained off-the-shelf monocular normal estimator [43] and \(\hat{\mathcal{N}}_j(u;M^k)\) the rendered unit normal map. On the
foreground support \(\Omega_j=\{u\mid S_j(u)=1\}\), the normal-map loss is \[\mathcal{L}_n^{(j)}(M^k)=
\frac{1}{|\Omega_j|}
\sum_{u\in\Omega_j}
\left(
1-\hat{\mathcal{N}}_j(u;M^k)^\top \mathcal{N}_j(u)
\right).\]
At iteration \(k\), the current mesh is \(M^k=(V^k,F^k)\). Given calibrated cameras \(\{K_j,T_{wc,j}\}_{j=1}^{N}\), we optimize \[\Phi(M^k)= \sum_{j=1}^{N} \left( \lambda_n\,\mathcal{L}_n^{(j)}(M^k)+ \lambda_s\,\mathcal{L}_{\mathrm{sil}}^{(j)}(M^k) \right)+ \mathcal{R}(M^k),\] where \(\mathcal{R}\) is a lightweight mesh regularizer, e.g.Laplacian smoothing.
Gradients of \(\Phi\) are backpropagated through the renderer to vertex positions, which are updated using the Isotropic Adam optimizer [44]. The optimizer also provides a relative vertex velocity \(\nu^k(i)\in[0,1]\) [44], which measures the normalized update magnitude at vertex \(i\) and serves as a local indicator of stabilization.
Our mesh optimization method adapts Continuous-Remeshing [44], originally designed for object-level remeshing, to urban scene mesh optimization. A fixed vertex budget tends to under-refine detailed regions, while allowing the mesh to grow freely during alternating optimization and remeshing can rapidly increase the number of vertices and faces, leading to unstable memory usage, occasional out-of-memory failures, and slower iterations. More importantly, uniform refinement spends vertices where they are not needed, whereas concentrating them in detail-rich regions improves approximation quality under the same budget and can also improve optimization stability by keeping the number of optimizable parameters smaller. To address this, instead of the edge-length consistency constraint in [44], we use a curvature-guided terminal target and a speed-aware slack schedule, which are better suited to large, spatially heterogeneous scenes.
Our key idea is to decouple where fine resolution is needed from when that refinement is released. To determine where, we estimate a screen-space normal-variation field from the predicted normal maps. For each view \(j\), let \(N^{(j)}(u,v)\in\mathbb{R}^{3\times 2}\) denote the image-space Jacobian of the unit normal field. We define the local normal-rotation rate as \(s^{(j)}(u,v)=\|N^{(j)}(u,v)\|_2\), and convert a normal-rotation tolerance \(\theta_0\) into a projected target edge length \[p_{\mathrm{tgt}}^{(j)}(u,v)= \mathrm{clip}\!\left(\frac{\theta_0}{s^{(j)}(u,v)},\,p_{\min},\,p_{\max}\right). \label{eq:main95pixel95target}\tag{1}\] Thus, regions with rapid normal variation receive smaller projected target edges. For a vertex \(i\) visible in view \(j\), with projection \(\pi^{(j)}(\mathbf{v}_i)\), depth \(z_i^{(j)}\), and focal length \(f^{(j)}\), this yields a world-space target \(L_{ij}=\frac{z_i^{(j)}}{f^{(j)}}\,p_{\mathrm{tgt}}^{(j)}(\pi^{(j)}(\mathbf{v}_i))\). Pooling across visible views \(\mathcal{V}_i\) gives the per-vertex curvature-guided baseline \(L_{\mathrm{ref\text{-}curv}}(i)=\operatorname*{median}_{j\in\mathcal{V}_i} L_{ij}\).
To determine when refinement is allowed, we introduce a nonnegative slack \(\zeta^k(i)\) on top of this baseline and update it using the relative vertex velocity: \[\zeta^{k}(i)\leftarrow \max\!\Big(0,\zeta^{k-1}(i)+g\big(\nu^k(i)-\nu_{\mathrm{ref}}\big)\Big). \label{eq:main95slack95update}\tag{2}\] The final per-vertex reference edge length is \[L_{\mathrm{ref}}^{k}(i)= \mathrm{clip}\!\Big(L_{\mathrm{ref\text{-}curv}}^{k}(i)+\zeta^k(i),\,L_{\min},\,L_{\max}\Big). \label{eq:main95ref95vertex}\tag{3}\] Edge targets are derived from the vertex field during remeshing; details, initialization, and implementation choices are given in [subsec:supp_curvature_field] [subsec:supp_slack] [subsec:supp_bounds] [subsec:supp_remeshing]. Intuitively, curvature specifies the terminal resolution, while velocity delays refinement until local motion has sufficiently stabilized.
After each vertex update, we invoke continuous remeshing under the current reference field \(L_{\mathrm{ref}}^k\): short edges are collapsed, long edges are split, and flips improve local triangle quality and valence regularity, using the same validity checks as in [44]; the exact rules are given in 8.5. The global bounds \(L_{\min}\) and \(L_{\max}\) are estimated automatically from camera intrinsics, robust depth statistics, and current mesh scale as described in 8.4.
The refinement therefore alternates render \(\rightarrow\) compare \(\rightarrow\) update \(\rightarrow\) remesh, and the final mesh is selected as the best snapshot \(M^\star=\arg\min_k \Phi(M^k)\).
Our adaptation of [44] using a curvature-guided terminal target and a speed-aware slack schedule concentrates triangles in
high-curvature regions while keeping flatter regions coarse, making remeshing more stable and memory-efficient for large heterogeneous urban scenes.
Mesh Stitching: To integrate two overlapping surface meshes \(M_1\) and \(M_2\), we use a decomposition-reconstruction strategy that preserves high-frequency
geometric detail while ensuring a watertight transition in the overlap region. For the meshes \(M_1\) and \(M_2\) with bounding volumes \(B_1\) and \(B_2\) obtained from the preceding partitioning stage (Sec. [sec:area95decomp95ranking]), the decomposition step
extracts the non-overlapping exterior regions by clipping each mesh against the other’s volume: \[M_{1 \setminus 2} = M_1 \setminus B_2, \quad M_{2 \setminus 1} = M_2 \setminus B_1\] where \(\setminus\) denotes clipping against the opposing bounding volume. To resolve the discontinuity at the intersection, we extract the overlap region \(\mathcal{O} = B_1 \cap B_2\). Let \(P_{int}\) be the set of vertices residing in the overlap region. We reconstruct the bridging surface \(M_{seam}\) by applying a Delaunay triangulation \(D\) over
the projected manifold, i.e., \(M_{seam} = D(P_{int})\).
The final merged mesh \(M_{final}\) is obtained by combining the preserved exterior surfaces with the reconstructed seam patch, followed by a stitching operation along their boundary vertices: \[M_{final} = \mathrm{Stitch}(M_{1 \setminus 2},\, M_{2 \setminus 1},\, M_{seam})\] where \(\mathrm{Stitch}(\cdot)\) denotes a boundary-aware topological refinement operator that aligns and welds corresponding boundary vertices within a distance threshold \(\epsilon\), removes duplicate or degenerate faces, and enforces manifold connectivity across the three mesh components. Please refer to 6 in supplementary, for mesh stitching visualization.
Our experiments cover five representative scenes drawn from two datasets: CUHK-LOWER, CUHK-UPPER, LFLS, SZIIT from GauU-Scene [45] and Residence from UrbanScene3D [46].
Our method was implemented on a (2 x 48 GB) NVIDIA A6000 Ada GPU system and all methods including ours were evaluated on the same system.
4pt
| Scene | Method | Precision \(\uparrow\) | Recall \(\uparrow\) | F1 \(\uparrow\) | Time (min) \(\downarrow\) |
|---|---|---|---|---|---|
| CUHK-LOWER | CityGS-v2 | 0.1312 | 0.0820 | 0.1009 | 340.90 |
| CityGS-X | 0.1327 | 0.0758 | 0.0965 | 75.00 | |
| Ours | 0.1420 | 0.0910 | 0.1110 | 95.00 | |
| CUHK-UPPER | CityGS-v2 | 0.2203 | 0.0953 | 0.1330 | 371.05 |
| CityGS-X | 0.2950 | 0.1510 | 0.2000 | 74.00 | |
| Ours | 0.2680 | 0.1120 | 0.1580 | 113.00 | |
| LFLS | CityGS-v2 | 0.0901 | 0.0666 | 0.0791 | 296.90 |
| CityGS-X | 0.0978 | 0.0599 | 0.0743 | 87.00 | |
| Ours | 0.0951 | 0.0777 | 0.0855 | 83.00 | |
| SZIIT | CityGS-v2 | 0.1513 | 0.0488 | 0.0738 | 299.07 |
| CityGS-X | 0.1849 | 0.0538 | 0.0833 | 94.00 | |
| Ours | 0.1925 | 0.0647 | 0.0968 | 92.00 |
2pt
| Method | Time (hrs)\(\downarrow\) | 2D Rep. Err. (px) \(\downarrow\) |
|---|---|---|
| MASt3R-COLMAP (all images) | 53.88 | 1.17 |
| MASt3R-GLOMAP (all images) | 19.45 | 1.856 |
| Ours (Clustering + MASt3R-COLMAP) | 2.74 | 1.38 |
As shown in Table 1, our full method achieves a strong trade-off between reconstruction quality and runtime across diverse large-scale scenes. It attains the best overall performance on CUHK-LOWER and SZIIT, where it achieves the highest Precision, Recall, and F1-score, and also performs best on LFLS in terms of Recall, F1-score, and runtime. Relative to CityGS-v2 [8], our method consistently improves both accuracy and efficiency, delivering higher geometric quality with substantially lower runtime on all evaluated scenes. Relative to CityGS-X [20], our method produces better reconstruction quality on most scenes while retaining comparable execution time. Although CityGS-X performs best on CUHK-UPPER, our approach remains competitive and continues to offer a significantly better runtime than CityGS-v2. Overall, these results demonstrate that our pipeline is an effective and scalable solution for city-scale surface reconstruction. In Table 2, we also demonstrate the effectiveness of our proposed Image to Sparse SfM stage. We achieve a significantly lower reprojection error compared to MASt3R-GLOMAP, while being substantially more efficient than MASt3R-COLMAP, which requires hours despite achieving a slightly lower error. This highlights that our approach provides a favorable trade-off between accuracy and computational cost, bridging the gap between highly accurate but slow incremental reconstruction and fast but less accurate global methods.
Qualitative comparisons in Fig. 2 further show that our method produces a watertight mesh surface with lower artifacts, smooth surface normals, and detailed geometry when compared to the latest methods of large-scale 3D reconstruction. We further evaluate multiple small-scale surface reconstruction methods [30]–[32] with our method. As shown in 3 (C), qualitative results demonstrate that our approach outperforms existing small-scale methods.
Our Proposed Image Clustering Method vs Existing Methods Ablation of different clustering methods in Table 3 shows that our
clustering method (SLPA) consistently outperforms both N-cut [33] and agglomerative clustering
[39] across different scenes. Specifically, our method achieves the lowest 2D reprojection error and significantly reduces the percentage of image
rejections during SfM compared to competing approaches. These results demonstrate that our clustering strategy produces more coherent and reliable image groupings, which directly contribute to improved reconstruction quality and stability.
Poisson Surface Reconstruction vs.Our Mesh Refinement. To evaluate the importance of the mesh refinement, we compare the mesh obtained directly from Screened Poisson reconstruction (before mesh refinement) with the output of our full
pipeline. As reported in Table 4, our proposed mesh refinement consistently improves Precision, Recall, and F1 on both scenes. 3 (A) demonstrates the clear improvement in surface details captured by the proposed mesh refinement stage when compared to the output from Poisson mesh reconstruction.
2pt
| Method | 2D Reproj. Err. \(\downarrow\) | % Img. Rej. \(\downarrow\) |
|---|---|---|
| Our Clustering (SLPA) | 1.38 | 4.373 |
| N-cut Clustering | 1.441 | 9.78 |
| Agglomerative Clustering | 1.563 | 10.27 |
2pt
| Scene | Method | P \(\uparrow\) | R \(\uparrow\) | F1 \(\uparrow\) |
|---|---|---|---|---|
| CUHK-LOWER | Ours (Poisson Only) | 0.131 | 0.049 | 0.071 |
| Ours (Poisson + Opt/Remesh) | 0.142 | 0.091 | 0.111 | |
| SZIIT | Ours (Poisson Only) | 0.1321 | 0.0163 | 0.029 |
| Ours (Poisson + Opt/Remesh) | 0.1925 | 0.0647 | 0.0968 |
3pt
| Scene-Name | Method | Precision \(\uparrow\) | Recall \(\uparrow\) | F1 \(\uparrow\) |
|---|---|---|---|---|
| CUHK-LOWER | Ours | 0.142 | 0.091 | 0.111 |
| ConRemesh | 0.1350 | 0.0632 | 0.0861 | |
| SZIIT | Ours | 0.1925 | 0.0647 | 0.0968 |
| ConRemesh | 0.1858 | 0.0317 | 0.0542 |
Continuous Remeshing vs.Our Method. We further compare our mesh refinement algorithm with the original formulation of Continuous Remeshing [44] under a fixed maximum vertex-budget cap, so that both methods operate under the same mesh-complexity constraint. 5 and 3 (B) shows that our method consistently achieves better results across both scenes. This is because near-uniform remeshing spends vertices unnecessarily in flat regions, leaving less budget for high-curvature areas and fine structures. In contrast, our method adaptively redistributes the same vertex budget using curvature-guided edge-length targets and a speed-aware refinement schedule, enabling more accurate reconstruction where detail is most needed.
We presented City-Mesh3R, an end-to-end scalable pipeline reconstructs simulation-ready 3D meshes from unordered images of large urban scenes. City-Mesh3R combines a distributed hierarchical sparse-to-dense design, a two-stage partitioning strategy for sparse SfM and dense reconstruction, and a curvature-aware adaptive remeshing module for high-fidelity surface recovery while maintaining topological consistency across independently processed clusters. Experimental results on urban-scale 3D datasets demonstrate a significant performance gain achieved by our method both in processing time, and improved 3D surface quality, important for applications in urban planning, infrastructure monitoring, etc.
This subsection provides the implementation details omitted from the main paper (Sec. [sec:area95decomp95ranking]) for constructing the spatial partitions used before camera selection.
Let \(\mathcal{P}=\{P_p\}\) denote the sparse 3D points reconstructed by COLMAP. We estimate a dominant support plane from \(\mathcal{P}\) via RANSAC plane fitting, yielding parameters \((\mathbf{n},d)\) such that \[\mathbf{n}^{\top}\mathbf{x}+d=0,\] where \(\mathbf{n}\in\mathbb{R}^3\) is a unit normal and \(d\in\mathbb{R}\) is the offset. This introduces only a weak geometric prior: unlike Manhattan-world alignment, we assume only one dominant support plane.
From the plane inliers, we construct a local planar frame \((\mathbf{o},\mathbf{u},\mathbf{v})\), where \(\mathbf{o}\) is a point on the plane and \(\mathbf{u},\mathbf{v}\) are orthonormal in-plane basis vectors. Each sparse point \(P_p\) is mapped to planar coordinates by \[\mathbf{r}_p=P_p-\mathbf{o}, \qquad u_p=\mathbf{r}_p^\top\mathbf{u}, \qquad v_p=\mathbf{r}_p^\top\mathbf{v}.\]
After fixing the plane normal, the in-plane basis remains ambiguous up to a rotation around \(\mathbf{n}\). To reduce this arbitrariness, we choose the in-plane rotation that minimizes the area of the axis-aligned bounding box of the projected sparse points. Let \((u_p^{(\phi)},v_p^{(\phi)})\) denote the planar coordinates after rotation by angle \(\phi\). We choose \[\phi^\star = \arg\min_{\phi} \left( \max_p u_p^{(\phi)}-\min_p u_p^{(\phi)} \right) \left( \max_p v_p^{(\phi)}-\min_p v_p^{(\phi)} \right).\]
Using the rotated planar coordinates, we define a regular \(R\times C\) grid over the support-plane extent. Let \[[u_{\min},u_{\max}] \times [v_{\min},v_{\max}]\] be the planar bounding box. The base grid boundaries are \[u_c = u_{\min} + \frac{c}{C}(u_{\max}-u_{\min}), \qquad v_r = v_{\min} + \frac{r}{R}(v_{\max}-v_{\min}),\] for \(c=0,\dots,C\) and \(r=0,\dots,R\).
For each grid cell \((r,c)\) with base rectangle \([u_c,u_{c+1}] \times [v_r,v_{r+1}]\), we enlarge its support around the cell center \[u^{(m)}_{r,c}=\frac{u_c+u_{c+1}}{2}, \qquad v^{(m)}_{r,c}=\frac{v_r+v_{r+1}}{2}.\] Given inflation factors \(\alpha_u\) and \(\alpha_v\), the enlarged side lengths are \[\Delta u'=(1+\alpha_u)(u_{c+1}-u_c), \qquad \Delta v'=(1+\alpha_v)(v_{r+1}-v_r).\] The inflated window becomes \[\hat{u}_0=u^{(m)}_{r,c}-\frac{\Delta u'}{2}, \qquad \hat{u}_1=u^{(m)}_{r,c}+\frac{\Delta u'}{2},\] \[\hat{v}_0=v^{(m)}_{r,c}-\frac{\Delta v'}{2}, \qquad \hat{v}_1=v^{(m)}_{r,c}+\frac{\Delta v'}{2},\] followed by clipping to the global planar extent. The point set assigned to partition \((r,c)\) is \[\mathcal{P}_{r,c} = \left\{ P_p \in \mathcal{P} \;\middle|\; \hat{u}_0 \le u_p \le \hat{u}_1,\; \hat{v}_0 \le v_p \le \hat{v}_1 \right\}.\]
These partition point sets are used by the camera-ranking procedure described next.
This subsection gives the full formulation of the partition-wise camera ranking summarized in the main paper (Sec. [sec:area95decomp95ranking]).
For each non-empty partition \((r,c)\), let \(\mathcal{P}_{r,c}\) denote the partition point set. We first form the candidate camera set as the union of all cameras that observe at least one partition point: \[\mathcal{I}^{\mathrm{cand}}_{r,c} = \bigcup_{p\in\mathcal{P}_{r,c}} \mathcal{I}(p).\] If \(\mathcal{P}_{r,c}=\emptyset\), or if no sufficiently co-visible camera pairs exist within the partition, we return an empty ranked set.
For each unordered camera pair \((i,j)\), let \(m_{ij}\) be the number of points in \(\mathcal{P}_{r,c}\) observed by both cameras, and let \(N_{r,c}=|\mathcal{P}_{r,c}|\). We retain only admissible pairs whose co-visibility ratio exceeds a threshold: \[\rho_{ij} = \frac{m_{ij}}{N_{r,c}}, \qquad \mathcal{A} = \{(i,j)\mid \rho_{ij}\ge \tau_{\mathrm{cov}}\}.\]
For each admissible pair \((i,j)\in\mathcal{A}\), we define a prior that captures viewpoint complementarity and local support. Let \(\theta_{ij}\) denote the relative rotation angle between cameras \(i\) and \(j\). We favor moderate relative orientation using \[P_{\mathrm{rot}}(i,j) = \exp\!\left( -\frac{(\theta_{ij}-\mu_{\mathrm{rot}})^2}{2\sigma_{\mathrm{rot}}^2} \right).\] We also map the pairwise co-visibility count to a normalized overlap term \[P_{\mathrm{overlap}}(i,j) = \frac{\log(1+m_{ij})}{\log(1+m_{\max})}, \qquad m_{\max}=\max_{(i,j)\in\mathcal{A}} m_{ij}.\] The resulting pair prior is \[\operatorname{prior}_{ij} = P_{\mathrm{rot}}(i,j)^{w_{\mathrm{rot}}} \cdot P_{\mathrm{overlap}}(i,j)^{w_{\mathrm{overlap}}}.\]
For each point \(p\in\mathcal{P}_{r,c}\), we consider only admissible pairs that observe it: \[\mathcal{A}(p)=\{(i,j)\in\mathcal{A}\mid i,j\in\mathcal{I}(p)\}.\] If \(\mathcal{A}(p)=\emptyset\), the point contributes no vote.
Let \(C_i\) and \(C_j\) be the camera centers, and let \(P_p\) be the 3D point. Using the normalized viewing rays from \(C_i\) and \(C_j\) toward \(P_p\), we compute the triangulation angle \(\alpha_{ij}(p)\) and define \[T_{ij}(p) = \exp\!\left( -\frac{(\alpha_{ij}(p)-\mu_{\mathrm{triang}})^2}{2\sigma_{\mathrm{triang}}^2} \right).\]
We also project \(P_p\) into images \(i\) and \(j\), compute the normalized radial distances \(r_i(p)\) and \(r_j(p)\) from the principal points, and define centrality weights \[C_i(p) = \exp\!\left( -\frac{r_i(p)^2}{2\sigma_{\mathrm{cent}}^2} \right), \qquad C_j(p) = \exp\!\left( -\frac{r_j(p)^2}{2\sigma_{\mathrm{cent}}^2} \right).\] These are combined symmetrically as \[C_{ij}(p)=\sqrt{C_i(p)\,C_j(p)}.\]
The final point-pair score is \[s_{p,ij} = \operatorname{prior}_{ij} \cdot T_{ij}(p)^{w_{\mathrm{triang}}} \cdot C_{ij}(p)^{w_{\mathrm{cent}}}.\]
Initialize \(S_i=0\) for all \(i\in\mathcal{I}^{\mathrm{cand}}_{r,c}\). To avoid over-counting many near-equivalent pairs for one point, we retain only the top-\(K\) scoring pairs \(\tilde{\mathcal{A}}(p)\subseteq\mathcal{A}(p)\) and normalize their scores: \[w_{p,ij} = \frac{s_{p,ij}}{\sum_{(i,j)\in\tilde{\mathcal{A}}(p)} s_{p,ij} + \varepsilon}.\] Each selected pair casts a unit-normalized vote split equally between its two cameras: \[S_i \leftarrow S_i + \tfrac{1}{2}w_{p,ij}, \qquad S_j \leftarrow S_j + \tfrac{1}{2}w_{p,ij}, \qquad (i,j)\in\tilde{\mathcal{A}}(p).\] Aggregating these votes over all points yields a partition-specific score for every candidate camera. We then sort cameras by \(S_i\) and retain the top \(M\): \[\mathcal{I}^{\mathrm{top}}_{r,c} = \operatorname{TopM}\!\left(\{S_i\}_{i\in\mathcal{I}^{\mathrm{cand}}_{r,c}}, M\right).\]
This section provides the detailed formulation for the dense reconstruction and surface initialization module summarized in [sec:dense95init].
For each partitioned area, we start from its sparse SfM reconstruction and retain the top-\(M\) cameras from the ranking stage. Their intrinsics and poses, \(\{K_n,R_n,t_n\}_{n=1}^{M}\), are kept fixed throughout. We then densify the partition using the depth predictions and pixel correspondences produced by MASt3R. Since MASt3R inference has already been run on image pairs in the original similarity graph, we reuse the cached outputs whenever possible; only pairs newly induced by the final partitioning require fresh inference.
Let \(\mathcal{V}\) denote the selected views and \(\mathcal{E}\) the set of view pairs for which MASt3R correspondences are available after this cache-or-infer step. For each view \(n\in\mathcal{V}\), let \(D_n\) denote the predicted depth map, and let \(\mathcal{M}^{n,m}\) be the set of matched pixels between views \(n\) and \(m\). For a correspondence \(u_c^n \leftrightarrow u_c^m\), we define the associated world-space dense points by back-projecting each pixel using the fixed SfM camera: \[X_c^n = \Pi_n^{-1}(u_c^n, D_n(u_c^n)), \qquad X_c^m = \Pi_m^{-1}(u_c^m, D_m(u_c^m)).\] Thus, \(X_c^n\) is the 3D point induced by pixel \(u_c^n\) from the predicted depth of view \(n\), expressed in the common SfM world frame. Equivalently, the collection \(\{X_c^n\}_c\) can be interpreted as the dense pointmap of view \(n\) in world coordinates, consistent with the pointmap representation used in MASt3R.
Although the cameras are fixed, these per-view dense reconstructions are not necessarily globally consistent: they may still suffer from residual scale drift across views and local depth noise. To address this, we adopt a two-stage alignment strategy inspired by MASt3R-SfM, but specialized to our setting by optimizing only scene geometry while keeping all camera parameters frozen.
We first assign each view a global scale \(s_n>0\) and rescale its dense points as \(\widetilde{X}_c^n = s_n X_c^n\). The scales are estimated by minimizing the discrepancy between matched 3D points: \[\label{eq:supp95coarse95scale} \{s_n^*\} = \arg\min_{\{s_n\}} \sum_{(n,m)\in\mathcal{E}} \sum_{c\in\mathcal{M}^{n,m}} q_c \, \| s_n X_c^n - s_m X_c^m \|_2^{\lambda_1},\tag{4}\] where \(q_c\) denotes the confidence of match \(c\), and \(0<\lambda_1\leq 1\) defines a robust penalty. To remove the gauge ambiguity, we normalize the estimated scales as \[s_n = s_n' / \min_k s_k'.\] This stage allows each view to expand or contract globally so that corresponding dense points become metrically compatible across overlapping views.
A single scale per view is insufficient to correct local depth errors. We therefore further refine the geometry by directly optimizing the dense 3D points while keeping the cameras fixed. Let \(\widehat{X}_c^n\) denote the refined 3D points. We optimize them by minimizing the bidirectional reprojection inconsistency: \[\label{eq:dense95refine} \begin{align} \{\widehat{X}_c^n\} = \arg\min_{\{\widehat{X}_c^n\}} \sum_{(n,m)\in\mathcal{E}} \sum_{c\in\mathcal{M}^{n,m}} q_c \big[ & \rho(u_c^n - \Pi_n(\widehat{X}_c^m)) \\ & + \rho(u_c^m - \Pi_m(\widehat{X}_c^n)) \big], \end{align}\tag{5}\] where \(\Pi_n(\cdot)\) is the projection operator of the fixed camera \(n\), and \(\rho(x)=\|x\|_2^{\lambda_2}\) with \(0<\lambda_2\leq 1\) is a robust reprojection penalty. This objective is analogous to bundle adjustment, except that only the scene geometry is optimized. The coarse stage removes dominant inter-view scale mismatch, while the refinement stage suppresses local geometric noise by enforcing multi-view reprojection consistency.
Both stages are optimized with Adam for a fixed number of iterations, yielding aligned dense depth maps for the partition.
The aligned depth maps are fused using TSDF integration, which aggregates multi-view evidence and regularizes local inconsistencies through volumetric averaging. From the TSDF grid, we extract a fused global point cloud with normals. These oriented points are then passed to Screened Poisson Surface Reconstruction to obtain an initial watertight mesh, which is used to initialize the subsequent differentiable-rendering-based mesh refinement stage.
In summary, this module converts fixed sparse-SfM cameras and mostly cached MASt3R predictions into a dense, globally aligned, and denoised surface initialization, without re-estimating camera poses.
This section provides the full formulation and implementation details omitted from [sec:mesh95refinement].
At iteration \(k\), the current mesh is \(M^k=(V^k,F^k)\). Given calibrated cameras \(\{K_j,T_{wc,j}\}_{j=1}^{N}\), we optimize \[\Phi(M^k)= \sum_{j=1}^{N} \left( \lambda_n\,\mathcal{L}_n^{(j)}(M^k)+ \lambda_s\,\mathcal{L}_{\mathrm{sil}}^{(j)}(M^k) \right)+ \mathcal{R}(M^k),\] where \(\mathcal{R}\) is a lightweight mesh regularizer, e.g.Laplacian smoothing.
For view \(j\), let \(S_j(u)\in[0,1]\) be the target silhouette and \(\hat{S}_j(u;M^k)\) the rendered silhouette. The silhouette loss is \[\mathcal{L}_{\mathrm{sil}}^{(j)}(M^k)= \frac{1}{|\Omega_{\mathrm{img}}|} \sum_{u\in\Omega_{\mathrm{img}}} \big(\hat{S}_j(u;M^k)-S_j(u)\big)^2.\] Let \(\mathcal{N}_j(u)\in\mathbb{S}^2\) be the predicted unit normal map in camera coordinates and \(\hat{\mathcal{N}}_j(u;M^k)\) the rendered unit normal map. On the foreground support \(\Omega_j=\{u\mid S_j(u)=1\}\), the normal-map loss is \[\mathcal{L}_n^{(j)}(M^k)= \frac{1}{|\Omega_j|} \sum_{u\in\Omega_j} \left( 1-\hat{\mathcal{N}}_j(u;M^k)^\top \mathcal{N}_j(u) \right).\]
Gradients of \(\Phi\) are backpropagated through the renderer to vertex positions, which are updated using the IsotropicAdam optimizer of Continuous Remeshing. As in the original framework, the optimizer also provides a relative vertex velocity \(\nu^k(i)\in[0,1]\), which measures the normalized update magnitude at vertex \(i\) and is used as a local indicator of stabilization.
The original controller tends to distribute resolution relatively uniformly, which is inefficient for large scenes: flat regions consume unnecessary triangles, while detailed regions remain under-sampled. We instead decouple where fine resolution is ultimately needed from when that refinement is allowed.
For each view \(j\), let \(\mathbf{n}^{(j)}(u,v)\in\mathbb{S}^2\) be the unit normal field, with image-space Jacobian \[N^{(j)}(u,v)= \begin{bmatrix} \partial_u \mathbf{n}^{(j)}(u,v) & \partial_v \mathbf{n}^{(j)}(u,v) \end{bmatrix} \in\mathbb{R}^{3\times 2}.\] We define \[M^{(j)}(u,v)=\big(N^{(j)}(u,v)\big)^\top N^{(j)}(u,v)\in\mathbb{R}^{2\times 2},\] and the worst-case local normal-rotation rate \[s^{(j)}(u,v)= \sqrt{\lambda_{\max}\!\left(M^{(j)}(u,v)\right)}= \|N^{(j)}(u,v)\|_2.\] Given a normal-rotation tolerance \(\theta_0\) (radians), the target projected edge length in pixels is \[p_{\mathrm{tgt}}^{(j)}(u,v)= \mathrm{clip}\!\left(\frac{\theta_0}{s^{(j)}(u,v)},\,p_{\min},\,p_{\max}\right).\] Hence, regions with rapid normal variation receive smaller projected target edges, while smooth regions permit coarser sampling.
We then lift this per-view pixel-domain target to the mesh surface. For a vertex \(i\) visible in view \(j\), let \(\pi^{(j)}(\mathbf{v}_i)\) be its projection and \(z_i^{(j)}\) its camera-space depth. Sampling the pixel target gives \(p_{ij}=p_{\mathrm{tgt}}^{(j)}(\pi^{(j)}(\mathbf{v}_i))\), which is converted to world units as \[L_{ij}=\frac{z_i^{(j)}}{f^{(j)}}\,p_{ij},\] where \(f^{(j)}\) is the focal length in pixels. Pooling across visible views \(\mathcal{V}_i\) yields the per-vertex curvature-guided baseline \[L_{\mathrm{ref\text{-}curv}}(i)=\operatorname*{median}_{j\in\mathcal{V}_i} L_{ij}.\]
In practice, remeshing acts on edge targets derived from the vertex field by endpoint averaging. This lets us maintain the reference field on vertices while applying split/collapse decisions on edges.
To control when refinement is released, we introduce a nonnegative slack \(\zeta^k(i)\) on top of the curvature-guided baseline. The slack is initialized as \[\zeta^0(i)\leftarrow \max\!\big(0,\;L_{\mathrm{curr}}^0(i)-L_{\mathrm{ref\text{-}curv}}^0(i)\big),\] where \(L_{\mathrm{curr}}^0(i)\) denotes the current local edge scale at initialization.
The slack is then updated using the optimizer’s relative vertex velocity: \[\zeta^{k}(i)\leftarrow \max\!\Big( 0,\; \zeta^{k-1}(i)+g\big(\nu^k(i)-\nu_{\mathrm{ref}}\big) \Big),\] where \(g>0\) is a gain and \(\nu_{\mathrm{ref}}\) is a reference speed. The resulting per-vertex reference edge length is \[L_{\mathrm{ref}}^{k}(i)= \mathrm{clip}\!\Big( L_{\mathrm{ref\text{-}curv}}^{k}(i)+\zeta^k(i),\, L_{\min},\, L_{\max} \Big).\] The corresponding per-edge reference length is obtained from the two incident vertex values. Intuitively, curvature specifies where fine resolution is eventually required, while velocity controls when the local mesh is allowed to approach that target.
The global bounds \(L_{\min}\) and \(L_{\max}\) are estimated automatically from camera intrinsics and a robust scene-depth statistic. Let \(f_x=K_{00}\), \(f_y=K_{11}\), and \(f_{\mathrm{iso}}=\sqrt{f_xf_y}\). Flatten all valid rendered depths of the initial mesh and let \(z_p\) be a robust percentile. Then \[\mathrm{wlpp}=\frac{z_p}{f_{\mathrm{iso}}},\qquad L_{\min}=\mathrm{wlpp}.\] If \(\mathcal{E}\) denotes the set of unique mesh edges and \[L_{\mathrm{med}}=\operatorname{median}_{(i,m)\in\mathcal{E}}\|\mathbf{v}_i-\mathbf{v}_m\|_2,\] we set \[L_{\max}= \max\big(k_1L_{\mathrm{med}},\;k_2L_{\min}\big),\] with \(k_1,k_2>1\). This prevents both sub-pixel triangles and overly coarse edges.
After each vertex update, we invoke the remesher under the current reference field \(L_{\mathrm{ref}}^k\). As in Continuous Remeshing, remeshing is embedded directly into the optimization loop and consists of local collapse, split, and flip operations. Let \(\ell(e)\) denote the current length of edge \(e\). Edges are collapsed when \[\ell(e) < \tau_c\,L_{\mathrm{ref}}^k(e),\] split when \[\ell(e) > \tau_s\,L_{\mathrm{ref}}^k(e),\] and flipped when doing so improves local triangle quality and valence regularity, subject to the same geometric validity checks as in the original framework. Thus, the remesher keeps local edge lengths near the target field while preserving a well-shaped triangulation.
The overall refinement alternates render \(\rightarrow\) compare \(\rightarrow\) update \(\rightarrow\) remesh, and the final mesh is selected as the best snapshot \(M^\star=\arg\min_k \Phi(M^k)\).
We assess mesh quality using the following metrics, which quantify complementary aspects of geometric validity, topological consistency, and structural regularity.
| Metric | Ours | MiLo | MeshSplatting | RadianceMesh |
| 488K V · 973K F | 1384K V · 2771K F | 1190K V · 3503K F | 483K V · 1101K F | |
| Aspect Ratio (AR) \(\downarrow\) | 3.025 | 15.652 | 5.224 | 3.147 |
| Angle bad ratio (ANG) \(\downarrow\) | 1.92% | 28.83% | 15.25% | 2.93% |
| Degenerate Triangle Ratio (DTR) \(\downarrow\) | 0.0000 | 0.0003 | 0.0000 | 0.0000 |
| Non-Manifold Edge Ratio (NME) \(\downarrow\) | 0.0039% | 0.0000% | 27.34% | 1.28% |
| Non-Manifold Vertex Ratio (NMV) \(\downarrow\) | 0.0027% | 0.0014% | 31.34% | 6.26% |
| Vertex Valence Deviation (VVD) \(\downarrow\) | 0.505 | 1.421 | 3.045 | 1.914 |
| Connected Components (CC) \(\downarrow\) | 9 | 918 | 340 | 232 |
| Interior Boundary Loops (IBL) \(\downarrow\) | 24 | 61 | 6554 | 58 |
For a triangular face \(f \in \mathcal{F}\), let \(h_f\) denote its longest edge length and \(r_f\) its inradius. The face-wise aspect ratio is defined as \[\mathrm{AR}(f) = \frac{h_f}{2r_f},\] and the mesh-level aspect ratio is given by \[\mathrm{AR} = \frac{1}{|\mathcal{F}|} \sum_{f \in \mathcal{F}} \mathrm{AR}(f).\] Lower values indicate better-shaped triangles.
Let \(\theta_{\min}(f)\) and \(\theta_{\max}(f)\) denote the minimum and maximum interior angles of face \(f\), respectively. Given angular thresholds \(\theta_{\mathrm{lo}}\) and \(\theta_{\mathrm{hi}}\), the angle bad ratio is defined as \[\mathrm{ANG} = \frac{ \left| \left\{ f \in \mathcal{F} \;\middle|\; \theta_{\min}(f) < \theta_{\mathrm{lo}} \;\text{or}\; \theta_{\max}(f) > \theta_{\mathrm{hi}} \right\} \right| }{ |\mathcal{F}| }.\] This metric measures the fraction of triangles with poor angular quality. Lower values are preferred.
Let \(A(f)\) denote the area of face \(f \in \mathcal{F}\), and let \(\varepsilon_A > 0\) denote a small numerical threshold. The degenerate triangle ratio is defined as \[\mathrm{DTR} = \frac{ |\{ f \in \mathcal{F} \mid A(f) \le \varepsilon_A \}| }{ |\mathcal{F}| }.\] This metric quantifies the proportion of triangles that are degenerate or nearly degenerate. Lower values are better.
Let \(\mathcal{E}\) denote the set of mesh edges, and let \(\deg(e)\) denote the number of incident faces of edge \(e \in \mathcal{E}\). The non-manifold edge ratio is defined as \[\mathrm{NME} = \frac{ |\{ e \in \mathcal{E} \mid \deg(e) > 2 \}| }{ |\mathcal{E}| }.\] Lower values indicate better topological consistency.
Let \(\mathcal{V}\) denote the set of mesh vertices. The non-manifold vertex ratio is defined as \[\mathrm{NMV} = \frac{ |\{ v \in \mathcal{V} \mid v \text{ is non-manifold} \}| }{ |\mathcal{V}| }.\] Lower values indicate fewer vertex-level topological defects.
For each vertex \(v \in \mathcal{V}\), let \(\mathrm{val}(v)\) denote its valence, i.e., the number of one-ring neighboring vertices. We define vertex valence deviation as \[\mathrm{VVD} = \frac{1}{|\mathcal{V}|} \sum_{v \in \mathcal{V}} \left| \mathrm{val}(v) - 6 \right|.\] This metric measures deviation from regular triangular connectivity. Lower values indicate a more regular tessellation.
Let \(\mathcal{C}\) denote the set of connected components of the mesh under standard edge-based connectivity. We define \[\mathrm{CC} = |\mathcal{C}|.\] Lower values indicate a more globally coherent mesh, with \(\mathrm{CC}=1\) corresponding to a single connected surface.
Let \(\mathcal{C}_{\max}\) be the largest connected component of the mesh, and let \(\mathcal{B}(\mathcal{C}_{\max})\) denote its boundary loops. For outdoor scene meshes, which are typically open surfaces, we treat the longest loop as the intended exterior boundary of the dominant surface, and define \[\mathrm{IBL} = \max\!\left( |\mathcal{B}(\mathcal{C}_{\max})| - 1,\; 0 \right).\] IBL counts the remaining boundary loops after excluding this primary exterior boundary. Lower values indicate fewer residual open-boundary artifacts, such as unintended holes, tears, or missing patches. It does not capture handle- or genus-type topological holes that do not induce boundaries.
As shown in 6, our method achieves the best overall mesh quality. It ranks first in aspect ratio (\(3.025\)), angle bad ratio (\(1.92\%\)), vertex valence deviation (\(0.505\)), connected components (\(9\)), and interior boundary loops (\(24\)), while also achieving zero degenerate triangles. This indicates better triangle quality, more regular tessellation, and markedly stronger global coherence than the baselines.
Although MiLo attains slightly lower non-manifold edge and vertex ratios, it remains highly fragmented (\(918\) connected components, \(61\) interior boundary loops). RadianceMesh also achieves zero degenerate triangles, but is substantially worse in connectivity and boundary quality (\(232\) connected components, \(58\) interior boundary loops). MeshSplatting further achieves zero degenerate triangles, yet suffers from severe topological defects, with very high non-manifold ratios (\(27.34\%\) NME, \(31.34\%\) NMV), \(340\) connected components, and \(6554\) interior boundary loops.
These results show that our method yields the most structurally clean and simulation-ready meshes among the compared approaches.
Although our pipeline is highly scalable and produces clean, simulation-ready meshes, a few practical constraints remain. Our method inherits any residual errors from upstream SfM or pretrained depth prediction models, and extremely texture-poor or reflective regions can still challenge local reconstruction quality. Seam merging is robust, but very fine details at partition boundaries may experience mild smoothing. Finally, as clusters are processed independently, the system currently does not leverage global structural priors such as repeated façades or urban symmetries. Overall, these limitations are localized and do not affect the method’s strong scalability, stability, and suitability for city-scale digital-twin applications.