> ## Documentation Index
> Fetch the complete documentation index at: https://aegean.ai/llms.txt
> Use this file to discover all available pages before exploring further.

# From splats to a watertight mesh

> Real-to-sim twin stage 3, geometry

[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/pantelis/eng-ai-agents/blob/main/notebooks/projects/cv/real-to-sim-twin/03-geometry/03-geometry.ipynb)

In this section you turn the posed, masked frames of stages 1 and 2 into a closed triangle mesh of your object, in meters. Stage 4 paints that mesh, and stage 5 computes its volume, mass distribution and collision shape from it.

You train two Gaussian splats on the same frames. The first is a 3D Gaussian splat (Kerbl et al., 2023), which renders photos well and is what stage 6 scores for appearance. The second is a 2D Gaussian splat (Huang et al., "2D Gaussian Splatting for Geometrically Accurate Radiance Fields", SIGGRAPH 2024), whose primitives are flat disks, so its rendered depth lies on a surface. Both come from the gsplat library (Apache-2.0) (Ye et al., 2024). You then render the 2D splat's depth from every training view and fuse the depth maps into a truncated signed distance field, the volumetric method of Curless and Levoy ("A volumetric method for building complex models from range images", SIGGRAPH 1996), using Open3D (MIT license). The zero level set of that field is the mesh.

Runtime budget on a Colab T4: about 7 minutes for 90 training frames once the session is set up, about 6 of them spent training the two splats. The first run in a new session takes about 19 minutes, of which about 12 minutes is one-time setup.

```python theme={null}
import json
import shutil
import time

import cv2
import numpy as np

T0 = time.time()
FALLBACK = ensure_inputs(STAGE)
STAGE1, STAGE2 = PROJECT / "stage1", PROJECT / "stage2"
REFERENCE = on_reference(STAGE)  # the reference stage2/ carries a .reference marker
frames = sorted((STAGE1 / "frames").glob("*.jpg"))
heldout = set((STAGE1 / "heldout.txt").read_text().split())
drifting = set(json.loads((STAGE2 / "drift.json").read_text())["flagged"])
train = [p for p in frames if p.name not in heldout and p.name not in drifting]
test = [p for p in frames if p.name in heldout]
masks = {p.name: cv2.imread(str(STAGE2 / "masks" / (p.stem + ".png")), cv2.IMREAD_GRAYSCALE) > 127 for p in frames}
SCAN_OBJECT = json.loads((PROJECT / "ground_truth" / "scan.json").read_text())["object"] if REFERENCE else None
print("object:", f"reference ({SCAN_OBJECT})" if REFERENCE else OBJECT)
print(f"{len(frames)} frames: {len(train)} for training, {len(test)} held out, {len(drifting)} left out for drifting masks")
```

```output theme={null}
stage2/ is missing or empty: running on the reference object (stage2.tar.gz).
```

```output theme={null}
object: reference (OmniObject3D boxed_beverage_001)
100 frames: 90 for training, 10 held out, 0 left out for drifting masks
```

## Two splats of the same object

Both splats train for 7,000 steps at half resolution on the training frames. Frames held out in stage 1 stay out, and so do frames whose masks the stage 2 drift check flagged. The loss is the masked loss of stage 2: every pixel outside the mask takes a background color drawn at random at each step, so the only way to match it is to leave that part of the image empty. Each splat starts from the sparse points of stage 1, and its spherical harmonics grow by one degree every 1,000 steps up to degree 3.

A 3D Gaussian is a blob with center $\mu$ and covariance $\Sigma = R S S^\top R^\top$. To render it, the rasterizer projects it to a 2D Gaussian with covariance $\Sigma' = J W \Sigma W^\top J^\top$, where $W$ is the view rotation and $J$ is the Jacobian of the projection at $\mu$, a local affine approximation. A pixel then reads the 2D Gaussian's value at its own position. The projection never computes where a ray meets the blob, and the blob has thickness along every axis, so there is no surface for a ray to meet.

A 2D Gaussian is a disk. It has a center $p$, two tangent directions $t_u$ and $t_v$ with scales $s_u$ and $s_v$, and its points are

$$
P(u, v) = p + s_u t_u\, u + s_v t_v\, v, \qquad \mathcal{G}(u, v) = \exp\!\left(-\tfrac{u^2 + v^2}{2}\right).
$$

To render it, the rasterizer intersects each pixel's ray with the disk's plane, solves for the $(u, v)$ where the ray hits it, and weights the disk by $\mathcal{G}(u, v)$. The depth of that hit is exact, and the disk's normal $t_u \times t_v$ is defined. After 3,000 steps the 2D splat adds two terms to its loss. The depth distortion term $\sum_{i,j} w_i w_j \lvert z_i - z_j \rvert$ pulls the disks a ray passes through toward one depth. The normal consistency term $\sum_i w_i (1 - n_i^\top N)$ asks each disk's normal $n_i$ to agree with the normal $N$ computed from the gradient of the rendered depth.

The cell prints each splat's training time, its number of Gaussians, and its PSNR over the object's pixels in the held-out frames.

