High-quality Surface Reconstruction using Gaussian Surfels

In plain words

Gaussian splatting represents a scene as many soft, blob-shaped pieces. Those blobs are thick, so they do not sit neatly on a surface, and it is unclear which way each one faces. This paper squashes every blob into a flat disk. The disk's face direction is then its normal, which gives the optimiser a clear target.

Flat disks alone do not learn to face the right way, because the photo error barely responds to the face direction. The paper adds a rule that the disk's facing direction must agree with the direction implied by the depth map the disks render. It also uses a normal estimate from a pretrained network to steer the disks in bright, shiny areas.

After training, the method renders depth and normal maps from many viewpoints, trims away stray points near object edges, and merges everything into one mesh with a standard surface-fitting step.

TL;DR

Replacing each 3D Gaussian with a zero-thickness disk gives each piece a well-defined normal. Depth-normal consistency, monocular normal priors, opacity and mask terms train the disks, and volumetric cutting plus screened Poisson produce the mesh. The authors report better surface accuracy than 3DGS and SuGaR on DTU and BlendedMVS, and a mean Chamfer distance that is worse than NeuS2 on DTU. Use this card if you need a fast, explicit surface reconstruction from posed images, and check the licence first (see the header).

Problem

Plain 3D Gaussian splatting has three problems for surface work. The Gaussians have non-zero thickness, so they do not align with a surface §1. The normal of each Gaussian is ambiguous, because its shortest axis can change during optimisation §1. Alpha blending biases depth near sharp edges and can place a surface point in front of the true surface §1 §3.4. Earlier fixes, such as SuGaR's flattening regulariser, need a minimal scale that may shift during training §1.

Method

The method has four parts: a surfel representation §3.2, an optimisation with a set of losses §3.3, a volumetric cutting step that removes erroneous depth points §3.4, and a Poisson meshing step §3.4. Implementation settings are in §3.5.

Gaussian surfels §3.2

1. Each surfel has a centre, a quaternion rotation, an opacity, spherical-harmonic colour coefficients and two in-plane scales. The third scale is set to zero, so the 3D covariance is a flat matrix built from the two in-plane scales and the rotation eq. 1 eq. 2.
2. The surfel normal is the third column of the rotation matrix §3.2.
3. For rendering, each surfel is projected to a 2D ellipse in image space, using the same ray-space approximation as 3DGS eq. 4 eq. 5 §3.2.
4. Colour is alpha-blended front to back from the projected ellipses eq. 3.
5. Depth and normal maps use the same blending weights. The blended result is divided by one minus the final transmittance, so that the weights sum to one. The authors state that this works better for depth and normals than adding a background value eq. 6 §3.2.
6. The depth contributed by each surfel is not its centre depth. It is the depth where the pixel ray meets the ellipse, approximated by a first-order Taylor expansion over the pixel offset eq. 7 §3.2 fig. 3.

Optimisation §3.3

The total loss has five named terms plus a mask term. The paper combines them as photometric, normal prior, opacity, depth-normal consistency and mask, with weights on the last three eq. 8 §3.3.

1. Photometric loss: a weighted sum of an L1 term and a D-SSIM term, with weights 0.8 and 0.2 eq. 9 §3.3.
2. Depth-normal consistency loss: one minus the dot product between the rendered normal and a normal computed from neighbouring points that the rendered depth back-projects to eq. 10 §3.3 fig. 4. The authors say this term is what provides a gradient to the normal direction, which the photometric loss cannot supply §3.3.
3. Normal prior loss: a term that pulls the rendered normal toward a normal map from a pretrained monocular network (Omnidata), with weight 0.04, plus an L1 penalty on the spatial gradient of the rendered normal, with weight 0.005 eq. 11 §3.3. The authors use it mainly in areas with highlights, where the photometric loss can place surface points wrongly §3.3.
4. Opacity loss: a term that peaks when a surfel's opacity is 0.5 and falls to near zero for opacities near 0 or 1, which pushes opacities towards 0 or 1 eq. 12 §3.3.
5. Mask loss: a binary cross-entropy between the accumulated alpha along each ray and the object mask §3.3. The equation for this term is commented out in the source, so it is described only in the text.
6. Trade-off weights: the opacity weight is 0.01, the mask weight is 1, and the consistency weight rises linearly from 0 to 0.1 during training §3.5.
7. Gradient rescaling: the derivative of the rendered normal with respect to the third rotation column is multiplied by 10, so that the normal receives a gradient comparable to the other axes §3.5.