```python theme={null}
import contextlib
import io

import cv2
import numpy as np
import torch
import torch.nn.functional as F
from pathlib import Path
from scipy.spatial import cKDTree


def load_views(stage1: Path, frames, masks=None, data_factor: int = 2):
    """Undistorted images (and masks) at 1/data_factor resolution, with their cameras.

    stage1/sparse/0 is already in meters in the world frame, so the splats are too.
    """
    from pycolmap import SceneManager

    with contextlib.redirect_stdout(io.StringIO()):  # the reader prints a note about the image folder
        model = SceneManager(str(Path(stage1) / "sparse" / "0"))
    model.load_cameras()
    model.load_images()
    model.load_points3D()
    by_name = {im.name: im for im in model.images.values()}
    views, skipped = [], []
    for p in frames:
        im = by_name.get(Path(p).name)
        if im is None:
            skipped.append(Path(p).name)
            continue
        cam = model.cameras[im.camera_id]
        K = np.array([[cam.fx, 0, cam.cx], [0, cam.fy, cam.cy], [0, 0, 1.0]])
        K[:2] /= data_factor
        dist = np.array([getattr(cam, k, 0.0) for k in ("k1", "k2", "p1", "p2")])
        img = cv2.cvtColor(cv2.imread(str(p)), cv2.COLOR_BGR2RGB)
        size = (img.shape[1] // data_factor, img.shape[0] // data_factor)
        img = cv2.undistort(cv2.resize(img, size, interpolation=cv2.INTER_AREA), K, dist)
        view = {"name": Path(p).name, "K": K, "image": img}
        if masks is not None:
            m = cv2.resize(masks[Path(p).name].astype(np.float32), size, interpolation=cv2.INTER_AREA)
            view["mask"] = cv2.undistort(m, K, dist) > 0.5
        viewmat = np.eye(4)
        viewmat[:3, :3], viewmat[:3, 3] = im.R(), im.tvec
        view["viewmat"] = viewmat
        views.append(view)
    if skipped:  # write_stage1 accepts a model with 80% of the frames registered
        print(f"{len(skipped)} of {len(frames)} frames have no pose in stage1/sparse/0 (stage 1 did not register "
              f"them) and are left out: {', '.join(skipped)}")
    return views, model.points3D, model.point3D_colors


def _ssim(x, y, size=11, sigma=1.5):
    """Structural similarity of two H x W x 3 images in [0, 1], Gaussian window."""
    g = torch.exp(-(torch.arange(size, device=x.device) - size // 2) ** 2 / (2 * sigma**2))
    g = (g / g.sum())[None, None]
    win = (g.transpose(-1, -2) @ g).expand(3, 1, size, size)
    x, y = x.permute(2, 0, 1)[None], y.permute(2, 0, 1)[None]
    blur = lambda z: F.conv2d(z, win, groups=3)
    mx, my = blur(x), blur(y)
    vx, vy, cxy = blur(x * x) - mx**2, blur(y * y) - my**2, blur(x * y) - mx * my
    c1, c2 = 0.01**2, 0.03**2
    return (((2 * mx * my + c1) * (2 * cxy + c2)) / ((mx**2 + my**2 + c1) * (vx + vy + c2))).mean()


def _rasterize(splat, viewmat, K, w, h, background, sh_degree, kind):
    """One render of a 3DGS or 2DGS splat. Returns color, alpha, depth (all [H, W, .]) and extras."""
    from gsplat import rasterization, rasterization_2dgs

    colors = torch.cat([splat["sh0"], splat["shN"]], 1)
    args = dict(means=splat["means"], quats=F.normalize(splat["quats"], dim=-1), scales=torch.exp(splat["scales"]),
                opacities=torch.sigmoid(splat["opacities"]), colors=colors, viewmats=viewmat, Ks=K, width=w,
                height=h, sh_degree=sh_degree, backgrounds=background[None], packed=False)
    if kind == "3dgs":
        out, alpha, info = rasterization(**args, render_mode="RGB+ED")
        return out[0, ..., :3], alpha[0, ..., 0], out[0, ..., 3], info, None
    out, alpha, normals, surf_normals, distort, median, info = rasterization_2dgs(**args, render_mode="RGB+ED")
    extras = {"normals": normals[0], "surf_normals": surf_normals[0], "distort": distort[0]}
    return out[0, ..., :3], alpha[0, ..., 0], median[0, ..., 0], info, extras


def train_splat(views, points, colors, kind: str, steps: int = 7000, seed: int = 0,
                reg_start: int = 3000, normal_weight: float = 0.05, distort_weight: float = 100.0):
    """Train a 3DGS (kind="3dgs") or 2DGS (kind="2dgs") splat on masked views.

    Each target pixel outside the object's mask is replaced by a random background color, drawn
    anew every step and also used behind the splat, so the splat empties out off the object.
    Loss: 0.8 L1 + 0.2 (1 - SSIM). For 2DGS, after reg_start steps, add the normal consistency
    and depth distortion terms of 2DGS. Spherical harmonics grow by one degree every 1,000 steps
    up to degree 3.
    """
    from gsplat.strategy import DefaultStrategy

    torch.manual_seed(seed)
    rng = np.random.default_rng(seed)
    dev = "cuda"
    data = [(torch.tensor(v["image"], dtype=torch.float32, device=dev) / 255,
             torch.tensor(v["mask"], dtype=torch.float32, device=dev)[..., None],
             torch.tensor(v["viewmat"], dtype=torch.float32, device=dev)[None],
             torch.tensor(v["K"], dtype=torch.float32, device=dev)[None]) for v in views]
    h, w = views[0]["image"].shape[:2]

    d, _ = cKDTree(points).query(points, k=4)  # one Gaussian per sparse point, sized by its neighbors
    scale0 = np.log(np.maximum(d[:, 1:].mean(1), 1e-6))
    n = len(points)
    centers = np.array([-v["viewmat"][:3, :3].T @ v["viewmat"][:3, 3] for v in views])
    scene_scale = 1.1 * np.linalg.norm(centers - centers.mean(0), axis=1).max()
    splat = torch.nn.ParameterDict({
        "means": torch.nn.Parameter(torch.tensor(points, dtype=torch.float32)),
        "scales": torch.nn.Parameter(torch.tensor(scale0, dtype=torch.float32)[:, None].repeat(1, 3)),
        "quats": torch.nn.Parameter(torch.rand(n, 4)),
        "opacities": torch.nn.Parameter(torch.logit(torch.full((n,), 0.1))),
        "sh0": torch.nn.Parameter(((torch.tensor(colors, dtype=torch.float32) / 255 - 0.5) / 0.28209479)[:, None]),
        "shN": torch.nn.Parameter(torch.zeros(n, 15, 3)),
    }).to(dev)
    lrs = {"means": 1.6e-4 * scene_scale, "scales": 5e-3, "quats": 1e-3, "opacities": 5e-2,
           "sh0": 2.5e-3, "shN": 2.5e-3 / 20}
    optimizers = {k: torch.optim.Adam([{"params": splat[k], "lr": lr, "name": k}], eps=1e-15)
                  for k, lr in lrs.items()}
    decay = torch.optim.lr_scheduler.ExponentialLR(optimizers["means"], gamma=0.01 ** (1 / steps))
    strategy = DefaultStrategy(verbose=False, key_for_gradient="means2d" if kind == "3dgs" else "gradient_2dgs")
    strategy.check_sanity(splat, optimizers)
    state = strategy.initialize_state(scene_scale=scene_scale)

    for step in range(steps):
        image, mask, viewmat, K = data[rng.integers(len(data))]
        background = torch.rand(3, device=dev)
        rendered, alpha, _, info, extras = _rasterize(splat, viewmat, K, w, h, background,
                                                      min(step // 1000, 3), kind)
        target = image * mask + background * (1 - mask)
        loss = 0.8 * (rendered - target).abs().mean() + 0.2 * (1 - _ssim(rendered, target))
        if kind == "2dgs" and step >= reg_start:
            surf = extras["surf_normals"] * alpha.detach()[..., None]
            loss = loss + normal_weight * (1 - (extras["normals"] * surf).sum(-1)).mean()
            loss = loss + distort_weight * extras["distort"].mean()
        strategy.step_pre_backward(splat, optimizers, state, step, info)
        loss.backward()
        for opt in optimizers.values():
            opt.step()
            opt.zero_grad(set_to_none=True)
        decay.step()
        strategy.step_post_backward(splat, optimizers, state, step, info, packed=False)
    return {k: v.detach() for k, v in splat.items()}


@torch.no_grad()
def render(splat, view, kind: str, background=(1.0, 1.0, 1.0)):
    """Render color (H x W x 3 in [0, 1]), opacity and depth in meters of a view, as tensors.

    For 3DGS the depth is the opacity-weighted mean depth along the ray. For 2DGS it is the
    median depth: the depth where the accumulated opacity along the ray crosses one half.
    """
    h, w = view["image"].shape[:2] if "image" in view else view["size"]
    dev = splat["means"].device
    color, alpha, depth, _, _ = _rasterize(
        splat, torch.tensor(view["viewmat"], dtype=torch.float32, device=dev)[None],
        torch.tensor(view["K"], dtype=torch.float32, device=dev)[None], w, h,
        torch.tensor(background, dtype=torch.float32, device=dev), 3, kind)
    return color.clamp(0, 1), alpha, depth


def render_depth(splat2d, view):
    """2DGS median depth (meters), opacity and color of a view, as numpy arrays."""
    color, alpha, depth = render(splat2d, view, "2dgs")
    return depth.cpu().numpy(), alpha.cpu().numpy(), (color.cpu().numpy() * 255).astype(np.uint8)


@torch.no_grad()
def masked_psnr(splat, views, kind: str) -> float:
    """Mean PSNR over the object's pixels (its mask) in the given views."""
    out = []
    for v in views:
        rendered, _, _ = render(splat, v, kind)
        image = torch.tensor(v["image"], dtype=torch.float32, device=rendered.device) / 255
        m = torch.tensor(v["mask"], device=rendered.device)
        out.append(float(-10 * torch.log10(((rendered - image)[m] ** 2).mean())))
    return float(np.mean(out))


def save_ply(splat, path) -> None:
    """Write a splat as a standard 3D Gaussian splatting PLY file."""
    from gsplat import export_splats

    export_splats(splat["means"], splat["scales"], F.normalize(splat["quats"], dim=-1), splat["opacities"],
                  splat["sh0"], splat["shN"], format="ply", save_to=str(path))
```

```python theme={null}
train_views, points, point_colors = load_views(STAGE1, train, masks, data_factor=2)
test_views, _, _ = load_views(STAGE1, test, masks, data_factor=2)
t = time.time()
_probe = {"means": torch.zeros(1, 3, device="cuda"), "quats": torch.ones(1, 4, device="cuda"),
          "scales": torch.zeros(1, 3, device="cuda"), "opacities": torch.zeros(1, device="cuda"),
          "sh0": torch.zeros(1, 1, 3, device="cuda"), "shN": torch.zeros(1, 15, 3, device="cuda")}
render(_probe, train_views[0], "3dgs"), render(_probe, train_views[0], "2dgs")  # the first calls build the GPU kernels
print(f"splatting kernels ready in {time.time() - t:.0f} s")

STEPS = 7000
splats, stats = {}, {}
for kind in ("3dgs", "2dgs"):
    t = time.time()
    splats[kind] = train_splat(train_views, points, point_colors, kind, steps=STEPS)
    stats[kind] = {"gaussians": len(splats[kind]["means"]), "train_s": time.time() - t,
                   "heldout_psnr": masked_psnr(splats[kind], test_views, kind)}
    print(f"{kind.upper()}: {STEPS:,} steps in {stats[kind]['train_s']:.0f} s, {stats[kind]['gaussians']:,} Gaussians, "
          f"held-out PSNR on the object {stats[kind]['heldout_psnr']:.2f} dB")
```

```output theme={null}
splatting kernels ready in 461 s
```

```output theme={null}
3DGS: 7,000 steps in 81 s, 54,361 Gaussians, held-out PSNR on the object 27.56 dB
```

```output theme={null}
2DGS: 7,000 steps in 100 s, 42,739 Gaussians, held-out PSNR on the object 26.91 dB
```

## Depth from each splat

A splat gives depth the same way it gives color, by blending along the ray. The 3D splat renders the expected depth $D = \sum_i w_i z_i / \sum_i w_i$, where $z_i$ is the depth of blob $i$'s center and $w_i$ its blending weight. Overlapping blobs at different depths average into a depth that need not lie on any of them. The 2D splat renders the median depth instead: the depth of the hit at which the accumulated opacity along the ray first passes one half.

Two numbers compare the two kinds of depth. The first is the disagreement between views. Take the object's pixels in one training view, lift them to 3D with their rendered depth, and project them into a view three frames later. If both depth maps describe the same surface, each lifted point lands at the depth the second view rendered there. The cell prints the median and the 90th percentile of the gap in millimeters, ignoring gaps above 2 cm, which belong to points hidden in the second view. The second number is the flatness of the 3D Gaussians: the median ratio of each Gaussian's smallest scale to its largest.

The figure shows a held-out frame, the two depth maps, and the normals computed from each depth map's gradient. A depth map that lies on a smooth surface gives smooth normals.