Volumetric cutting and meshing §3.4

1. Render depth and normal maps at every training view §3.4.
2. Build a 512 by 512 by 512 voxel grid over the bounding box. For each surfel, add its Gaussian weight times its opacity into the voxels it covers. The integral is approximated by the value at the voxel centre §3.4.
3. Mark voxels with accumulated value below lambda = 1 as empty. Remove back-projected depth points that fall in empty voxels §3.4 fig. 5 fig. 6.
4. Fuse the remaining oriented points from all views and run screened Poisson reconstruction with tree depth 10 §3.4.

# Volumetric cutting and Poisson meshing  [§3.4]
# Input:  trained surfels, training cameras
# Output: triangle mesh

def extract_mesh(surfels, cameras):
    # Accumulate weighted opacity of all surfels into a 512^3 grid  [§3.4]
    grid = zeros(512, 512, 512)
    for surfel in surfels:
        grid += gaussian_weight(surfel, voxel_centres) * surfel.opacity
    points = []
    for camera in cameras:
        # Render depth and normal maps for this view  [§3.2] [eq. 6]
        depth, normal = render_depth_and_normal(surfels, camera)
        # Back-project each pixel to an oriented 3D point  [§3.4]
        candidates = backproject(depth, normal, camera)
        # Keep points whose voxel has accumulated value above lambda = 1  [§3.4]
        keep = grid_value(grid, candidates) > 1.0
        points.append(candidates[keep])
    # Screened Poisson on the fused oriented points, tree depth 10  [§3.4]
    return screened_poisson(concat(points), depth=10)

Training loop §3.5

1. Use Adam. Each iteration draws one training image and uses all of its pixels §3.5.
2. Run 15k iterations. Start at a lower resolution as a warm-up, without a stated schedule §3.5.
3. Learning rates (exponential decay for position): see the table below §3.5.
4. Density control follows 3DGS: split, clone and prune §3.5.
5. Prune surfels that receive no gradient over N iterations, where N is the number of training images, on the reasoning that such surfels are invisible from all views §3.5.
6. Initialise surfels at random positions and rotations inside an approximate bounding box. SfM points were tested and sped up only the first steps §3.5.

# One training iteration  [§3.3] [§3.5]
# Input:  surfels, one posed image with its mask and monocular normal map
# Output: scalar loss, back-propagated before the optimiser step

def training_step(surfels, camera, image, mask, mono_normal, step):
    # Splat surfels into colour, opacity, normal and depth maps  [§3.2]
    color, opacity, normal, depth = render(surfels, camera)
    # Compute a normal from the rendered depth by neighbouring points  [§3.3]
    depth_normal = normal_from_depth(depth, camera)
    # Photometric term: 0.8 L1 + 0.2 D-SSIM  [eq. 9]
    loss_photo = 0.8 * l1(color, image) + 0.2 * dssim(color, image)
    # Monocular normal prior and gradient smoothness  [eq. 11]
    loss_normal = 0.04 * (1 - dot(normal, mono_normal)).mean()
    loss_normal += 0.005 * l1(normal_gradient(normal), 0)
    # Make rendered normal agree with depth-derived normal  [eq. 10]
    loss_consistency = (1 - dot(normal, depth_normal)).mean()
    # Push each surfel's opacity towards 0 or 1  [eq. 12]
    loss_opacity = exp(-(surfels.opacity - 0.5) ** 2 / 0.05).mean()
    # Binary cross-entropy between accumulated alpha and mask  [§3.3]
    loss_mask = bce(accumulated_alpha(surfels, camera), mask)
    # Weighted sum with opacity 0.01, mask 1, consistency ramp 0 to 0.1  [§3.5]
    lambda_c = 0.1 * min(step / total_steps, 1.0)
    loss = loss_photo + loss_normal + 0.01 * loss_opacity
    loss = loss + lambda_c * loss_consistency + 1.0 * loss_mask
    return loss

Hyperparameters

The paper gives most values in §3.5. Entries marked "code" come from the official repository at the commit below and may differ from the paper (see Paper vs code).

Results