```python theme={null}
import numpy as np
import open3d as o3d


def backproject(depth, K, viewmat, valid):
    """World points (N x 3) of the valid pixels of a depth map."""
    v, u = np.nonzero(valid)
    z = depth[v, u]
    cam = np.stack([(u + 0.5 - K[0, 2]) / K[0, 0] * z, (v + 0.5 - K[1, 2]) / K[1, 1] * z, z], 1)
    R, t = viewmat[:3, :3], viewmat[:3, 3]
    return (cam - t) @ R  # R^T (x - t), row by row


def drop_too_deep(depths, valid, cameras, tol_m: float, step: int = 2):
    """Drop depth pixels that lie behind the surface the neighboring cameras see.

    depths, valid: per camera, an H x W depth map in meters and the mask of pixels to trust.
    cameras are in orbit order. Each valid pixel of camera i is lifted to 3D and projected into
    cameras i - step and i + step. If, in every neighbor where it lands on a valid pixel, it is
    more than tol_m behind that neighbor's depth, the neighbor sees empty space in front of it:
    the depth of camera i is too deep there, and the pixel is dropped. Returns the new masks.
    """
    out = []
    for i, cam in enumerate(cameras):
        X = backproject(depths[i], cam["K"], cam["viewmat"], valid[i])
        behind = np.zeros(len(X), int)
        seen = np.zeros(len(X), int)
        for j in (i - step, i + step):
            if not 0 <= j < len(cameras):
                continue
            V = cameras[j]["viewmat"]
            c = X @ V[:3, :3].T + V[:3, 3]
            front = c[:, 2] > 1e-6
            uv = c[front] @ cameras[j]["K"].T
            u, v = np.floor(uv[:, 0] / uv[:, 2]).astype(int), np.floor(uv[:, 1] / uv[:, 2]).astype(int)
            h, w = depths[j].shape
            idx = np.flatnonzero(front)
            ok = (u >= 0) & (u < w) & (v >= 0) & (v < h)
            idx, u, v = idx[ok], u[ok], v[ok]
            on = valid[j][v, u]
            idx, u, v = idx[on], u[on], v[on]
            seen[idx] += 1
            behind[idx] += c[idx, 2] > depths[j][v, u] + tol_m
        keep = valid[i].copy()
        ys, xs = np.nonzero(valid[i])
        keep[ys, xs] = ~((seen > 0) & (behind == seen))
        out.append(keep)
    return out


def object_box(points, lo_pct: float = 0.5, hi_pct: float = 99.5):
    """Axis-aligned box (lo, hi) of a point cloud, ignoring the extreme 0.5% on each side."""
    return np.percentile(points, lo_pct, axis=0), np.percentile(points, hi_pct, axis=0)


def is_watertight(mesh) -> bool:
    """True for a closed, consistently oriented surface: every edge has exactly two triangles,
    every vertex a single fan of triangles, and the two triangles at each edge traverse it in
    opposite directions, so all normals point the same way (out or in). Open3D's is_watertight
    also tests every pair of triangles for crossings, which takes minutes on a mesh of this size;
    marching cubes on one grid does not produce them."""
    if len(mesh.triangles) == 0 or not (mesh.is_edge_manifold(allow_boundary_edges=False)
                                        and mesh.is_vertex_manifold()):
        return False
    t = np.asarray(mesh.triangles)
    directed = np.concatenate([t[:, [0, 1]], t[:, [1, 2]], t[:, [2, 0]]])
    return len(np.unique(directed, axis=0)) == len(directed)  # each directed edge used once


def near_grid_border(mesh, origin, length: float, voxel_m: float) -> bool:
    """True if the mesh comes within two voxels of the grid's faces.

    The closing step forces the outermost voxels outside, so a surface that reaches them was
    cut off there: some region was outside every camera's view and was filled in as inside.
    """
    v = np.asarray(mesh.vertices)
    lo, hi = np.asarray(origin), np.asarray(origin) + length
    return bool(len(v) and (np.any(v - lo < 2 * voxel_m) or np.any(hi - v < 2 * voxel_m)))


def crossing_fraction(mesh, n: int = 48) -> float:
    """Share of rays through the mesh that cross its surface more than twice.

    An n x n grid of rays runs along each of x, y and z through the mesh's bounding box, kept
    10% of the box away from its faces so that roughness right at the walls counts less. A ray that
    enters and leaves one shell crosses it twice. More crossings mean a pocket inside the
    object, a deep dent or fold in its wall, or a part of the object that is not convex.
    """
    scene = o3d.t.geometry.RaycastingScene()
    scene.add_triangles(o3d.t.geometry.TriangleMesh.from_legacy(mesh))
    bb = mesh.get_axis_aligned_bounding_box()
    lo, hi = np.asarray(bb.min_bound), np.asarray(bb.max_bound)
    s = np.linspace(0.1, 0.9, n)
    counts = []
    for axis in range(3):
        a, b = [i for i in range(3) if i != axis]
        A, B = np.meshgrid(lo[a] + s * (hi[a] - lo[a]), lo[b] + s * (hi[b] - lo[b]), indexing="ij")
        rays = np.zeros((A.size, 6), np.float32)
        rays[:, a], rays[:, b] = A.ravel(), B.ravel()
        rays[:, axis] = lo[axis] - 0.01 * (hi[axis] - lo[axis]) - 1e-3
        rays[:, 3 + axis] = 1.0
        counts.append(scene.count_intersections(o3d.core.Tensor(rays)).numpy())
    c = np.concatenate(counts)
    c = c[c > 0]
    return float((c > 2).mean()) if len(c) else 0.0


def mesh_volume(mesh) -> float:
    """Enclosed volume of a closed, consistently oriented mesh, in cubic meters (divergence theorem)."""
    v = np.asarray(mesh.vertices)[np.asarray(mesh.triangles)]
    return float(abs(np.einsum("ij,ij->i", v[:, 0], np.cross(v[:, 1], v[:, 2])).sum()) / 6)
```

```output theme={null}
Jupyter environment detected. Enabling Open3D WebVisualizer.
[Open3D INFO] WebRTC GUI backend enabled.
[Open3D INFO] WebRTCWindowSystem: HTTP handshake server disabled.
```