Author-reported. Setup: DTU has 15 scenes with 49 or 64 images at 1600 by 1200 pixels, with object masks §4. BlendedMVS has 18 scenes with 24 to 143 images at 768 by 576 pixels, with masks §4. Baselines are NeuS, INSR and NeuS2, which are implicit methods, and 3DGS and SuGaR, which are Gaussian-point methods. All baselines use the original authors' code and hyperparameters §4. For 3DGS, the authors extracted a mesh with a similar pipeline but took normals from the rendered depth. For SuGaR, they used the lambda = 0.3 level set of its density function and Poisson depth 10 §4. Times are training minutes tab. 1 tab. 2.

DTU, mean Chamfer distance in mm (lower is better) tab. 1

BlendedMVS, mean Chamfer distance in mm (lower is better) tab. 2

DTU rendering quality, PSNR and SSIM tab. 3

Ablation on DTU tab. 4

The paper says the method has a larger Chamfer distance than NeuS2 on DTU and attributes this to alpha-blending bias in per-view depth §4. The paper also says it is better than 3DGS and SuGaR on both datasets §4.

Paper vs code

The repository is turandai/gaussian_surfels at the commit in the frontmatter. Its owner's profile name is "pinxuan turandai", which matches author Pinxuan Dai. The README cites this paper and links the authors' project page README.md. The code builds on the 3DGS codebase and borrows IDR data loading README.md.

• Mask loss. Paper: binary cross-entropy on accumulated alpha, weight 1 §3.3 §3.5. Code: opacity multiplied by one minus the max-pooled mask, weight 0.1 train.py, train.py.
• Consistency weight. Paper: linear ramp from 0 to 0.1 §3.5. Code: 0.01 plus 0.1 times min(2 * iteration / iterations, 1), so it starts at 0.01 and reaches 0.11 at half of training train.py.
• Opacity loss. Paper: exp(-(o - 0.5)^2 / 0.05), weight 0.01 eq. 12 §3.5. Code: the same exponential with 20 in place of 1 / 0.05, multiplied by a mask that gives 0.01 weight for opacities between 0.01 and 0.5 and weight 1 between 0.5 and 0.99, and zero otherwise train.py, train.py.
• Normal prior. Paper: weight 0.04 on the normal term and 0.005 on the L1 gradient penalty eq. 11. Code: the normal term weight decays from 0.04 to 0.02 over training train.py. The smoothness term is a sum of differences to four neighbours, with an L1 norm over channels, not a gradient image_utils.py.
• Normal gradient factor. Paper: multiply the derivative with respect to the third rotation column by 10 §3.5. Code: multiply the gradient on each surfel's normal channel by 10 in the rasterizer backward pass backward.cu. The code does not show the chain to the rotation in this file, so this is a reading, not a verified equivalence.
• Surfel normal sign. Code flips the normal of back-facing surfels by default [code: arguments/init.py#L64], forward.cu. The paper does not mention this.
• Optimiser. Paper: "ADMM training procedure" §3.5. Code: plain Adam on a per-iteration loss, with no ADMM step found gaussian_model.py, train.py.
• Resolution schedule. Paper: lower-resolution warm-up, schedule not stated §3.5. Code: factor 4 at iteration 1, factor 2 at iteration 2002 and factor 1 at iteration 5002 train.py. The factor is used only when the camera stack is next refilled, so a change takes effect at the next refill (inferred from the code, not run).
• Initialisation. Paper: random positions and rotations inside an approximate bounding box §3.5. Code: uses points3d.ply (SfM) when it exists. Otherwise it generates one million random points on a sphere of radius 0.5 in the normalised frame dataset_readers.py. With masks enabled, the code also carves points by the visual hull before training train.py. The paper does not mention the visual-hull step.
• Scales and opacity at start. Not stated in the paper. Code sets scales from the distance to the nearest neighbours gaussian_model.py, sets the third scale to a very negative log value so it is effectively zero gaussian_model.py, and sets initial opacity to 0.1 gaussian_model.py.
• Position learning rate. Paper: 1.6e-4 decaying to 1.6e-6 §3.5. Code multiplies the position rate and its final value by the scene radius gaussian_model.py, and the scene radius is set in [code: scene/init.py#L75]. A delay multiplier of 0.01 is set [code: arguments/init.py#L85]. None of this is in the paper.
• Opacity reset. Not in the paper. Code multiplies every opacity by 0.12 every 1000 iterations, rather than capping it train.py, gaussian_model.py. The commented-out capped version is in the same file gaussian_model.py.
• Background during training. Not in the paper. The default is a random background colour per iteration [code: arguments/init.py#L55], train.py.
• Spherical harmonics. Not in the paper. Maximum degree 3, raised one level every 1000 iterations [code: arguments/init.py#L49], train.py.
• Density control. Paper: follows 3DGS, no thresholds stated §3.5. Code: starts at iteration 200, gradient threshold 0.0001, and prunes opacity below 0.1 [code: arguments/init.py#L98-L100], train.py, train.py. The no-gradient prune runs every max(100, number of images) iterations gaussian_model.py, train.py.
• Volumetric grid bounds. Paper: bounding box, no margin stated §3.4. Code: expands the box by 10 percent on each side and sets the voxel size from the longest side divided by 512 gaussian_model.py.
• Cutting threshold. Paper: lambda = 1 §3.4. Code: the grid is built with cutoff 0.0 and points are kept when the grid value is above 1 render.py, render.py. The meaning of the 0.0 argument is not traced in the repository (the kernel is in a compiled extension not included here).
• Depth range for meshing. Not in the paper. Code keeps only depths between 0 and 20 render.py.
• Poisson pruning. Not in the paper. Code prunes the Poisson mesh with threshold 3e-5 render.py. The tree depth 10 matches the paper §3.4 and the README README.md.
• Monocular normals. Paper: Omnidata normal maps §3.3. Code: loads precomputed normals from a folder next to the images and uses them when mono_normal is set, which is the default dataset_readers.py, [code: arguments/init.py#L65].

This is the official code: its owner matches author Pinxuan Dai and its README cites this paper. The arXiv abstract page lists no code link.

Pitfalls

• Normal direction is learned indirectly. The photometric loss has no gradient with respect to the third rotation column. Only the depth-normal consistency loss and the normal prior move the normal §3.3. If you drop either, expect the normals to drift.
• Reduction of the opacity loss is not stated. The paper gives the per-surfel term but not whether it is averaged over surfels or pixels eq. 12. The code averages over surfels and applies a mask train.py.
• Several optimiser details are missing. The paper does not state the spherical-harmonic degree, the initial opacity, the initial scales, the density-control thresholds or the learning-rate delay. The code uses values listed in Paper vs code.
• Opacity reset is multiplicative in the code. It scales all opacities by 0.12 every 1000 iterations, so a surfel that starts at 0.9 drops to about 0.11 gaussian_model.py. Check this before copying the training schedule.
• Cutting depends on scene scale. The voxel size comes from the bounding box, so a different scene size changes what counts as an outlier. The threshold of 1 is relative to that grid §3.4 gaussian_model.py.
• Specular regions and weak texture. The authors say the method cannot guarantee accuracy on strong specular reflections, and that very weak textures can show a global shift §5, Limitations.
• Cost. Reported training time is 6.67 minutes on DTU and 3.82 minutes on BlendedMVS tab. 1 tab. 2, on an RTX 4090 §3.5. The volumetric cut and Poisson steps add to this, and the paper does not time them separately. The 512 cubed grid alone needs about 0.5 GB at float32 (computed from its size, not measured).
• Evaluation caveat. BlendedMVS uses all images for training, so there is no held-out test set for PSNR there §4. DTU geometry follows NeuS2 in holding out 7 to 8 images, and the rendering test holds out 12.5% or 50% of images §4.
• Licence of the code. The repository has no LICENSE file. The training script and the rasteriser carry an Inria notice that the software is free for non-commercial, research and evaluation use train.py. Do not copy code from this repository into a commercial product without a clear licence. This card cites code by path and line only and quotes none of it.
• Licence of the normal prior. The normal prior uses Omnidata §3.3. The wiki's monocular-depth-estimation concept page states that Omnidata's data and weights are non-commercial only. Check that before reusing the normal prior commercially.
• Figures. No paper figure is reproduced here. They are described in words and cited by number: fig. 3 for ray-ellipse depth, fig. 4 for the depth-normal consistency, fig. 5 for the cutting problem, fig. 6 for cutting comparisons §3.2 §3.4.

Lineage

• Builds on: 3D Gaussian Splatting (Kerbl et al.). The rasteriser, learning-rate schedule and density control come from it §3.2 §3.5 README.md. No card yet.
• Builds on: surfels (Pfister et al., 2000) for the surface-aligned representation §1. No card yet.
• Compared against: NeuS (implicit surface), NeuS2, INSR, SuGaR and 3DGS §4. NeuS has a card; the others do not yet.
• Related card: 2D Gaussian Splatting is a flat-disk Gaussian method. This paper does not compare against it §4.
• Superseded by: not known from content in this wiki.