```python theme={null}
import numpy as np
import torch


def depth_agreement(splat, views, kind: str, pairs: int = 20) -> dict:
    """How well the rendered depth of one view predicts the rendered depth of another.

    For pairs of views a few frames apart, the object's pixels of view a are lifted to 3D with
    their depth and projected into view b. Where they land on the object in b and are not
    hidden, the gap |z_a->b - depth_b| is a multi-view disagreement: zero for a depth map
    that describes one surface. Returns its median and 90th percentile in millimeters.
    """
    gaps = []
    idx = np.linspace(0, len(views) - 4, pairs).astype(int)
    for i in idx:
        a, b = views[i], views[i + 3]
        _, alpha_a, depth_a = render(splat, a, kind)
        _, alpha_b, depth_b = render(splat, b, kind)
        depth_a, depth_b = depth_a.cpu().numpy(), depth_b.cpu().numpy()
        ok_a = a["mask"] & (alpha_a.cpu().numpy() > 0.5)
        X = backproject(depth_a, a["K"], a["viewmat"], ok_a)
        cam = X @ b["viewmat"][:3, :3].T + b["viewmat"][:3, 3]
        cam = cam[cam[:, 2] > 1e-6]  # only points in front of camera b project into it
        uv = cam @ b["K"].T
        u, v = (uv[:, 0] / uv[:, 2]).astype(int), (uv[:, 1] / uv[:, 2]).astype(int)
        h, w = depth_b.shape
        inside = (u >= 0) & (u < w) & (v >= 0) & (v < h)
        u, v, z = u[inside], v[inside], cam[inside, 2]
        on = b["mask"][v, u] & (alpha_b.cpu().numpy()[v, u] > 0.5)
        gap = np.abs(z[on] - depth_b[v[on], u[on]])
        gaps.append(gap[gap < 0.02])  # beyond 2 cm the point is hidden in view b, not a disagreement
    g = 1000 * np.concatenate(gaps)
    return {"median_mm": float(np.median(g)), "p90_mm": float(np.percentile(g, 90))}


def flatness(splat) -> float:
    """Median ratio of each Gaussian's smallest scale to its largest (0 for a flat disk)."""
    s = torch.exp(splat["scales"]).cpu().numpy()
    return float(np.median(s.min(1) / s.max(1)))
```

```python theme={null}
agreement = {}
for kind in ("3dgs", "2dgs"):
    agreement[kind] = depth_agreement(splats[kind], train_views, kind)
    print(f"{kind.upper()}: depth disagreement between views, median {agreement[kind]['median_mm']:.2f} mm, "
          f"90th percentile {agreement[kind]['p90_mm']:.2f} mm")
print(f"3DGS Gaussians: median smallest-to-largest scale ratio {flatness(splats['3dgs']):.3f}; "
      "2DGS Gaussians are flat disks, so the ratio there is 0 by construction")
```

```output theme={null}
3DGS: depth disagreement between views, median 0.26 mm, 90th percentile 1.32 mm
```

```output theme={null}
2DGS: depth disagreement between views, median 0.12 mm, 90th percentile 0.75 mm
3DGS Gaussians: median smallest-to-largest scale ratio 0.099; 2DGS Gaussians are flat disks, so the ratio there is 0 by construction
```

![Output from cell 8](https://mintcdn.com/aegeanaiinc/kXOaaxR5t9cRJZkL/aiml-common/projects/cv/real-to-sim-twin/03-geometry/images/cell_8_output_1.jpg?fit=max\&auto=format\&n=kXOaaxR5t9cRJZkL\&q=85\&s=52a7d58edb5b7a36b501e8fed57f2937)

## Fusing depth into a signed distance field

A truncated signed distance field stores, at every voxel $x$ of a grid, an estimate of the signed distance to the nearest surface, positive in front of it and negative behind it. Each depth map gives one estimate along its own rays. For a camera with depth map $D$ that sees $x$ at pixel $\pi(x)$ and depth $z(x)$,

$$
f(x) = \min\!\left(1, \frac{D(\pi(x)) - z(x)}{\mu}\right) \quad \text{if } D(\pi(x)) - z(x) > -\mu,
$$

and the camera says nothing about $x$ when $x$ lies more than the truncation distance $\mu$ behind the surface it sees. The field is the running weighted mean of these estimates over all cameras, $F(x) = \sum_k w_k f_k(x) / \sum_k w_k$, and the surface is its zero level set, extracted as triangles by marching cubes. Averaging many depth maps this way cancels much of their noise when they agree on one surface.

The masks enter the fusion twice. Inside the mask, eroded by two pixels so that edge pixels do not count, the 2D splat's median depth is a measurement wherever the splat is opaque. Outside the mask the ray misses the object, so you set its depth far beyond the grid: every voxel along that ray gets $f = 1$, empty. This carves away everything outside the object's visual hull.

A grid filled this way still has holes. No camera saw the bottom of the object, which rests on the support plane $z = 0$, and no camera saw the inside. After fusion you close the field. A voxel near the object that no camera measured lies behind the surface in every view, so it is inside ($F = -1$). Far from the object the same rule goes wrong. The grid is a cube sized by the object's longest side, so above a long, low object, such as a hammer lying flat, it holds a large block of space outside every camera's view. Set to inside, that block would become the largest solid in the grid and replace the object. A voxel that no camera measured is therefore inside only within the object's box, grown by three voxels on every side, and outside beyond it. Voxels below the support plane and on the border of the grid are outside ($F = +1$). Last, empty voxels that cannot be reached from the border of the grid through other empty voxels are enclosed by the surface, so they are filled in as inside too. Marching cubes on a grid that is outside all around its border then gives one closed outer shell with nothing inside it. The bottom of the mesh is therefore a flat cap on the support plane. Its shape is assumed, since no camera measured it.

The voxel size is $1/256$ of the diagonal of the object's box, measured from the 2D splat's masked depth in the cell below, and the truncation distance is ten voxels. A finer grid follows the depth noise more closely. A coarser one rounds off edges. The truncation matters for a depth map that lands too deep on part of a face. The voxels between the true face and that depth look empty to that camera, while every other camera sees them more than $\mu$ behind its own surface and says nothing, so the one bad camera digs a dent into the wall. A wider band lets the cameras that see the face correctly outvote it, at the cost of rounding parts of the object thinner than about $2\mu$.

Before fusing, each depth map is also checked against the cameras two frames before and after it. A pixel whose 3D point lies more than five voxels behind the surface that every such neighbor sees at that spot is dropped, because the neighbors see empty space in front of it. The check cannot help when neighboring views share the same depth error, as they can on a face that only a few frames see. It finds nothing to drop, and where the shared error is deeper than the band, it stays in the mesh as a dent. On the reference object such dents are 2 to 5 mm deep and sit mostly on one of the two narrow faces.

Write `fuse_tsdf`. It renders the 2D splat's depth in every training view, builds the field on a uniform grid, closes it, and returns the largest connected piece of the extracted mesh. That piece must be watertight: every edge shared by exactly two triangles and every vertex surrounded by a single fan of them. Stage 5 computes the volume and inertia from it, and both need a closed surface.

```python theme={null}
box_points = []
for v in train_views:
    depth, alpha, _ = render_depth(splats["2dgs"], v)
    box_points.append(backproject(depth, v["K"], v["viewmat"], v["mask"] & (alpha > 0.5))[::20])
box_lo, box_hi = object_box(np.concatenate(box_points))
diagonal = float(np.linalg.norm(box_hi - box_lo))
VOXEL_M = diagonal / 256
print(f"object box {1000 * (box_hi - box_lo).round(4)} mm, diagonal {1000 * diagonal:.1f} mm: "
      f"voxel {1000 * VOXEL_M:.2f} mm")
```

```output theme={null}
object box [ 47.6  67.5 117.9] mm, diagonal 144.0 mm: voxel 0.56 mm
```

```python theme={null}
import cv2
import numpy as np
import open3d as o3d
from scipy import ndimage

MAX_UNSEEN_AT_BOX = 0.05  # share of never-seen voxels on one face of the object box that means a cut-off


def fuse_tsdf(splat2d, cameras, masks, voxel_m: float, floor_z: float = 0.0,
              info: dict | None = None) -> o3d.geometry.TriangleMesh:
    """Fuse the masked 2DGS depth of every camera into a TSDF and extract a watertight mesh.

    cameras: views with "name", "K" (3 x 3) and "viewmat" (world to camera, 4 x 4).
    masks: {name: boolean mask at the camera's resolution}. floor_z: height of the support plane.
    info: if given, a dict that receives "unseen_outside", "unseen_at_box" and "cut_off" (step 3
    and 4); the check reads "cut_off".

    1. For each camera, render_depth(splat2d, cam) gives the 2DGS median depth, opacity and
       color. Keep the depth where the mask, eroded by a few pixels, is set and the opacity is
       above 0.5, then drop the pixels that drop_too_deep(depths, kept, cameras, 5 * voxel_m)
       rejects. Where the mask is not set, use a very large depth (say 100 m): that ray
       misses the object, so the TSDF marks everything along it empty.
       Elsewhere set the depth to 0 (no measurement).
    2. Size the grid from the object's box: backproject the kept depths, take object_box of the
       points, lower its bottom to floor_z, and add a margin of a few voxels. Build an
       o3d.pipelines.integration.UniformTSDFVolume with voxel size voxel_m (length = resolution
       * voxel_m, origin at the box corner) and sdf_trunc of about 10 voxels, and integrate one
       RGBDImage per camera (depth_scale=1, a depth_trunc above your large depth) with the
       camera's intrinsics and viewmat as the extrinsic.
    3. Close the surface: read the grid with extract_volume_tsdf() (an N x 2 array of value and
       weight, x-major). Set never-seen voxels (weight 0) to -1 if their center lies within the
       object box (from step 2, bottom lowered to floor_z) grown by 3 voxels on every side, and to
       +1 outside it: the grid is a cube sized by the longest side, so above a long, low object it
       holds space that no camera sees and that is empty. Count the second kind into
       info["unseen_outside"]. Set voxels below floor_z and the grid's border to +1, and give
       every voxel a weight of at least 1. Then fill pockets:
       label the face-connected regions of empty voxels (value > 0) with
       scipy.ndimage.label, and set every empty voxel outside the region that holds the
       border to -1, so the surface is one outer shell. Write the grid back with
       inject_volume_tsdf.
    4. extract_triangle_mesh(). The surface was cut off where no camera saw it if, on any side
       or the top of the grown box (not its bottom), more than MAX_UNSEEN_AT_BOX of the voxels in
       the box's outermost layer were never seen (store the largest share in
       info["unseen_at_box"]), or if near_grid_border(mesh, origin, length, voxel_m) is True. Set
       info["cut_off"] to that, so the check can report it. Keep the largest
       connected component and check it with is_watertight. If that fails, try fill_holes() on
       the tensor mesh
       (o3d.t.geometry.TriangleMesh.from_legacy(mesh).fill_holes().to_legacy()).
    """
    raise NotImplementedError
```

```python theme={null}
t = time.time()
fuse_info = {}
mesh = fuse_tsdf(splats["2dgs"], train_views, {v["name"]: v["mask"] for v in train_views}, VOXEL_M,
                 info=fuse_info)
watertight = is_watertight(mesh)
extent = mesh.get_axis_aligned_bounding_box().get_extent()
volume_ml = mesh_volume(mesh) * 1e6 if watertight else float("nan")
print(f"mesh: {len(mesh.vertices):,} vertices, {len(mesh.triangles):,} triangles in {time.time() - t:.0f} s")
print(f"extent {1000 * extent.round(4)} mm, watertight: {watertight}, enclosed volume {volume_ml:.0f} mL")
```

```output theme={null}
mesh: 131,948 vertices, 263,900 triangles in 17 s
extent [ 47.7  68.4 119.3] mm, watertight: True, enclosed volume 295 mL
```

## Splats and mesh in the viewer

The viewer shows the centers of both splats, 50,000 of each at most, with their base colors, and the mesh, all in the world frame of stage 1. Turn the splats off to see the mesh alone, and look at its bottom: it is the flat cap on the support plane.

```python theme={null}
import open3d as o3d


def splat_points(splat, limit: int = 50_000, seed: int = 0):
    """Up to `limit` splat centers with their base colors, as numpy arrays."""
    means = splat["means"].cpu().numpy()
    rgb = (np.clip(0.5 + 0.28209479 * splat["sh0"][:, 0].cpu().numpy(), 0, 1) * 255).astype(np.uint8)
    keep = np.random.default_rng(seed).permutation(len(means))[:limit]
    return means[keep], rgb[keep]


def log_stage3(splats, mesh):
    for kind, splat in splats.items():
        pts, rgb = splat_points(splat)
        rr.log(f"world/splat_{kind}", rr.Points3D(pts, colors=rgb, radii=0.0005), static=True)
    rr.log("world/mesh", rr.Mesh3D(vertex_positions=np.asarray(mesh.vertices),
                                   triangle_indices=np.asarray(mesh.triangles),
                                   vertex_normals=np.asarray(mesh.vertex_normals),
                                   vertex_colors=(np.asarray(mesh.vertex_colors) * 255).astype(np.uint8)),
           static=True)


log_stage3(splats, mesh)
```

## Check

On any object the mesh must be watertight. The check also casts a fixed grid of rays along $x$, $y$ and $z$ through the mesh's box, kept a tenth of the box away from its faces, and counts how often each ray crosses the surface. A ray through one shell crosses it twice. A larger share of rays with more crossings means pockets or dents in the walls. A ray that runs close to a dented wall passes in and out of it, so even the reference keeps a few percent; on the reference object it fails above a limit set from the reference run, and on your object it is a warning, because an object that is not convex, such as a mug with a handle, crosses more often for real. On the reference object the check also compares the mesh with the OmniObject3D scan. The scan is in millimeters (inferred: the dataset does not state its unit) and the mesh in meters, and stage 1 recorded the transform between the scan and the world frame. Both surfaces are sampled with 200,000 points each, and the symmetric Chamfer distance is

$$
d_{\mathrm{CD}} = \frac{1}{2}\left( \frac{1}{|A|} \sum_{a \in A} \min_{b \in B} \lVert a - b \rVert + \frac{1}{|B|} \sum_{b \in B} \min_{a \in A} \lVert a - b \rVert \right),
$$

where $A$ samples the mesh and $B$ the scan. The first term is the mesh's accuracy and the second the scan's completeness. The check passes when $d_{\mathrm{CD}}$ is at most 1.5 times the value of the reference carton's run, a fixed 1.24 mm.

The Chamfer distance measures the mesh only if the scan sits where the object is. Stage 1 placed the scan by fitting it to the sparse points, and on a long object seen mostly from above that fit can slide the scan along the object's length while every stage 1 check passes. So the check also fits the scan to the mesh directly, with a shift, a rotation and a scale found by iterative closest points, and prints how far the scan's center moves and the Chamfer distance after the fit next to the one before it. On the reference object a move of more than 10 mm fails. The fusion also reports when space that no camera measured reaches a side or the top of the object's box: part of the object was then outside the frame in every view, and the mesh there is filled in rather than measured. That fails on the reference object and is a warning on yours. On your own object there is no scan. A splat below 20 dB on the held-out frames or a mesh smaller than 5 mm or larger than 1 m prints a warning; the first usually means bad poses or masks, the second a scale problem in stage 1.

```python theme={null}
import open3d as o3d
from pathlib import Path


def write_stage3(out: Path, splats: dict, mesh) -> None:
    """Write stage3/ in one step: the two splats and the mesh, in meters."""
    with atomic_dir(out) as d:
        save_ply(splats["3dgs"], d / "splat.ply")
        save_ply(splats["2dgs"], d / "splat2d.ply")
        if not o3d.io.write_triangle_mesh(str(d / "mesh.ply"), mesh):
            raise RuntimeError("mesh.ply could not be written; stage3/ was not written")
```

```python theme={null}
import numpy as np
import open3d as o3d
from pathlib import Path
from scipy.spatial import cKDTree


def chamfer_mm(mesh, scan_path: Path, scan_from_world, n: int = 200_000, seed: int = 0) -> dict:
    """Symmetric Chamfer distance between the mesh (meters) and the scan (millimeters, inferred).

    Both surfaces are sampled uniformly with n points, and the scan is brought into the world
    frame with the inverse of scan_from_world. Accuracy is the mean distance from mesh points to
    the scan, completeness from scan points to the mesh, and the Chamfer distance their mean.
    """
    o3d.utility.random.seed(seed)
    scan = o3d.io.read_triangle_mesh(str(scan_path))
    scan.transform(np.linalg.inv(np.asarray(scan_from_world)))  # millimeters to world meters
    a = np.asarray(mesh.sample_points_uniformly(n).points)
    b = np.asarray(scan.sample_points_uniformly(n).points)
    acc = cKDTree(b).query(a)[0] * 1000
    comp = cKDTree(a).query(b)[0] * 1000
    return {"accuracy_mm": float(acc.mean()), "completeness_mm": float(comp.mean()),
            "chamfer_mm": float((acc.mean() + comp.mean()) / 2)}


def register_scan(mesh, scan_path: Path, scan_from_world, n: int = 200_000, seed: int = 0,
                  max_distance_m: float = 0.02) -> dict:
    """Re-register the scan to the mesh, starting from the stage 1 placement.

    The scan's points are moved so their centroid matches the mesh's, then point-to-point ICP with
    a scale factor fits them to the mesh. Returns offset_mm (how far the scan's center moved from
    its stage 1 placement), scale, rotation_deg, and scan_from_world for the new placement, which
    chamfer_mm accepts.
    """
    o3d.utility.random.seed(seed)
    scan = o3d.io.read_triangle_mesh(str(scan_path))
    world_from_scan = np.linalg.inv(np.asarray(scan_from_world, dtype=float))
    scan.transform(world_from_scan)
    a = mesh.sample_points_uniformly(n)
    b = scan.sample_points_uniformly(n)
    center = np.asarray(b.points).mean(0)
    init = np.eye(4)
    init[:3, 3] = np.asarray(a.points).mean(0) - center
    reg = o3d.pipelines.registration.registration_icp(
        b, a, max_distance_m, init,
        o3d.pipelines.registration.TransformationEstimationPointToPoint(with_scaling=True),
        o3d.pipelines.registration.ICPConvergenceCriteria(max_iteration=200))
    T = np.asarray(reg.transformation)
    scale = float(np.cbrt(np.linalg.det(T[:3, :3])))
    turn = np.degrees(np.arccos(np.clip((np.trace(T[:3, :3] / scale) - 1) / 2, -1, 1)))
    return {"offset_mm": float(np.linalg.norm(T[:3, :3] @ center + T[:3, 3] - center) * 1000),
            "scale": scale, "rotation_deg": float(turn), "fitness": float(reg.fitness),
            "scan_from_world": np.linalg.inv(T @ world_from_scan)}


MAX_CROSSING_FRACTION = 0.12  # the reference mesh's value plus a margin
MAX_SCAN_OFFSET_MM = 10.0  # how far re-registration may move the scan from its stage 1 placement


def stage3_problems(mesh, stats: dict, chamfer: dict | None, reference_chamfer_mm: float | None,
                    crossings: float | None = None, cut_off: bool = False, offset_mm: float | None = None):
    """Return (failures, warnings) for the stage 3 check. chamfer is None off the reference object.

    crossings is crossing_fraction(mesh): above MAX_CROSSING_FRACTION it fails on the reference
    object and warns on yours. cut_off is fuse_tsdf's flag for a surface cut off where no camera
    saw it: it fails on the reference object and warns on yours. offset_mm is register_scan's
    offset of the scan from its stage 1 placement (reference object only): above
    MAX_SCAN_OFFSET_MM it fails.
    """
    failures, warnings = [], []
    if len(mesh.triangles) == 0:
        failures.append("the mesh is empty: check the masks and the 2DGS depth")
        return failures, warnings
    if not is_watertight(mesh):
        failures.append("the mesh is not watertight: stage 5 needs a closed surface for its volume and inertia")
    for kind in ("3dgs", "2dgs"):
        if stats[kind]["heldout_psnr"] < 20:
            warnings.append(f"the {kind.upper()} splat reaches only {stats[kind]['heldout_psnr']:.1f} dB on the "
                            "held-out frames: check the poses and the masks")
    extent = mesh.get_axis_aligned_bounding_box().get_extent()
    if extent.max() > 1.0 or extent.min() < 0.005:
        warnings.append(f"the mesh spans {1000 * extent.round(4)} mm: check the stage 1 scale")
    if crossings is not None and crossings > MAX_CROSSING_FRACTION:
        message = (f"{100 * crossings:.1f}% of the rays through the mesh cross its surface more than twice "
                   f"(at most {100 * MAX_CROSSING_FRACTION:.1f}%): the mesh has pockets or folded walls, which "
                   "stage 5 would count in the volume and inertia")
        (failures if chamfer is not None else warnings).append(message)
    if cut_off:
        message = ("the surface is cut off where no camera saw it (the edge of the object's box or of the grid): "
                   "part of the object may be outside the frame in every view, and the mesh there is filled in, "
                   "not measured")
        (failures if chamfer is not None else warnings).append(message)
    if offset_mm is not None and offset_mm > MAX_SCAN_OFFSET_MM:
        failures.append(f"re-registering the scan to the mesh moves it {offset_mm:.1f} mm from its stage 1 placement "
                        f"(at most {MAX_SCAN_OFFSET_MM:.0f} mm): stage 1 placed the world wrongly, and the Chamfer "
                        "distance measures that, not the mesh")
    if chamfer is not None and reference_chamfer_mm is not None:
        limit = 1.5 * reference_chamfer_mm
        if chamfer["chamfer_mm"] > limit:
            failures.append(f"Chamfer distance to the scan {chamfer['chamfer_mm']:.2f} mm is above "
                            f"{limit:.2f} mm (1.5 times the reference run's {reference_chamfer_mm:.2f} mm)")
    return failures, warnings
```

```python theme={null}
REFERENCE_CHAMFER_MM = 1.24  # the solutions run on the reference carton, a fixed number
chamfer = placed = None
if REFERENCE:
    scale = json.loads((STAGE1 / "scale.json").read_text())
    SCAN = PROJECT / "ground_truth" / "scan" / "Scan.obj"
    chamfer = chamfer_mm(mesh, SCAN, scale["scan_from_world"])
    placed = register_scan(mesh, SCAN, scale["scan_from_world"])
    chamfer_placed = chamfer_mm(mesh, SCAN, placed["scan_from_world"])
    print(f"against the scan (millimeters, unit inferred): accuracy {chamfer['accuracy_mm']:.2f} mm, "
          f"completeness {chamfer['completeness_mm']:.2f} mm, Chamfer {chamfer['chamfer_mm']:.2f} mm "
          f"(limit {1.5 * REFERENCE_CHAMFER_MM:.2f} mm)")
    print(f"scan re-registered to the mesh: its center moves {placed['offset_mm']:.1f} mm from the stage 1 placement "
          f"(at most {MAX_SCAN_OFFSET_MM:.0f} mm), scale x{placed['scale']:.3f}, turn {placed['rotation_deg']:.2f} "
          f"degrees; Chamfer {chamfer_placed['chamfer_mm']:.2f} mm after, {chamfer['chamfer_mm']:.2f} mm with the "
          "stage 1 placement")
crossings = crossing_fraction(mesh)
print(f"{len(mesh.triangles):,} triangles, watertight: {is_watertight(mesh)}, rays crossing more than twice: "
      f"{100 * crossings:.1f}%; stage time {(time.time() - T0) / 60:.1f} min")
failures, warnings = stage3_problems(mesh, stats, chamfer, REFERENCE_CHAMFER_MM if REFERENCE else None, crossings,
                                     cut_off=fuse_info.get("cut_off", False),
                                     offset_mm=placed["offset_mm"] if placed else None)
for w in warnings:
    print("warning:", w)
assert not failures, "; ".join(failures)  # nothing is written when a gate fails
write_stage3(PROJECT / "stage3", {k: splats[k] for k in ("3dgs", "2dgs")}, mesh)
if REFERENCE:
    (PROJECT / "stage3" / ".reference").touch()  # stage 4 then knows it runs on the reference
print("passed")
```

```output theme={null}
against the scan (millimeters, unit inferred): accuracy 1.29 mm, completeness 1.21 mm, Chamfer 1.25 mm (limit 1.86 mm)
scan re-registered to the mesh: its center moves 2.3 mm from the stage 1 placement (at most 10 mm), scale x1.013, turn 1.88 degrees; Chamfer 1.04 mm after, 1.25 mm with the stage 1 placement
```

```output theme={null}
263,900 triangles, watertight: True, rays crossing more than twice: 4.3%; stage time 11.4 min
```

```output theme={null}
passed
```

### Question 1. Depth from a blob and from a disk

Explain why the expected depth of a 3D Gaussian splat need not lie on any surface, and name the change in how a 2D Gaussian is rendered that makes its depth a point on a surface. Point to a region of your depth figure where the two normal maps differ.

*Write your answer here.*

### Question 2. Fusing depth into a signed distance field

Explain what one depth map contributes to the field $F(x)$, including the rays that fall outside the mask, and why the truncation distance $\mu$ is needed. Point to the printed voxel size and the extent of your mesh.

*Write your answer here.*

### Question 3. Closing the field

After fusion, some voxels were never measured by any camera. Explain how the fusion decides whether such a voxel is inside or outside the object, and describe an object for which that rule gives the wrong shape. Point to the part of your mesh in the viewer whose shape was assumed rather than measured.

*Write your answer here.*

## References

* Kerbl, B., Kopanas, G., Leimkuehler, T., Drettakis, G. (2023). *3D Gaussian splatting for real-time radiance Field rendering*. ACM Trans. Graph.
* Ye, V., Li, R., Kerr, J., Turkulainen, M., Yi, B., et al. (2024). *gsplat: An Open-Source Library for Gaussian Splatting*.

***

Found an error on this page? [Report an issue](https://github.com/aegean-ai/eaia-feedback/issues/new?template=page-issue.yml\&page=/aiml-common/projects/cv/real-to-sim-twin/03-geometry/03-geometry).


This documentation is built and hosted on [Mintlify](https://mintlify.com), a developer documentation platform.