Point Clouds and 3D Modelling (NPM3D)

Contents

• Introduction and overview
• 0. Review: geometry, linear algebra and probability for 3D data
• 1. 3D perception systems
• 2. Processing and operators
• 3. Registration
• 4. Segmentation of point clouds
• 5. Modelling by primitives
• 6. Reconstruction of curves and surfaces
• 7. Rendering point clouds and meshes
• 8. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with an M1 in mathematics, computer science or physics who is fluent in linear algebra and Python, and who wants to turn raw 3D scans into measurements, maps and models.

Empirical first: put a 16-beam rotating laser scanner 2 m above a flat road. Every range it returns is accurate to a few centimetres, yet the seven rings it draws on the ground within 40 m lie between 1.2 m and 15.3 m apart (Code 1.7). Draw those points one pixel each and 86 % of the object you look at is a hole, with the back of the object showing through (Code 7.7). Point clouds are precise and sparse, unordered and uneven: this course builds the operators that cope with both facts.

The course follows the whole chain from the sensor to the model. It explains how 3D points are acquired (ch. 1), which local operators extract normals, curvature and descriptors (ch. 2), how two scans are aligned (ch. 3), how a cloud is cut into meaningful parts (ch. 4), fitted with geometric primitives (ch. 5), turned into a surface (ch. 6) and drawn on a screen (ch. 7).

Its image-side companion is 3D Computer Vision (camera geometry, stereo, multi-view reconstruction). Meshes and learning on shapes are developed in Geometry Processing and Geometric Deep Learning, and LiDAR-based localisation reappears in Robotics.

What you will be able to do

• Explain how triangulation scanners, time-of-flight sensors and mobile mapping systems produce points, and predict their sampling density and their errors.
• Build neighbourhoods with kd-trees or voxel grids, subsample a cloud, and estimate normals and curvature by local PCA at a justified scale.
• Compute covariance and histogram descriptors and use them to classify points.
• Derive and implement the closed-form rigid alignment of two point sets, the ICP algorithm and its point-to-plane variant, and prove their basic properties.
• Make registration robust (trimming, M-estimators) and global (feature matching with RANSAC), and size RANSAC with the right number of iterations.
• Segment a cloud by region growing, graph cuts, random forests on multi-scale features or a point-based neural network (PointNet, KPConv).
• Detect planes, spheres and cylinders with RANSAC and the Hough transform.
• Reconstruct curves and surfaces from samples with Delaunay-based methods, signed distances, marching cubes and Poisson reconstruction.
• Render point clouds and meshes with a z-buffer, splats and level-of-detail hierarchies, and explain holes, aliasing and show-through.
• Read a recent research paper on point clouds, reimplement its algorithm and evaluate it, as the course project asks.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com sets the objective of giving a panorama of the concepts and techniques for the acquisition, processing and visualisation of 3D point clouds, together with their mathematical and algorithmic foundations.

It lists seven themes (“Thèmes abordés”), matched below. The Université Paris-Saclay course detail gives the volume (16 h of lectures, 16 h of labs) and the nature of the assessment but no programme items. For the summary and schedule of the sessions, master-mva.com points to the NPM3D page on caor.minesparis.psl.eu, which these notes do not reproduce.

Beyond the syllabus, this page adds a review chapter (ch. 0), a research chapter for self-study on learned registration, neural implicit surfaces and Gaussian splatting (ch. 8), and a written-exam proposal.

Where this course fits in the programme

Official prerequisites: geometry, linear algebra and probability; programming in C, C++, Python or Matlab is useful for the project. The course page on master-mva.com lists Jean-Emmanuel Deschaud and François Goulette (MINES ParisTech) and Tamy Boubekeur (Telecom ParisTech) as instructors. It also says that external auditors may attend, without any grade, assessment or certificate.

• Refreshers — Linear Algebra, Analysis and Optimization for the SVD, the spectral theorem and least squares; Probability and Statistics for noise models and the binomial arguments behind RANSAC; Scientific Python, NumPy and PyTorch for the labs.
• Semester 1 courses that prepare it — 3D Computer Vision (projective geometry, stereo), Geometric Data Analysis (distances between shapes), Geometry Processing and Geometric Deep Learning (meshes, learning on shapes), Deep Learning (the networks of ch. 4) and Robotics (localisation).
• Related courses of semester 2 — Geometry and Shape Spaces, Deformable Models and Geodesic Methods for Image Analysis, Remote Sensing Data (airborne LiDAR) and Graphs in Machine Learning (neighbourhood graphs).
• Afterwards — the Research Internship, for which the paper-study project of this course is a rehearsal.

Course outline

• Review: geometry, linear algebra and probability for 3D data — rigid motions, rotations, the spectral theorem, least squares and the probability tools used later.
• 3D perception systems — how triangulation, time of flight and mobile mapping produce points, and what their errors look like.
• Processing and operators — neighbourhoods, subsampling, normals, curvature, descriptors and denoising.
• Registration — aligning scans: closed form, ICP, point-to-plane, robust and global methods.
• Segmentation of point clouds — from region growing and graph cuts to random forests and point-based neural networks.
• Modelling by primitives — planes, spheres and cylinders found by RANSAC and the Hough transform.
• Reconstruction of curves and surfaces — Delaunay-based meshing, implicit functions, marching cubes and Poisson reconstruction.
• Rendering point clouds and meshes — z-buffer, shading, surfels, splatting and level of detail.
• Research directions and open problems — learned registration and odometry, large-scale segmentation, neural surfaces and Gaussian splatting (self-study).

Session plan

32 h in class = 8 sessions of 4 h.

The volume is the Université Paris-Saclay programme table (16 h of lectures and 16 h of labs). The split into eight sessions of 2 h of lecture followed by 2 h of lab is a proposal; the official schedule is on the course page linked from master-mva.com.

Level 1 = foundations, 2 = core methods, 3 = advanced methods close to research; the level never decreases from one session to the next.

Self-study, outside the 32 h (about 4 h in total): ch. 8 Research directions and open problems — learned registration and LiDAR odometry, large-scale segmentation, neural implicit surfaces, Gaussian splatting (≈ 4 h, with Ex. 19).

How to study this module

• Read ch. 0 first even if the material looks familiar: the SVD argument and the Rayleigh quotient return in chapters 2, 3, 5 and 6.
• After each session, run the code of the chapter and change one parameter (the neighbourhood radius, the initial rotation, the RANSAC threshold) to see what breaks.
• Do the exercises listed in the session plan before the next session; the ★★★ ones are good starting points for the project.
• Python environment: NumPy, SciPy (scipy.spatial for kd-trees, Delaunay and Voronoi), scikit-learn and PyTorch on CPU; every code block runs in a few seconds on a laptop with synthetic data.
• For visual checks, export clouds as PLY or text files and open them in a point cloud viewer such as CloudCompare. Real data come as PLY, LAS/LAZ (airborne and mobile LiDAR) or E57 files; the open-source libraries PCL (Rusu & Cousins, 2011) and Open3D (Zhou, Park & Koltun, 2018) implement most operators of this course.
• Start choosing a paper for the project around session 4: the list of solo subjects gives the scale to aim at.

0. Review: geometry, linear algebra and probability for 3D data

Every later chapter moves points with rigid motions, extracts directions from small covariance matrices, solves least-squares problems and counts the chances of random draws. This review states these tools once, with the proofs that later chapters rely on.

For a unit vector k, let [k]_\times be the matrix of the cross product, so that [k]_\times v = k \times v: [k]_\times = \begin{pmatrix} 0 & -k_3 & k_2 \\ k_3 & 0 & -k_1 \\ -k_2 & k_1 & 0 \end{pmatrix}.

Small rotations are the workhorse of iterative registration (ch. 3). For a rotation vector \omega = \theta k with small \theta, Rodrigues’ formula gives R \approx I + [\omega]_\times , \qquad R x \approx x + \omega \times x .

This theorem is behind the normal of a plane fitted to points (ch. 2), the plane primitive (ch. 5) and the covariance descriptors (ch. 2).

import numpy as np

def rodrigues(axis, theta):
    """Rotation of angle theta about the unit axis k (Rodrigues' formula)."""
    k = np.asarray(axis, float) / np.linalg.norm(axis)
    K = np.array([[0, -k[2], k[1]], [k[2], 0, -k[0]], [-k[1], k[0], 0]])
    return np.eye(3) + np.sin(theta) * K + (1 - np.cos(theta)) * K @ K

R = rodrigues([1, 1, 0], np.pi / 3)
print(np.round(R, 4))
print("orthogonal:", np.allclose(R.T @ R, np.eye(3)), "| det =", round(np.linalg.det(R), 6))
print("angle from the trace (deg):", round(np.degrees(np.arccos((np.trace(R) - 1) / 2)), 4))

# Nearest rotation to a noisy matrix: projection onto SO(3) by the SVD
rng = np.random.default_rng(0)
M = R + 0.05 * rng.standard_normal((3, 3))
U, s, Vt = np.linalg.svd(M)
S = np.diag([1, 1, np.sign(np.linalg.det(U @ Vt))])
R_hat = U @ S @ Vt
print("||M - R||_F =", round(np.linalg.norm(M - R), 4),
      "| ||R_hat - R||_F =", round(np.linalg.norm(R_hat - R), 4))

1. 3D perception systems

Before processing points we must know where they come from. This chapter describes the main families of 3D sensors, derives the accuracy and the sampling density of each, and lists the artefacts that later chapters have to survive.

Sensors fall into two physical principles. Triangulation measures an angle seen from two places a known distance apart. Time of flight measures how long light takes to go to the surface and back.

Both can be passive (cameras only, as in photogrammetry and multi-view stereo) or active (the sensor emits a laser stripe, a projected pattern or light pulses). Active sensors work on textureless surfaces and in the dark, which is why they dominate 3D scanning.

Time of flight gives the range along a known direction. A laser scanner sweeps that direction with rotating mirrors or a rotating head and returns, for every shot, a range r, an azimuth \theta and an elevation \varphi.

The range image is a useful by-product: neighbours in the image are usually neighbours in space, which gives cheap normals and segmentation on raw scans.

import numpy as np
rng = np.random.default_rng(0)

# A 16-beam rotating scanner, 2 m above a flat ground, in front of a wall and a sphere.
h = 2.0
elev = np.radians(np.linspace(-15, 15, 16))           # beam elevations
azim = np.radians(np.arange(0, 360, 0.2))             # azimuth step 0.2 degrees
E, A = np.meshgrid(elev, azim, indexing="ij")
D = np.stack([np.cos(E) * np.cos(A), np.cos(E) * np.sin(A), np.sin(E)], -1).reshape(-1, 3)
o = np.array([0.0, 0.0, h])

def hit_plane(n, c):          # plane n.x = c ; returns ray parameter t (inf if missed)
    den = D @ n
    with np.errstate(divide="ignore", invalid="ignore"):
        t = (c - o @ n) / den
    return np.where((den != 0) & (t > 0), t, np.inf)

def hit_sphere(center, rad):
    oc = o - center
    b = D @ oc
    disc = b**2 - (oc @ oc - rad**2)
    t = -b - np.sqrt(np.maximum(disc, 0))
    return np.where((disc > 0) & (t > 0), t, np.inf)

t = np.minimum.reduce([hit_plane(np.array([0, 0, 1.0]), 0.0),     # ground z = 0
                       hit_plane(np.array([1.0, 0, 0]), 12.0),     # wall x = 12
                       hit_sphere(np.array([5.0, 3.0, 1.0]), 1.0)])
keep = t < 100.0                                                   # maximal range 100 m
r = t[keep] + 0.02 * rng.standard_normal(keep.sum())               # range noise, 2 cm
P = o + r[:, None] * D[keep]
print("rays:", len(D), "| returns:", len(P))

# Ground rings: beam k hits the ground at horizontal distance h / tan(-elevation)
down = -elev[elev < 0]
rings = h / np.tan(down)
print("ring radii (m):", np.round(np.sort(rings), 2))
near = np.sort(rings)[np.sort(rings) < 40]
print("gaps between consecutive rings within 40 m (m):", np.round(np.diff(near), 2))

Static scans cover one viewpoint. Mobile mapping puts the scanner on a vehicle, a backpack or a drone and places every shot in a world frame, using a satellite receiver (GNSS) and an inertial measurement unit (IMU) whose fusion gives the pose of the platform at each instant.

With a heading error of 0.01° (1.75 \times 10^{-4} rad), a point at 50 m moves by 8.7 mm. At 10 m/s, a synchronisation error of 1 ms moves every point by 1 cm.

These numbers explain why mobile maps are refined by registering overlapping passes (ch. 3) rather than trusted from GNSS and IMU alone, which can drift by decimetres when satellites are masked by buildings.

2. Processing and operators

A raw cloud is a list of coordinates. To reason about surfaces we need, for every point, its neighbours and the local shape they describe. This chapter builds neighbourhoods efficiently, thins clouds, and computes normals, curvature and descriptors, with attention to the choice of scale.

Neighbourhoods and spatial data structures

Brute force costs O(n) per query, hence O(n^2) for all points: hopeless for clouds of 10^7 points. Spatial data structures reduce this cost.

Two other structures are common. An octree cuts space (not the points) into eight cubes recursively; it is the natural structure for multi-resolution processing and rendering (ch. 7). A voxel hash maps each point to the integer cell \lfloor x / s \rfloor; a ball query of radius r \le s then reads the 27 cells around the query.

Subsampling

Clouds are often far denser than needed, and denser near the sensor than far away. Subsampling reduces the cost and evens out the density.

FPS is used inside PointNet++ (ch. 4) and has a clean guarantee. The k-center problem asks for k centres c_1, \dots, c_k minimising the covering radius \rho(C) = \max_{x \in \mathcal{X}} \min_{j} \|x - c_j\| .

Normals by local PCA

Most operators need the tangent plane of the underlying surface at each point. It is estimated by fitting a plane to the neighbourhood.

Let \mathcal{N} be a neighbourhood of m points, with centroid and covariance \bar x = \frac{1}{m} \sum_{x_i \in \mathcal{N}} x_i, \qquad C = \frac{1}{m} \sum_{x_i \in \mathcal{N}} (x_i - \bar x)(x_i - \bar x)^\top .

The eigenvector is defined up to sign. A consistent orientation is needed for shading (ch. 7) and for Poisson reconstruction (ch. 6).

Curvature and the choice of scale

The other two eigenvalues measure how far the neighbourhood is from a plane.

The next computation makes that dependence explicit and gives a rule for choosing the radius.

Small radii are dominated by noise, large radii by curvature. For \sigma = 1 cm and a radius of curvature of 0.5 m (\kappa = 2 m⁻¹), the rule gives r^\star \approx 2.63 \times \sqrt{0.005} \approx 0.19\ \text{m}.

Mitra and Nguyen (2003) derive rigorous bounds of the same nature for the error of PCA normals, with density as a third factor.

Descriptors

The eigenvalues describe the dimensionality of the neighbourhood: one dominant direction for a pole or a cable, two for a wall, three for foliage.

Demantké and coauthors (2011) choose, for each point, the radius at which the entropy of dimensionality features of this type (built from the square roots of the eigenvalues) is minimal: the scale at which the neighbourhood looks most clearly 1D, 2D or 3D.

Covariance features ignore how normals vary inside the neighbourhood. Histogram descriptors capture this, and are used to match points between scans (ch. 3).

Outliers and noise

import numpy as np
from scipy.spatial import cKDTree

rng = np.random.default_rng(0)
# A toy street: ground (label 0), a thin pole (1), a bush (2)
g = np.c_[rng.uniform(-3, 3, (3000, 2)), 0.01 * rng.standard_normal(3000)]
th, z = rng.uniform(0, 2 * np.pi, 800), rng.uniform(0, 3, 800)
p = np.c_[0.08 * np.cos(th) + 1.0, 0.08 * np.sin(th) - 1.0, z]
b = rng.normal(0, 1, (1200, 3)); b /= np.linalg.norm(b, axis=1, keepdims=True)
b = b * rng.uniform(0, 0.6, (1200, 1)) ** (1 / 3) + [-1.0, 1.0, 0.7]
X = np.r_[g, p, b]
y = np.r_[np.zeros(3000), np.ones(800), 2 * np.ones(1200)].astype(int)
tree = cKDTree(X)

def pca_features(neigh_lists):
    """Linearity, planarity, sphericity and normal from each local covariance."""
    F, N = np.zeros((len(X), 3)), np.zeros((len(X), 3))
    for i, nb in enumerate(neigh_lists):
        Q = X[nb] - X[nb].mean(0)
        lam, vec = np.linalg.eigh(Q.T @ Q / len(nb))      # ascending order
        l3, l2, l1 = np.maximum(lam, 1e-12)
        F[i] = [(l1 - l2) / l1, (l2 - l3) / l1, l3 / l1]
        N[i] = vec[:, 0]                                  # smallest eigenvector
    return F, N

for name, nbs in [("k = 20", tree.query(X, k=20)[1]),
                  ("r = 0.5 m", tree.query_ball_point(X, r=0.5))]:
    F, N = pca_features(nbs)
    print(name)
    for c, cls in enumerate(["ground", "pole", "bush"]):
        print(f"  {cls:6s} lin/plan/sph =", np.round(F[y == c].mean(0), 2))
    err = np.degrees(np.arccos(np.abs(N[y == 0] @ [0, 0, 1])))
    print("  ground normals, median angular error:", round(np.median(err), 2), "deg")

# Voxel-grid subsampling: one barycentre per occupied 0.2 m cell
keys = np.floor(X / 0.2).astype(int)
_, inv, cnt = np.unique(keys, axis=0, return_inverse=True, return_counts=True)
sub = np.zeros((len(cnt), 3)); np.add.at(sub, inv.ravel(), X); sub /= cnt[:, None]
print("points:", len(X), "-> after the 0.2 m grid:", len(sub))

3. Registration

A scanner sees a scene from one place. To model a building, a street or an object we must bring many scans into one frame, and a moving robot must estimate its motion from successive scans. This chapter solves the alignment problem: first with known correspondences in closed form, then by iterating closest points, then robustly and globally.

Alignment with known correspondences

Write \bar p, \bar q for the centroids, p_i' = p_i - \bar p and q_i' = q_i - \bar q for the centred points, and H = \sum_i p_i' \, q_i'^\top \in \mathbb{R}^{3 \times 3}, \qquad H = U \Sigma V^\top \ \text{(SVD)}.

The reflection case is not exotic: it happens when the points are nearly coplanar (\sigma_3 \approx 0) or very noisy. Forgetting the correction then returns a mirror image.

Iterative closest point

Without correspondences, the natural idea is to guess them: match each source point to its closest target point, align, and repeat.

The theorem guarantees convergence, not correctness. ICP converges to a local minimum, and its basin of attraction is limited: a rotation of a few tens of degrees or a small overlap is enough to lock it into a wrong alignment. ICP is therefore a fine registration method, started from a coarse estimate (GNSS/IMU, odometry, or the global methods below).

Point-to-plane ICP

Matching to the closest point penalises sliding along the surface, although two samplings of the same plane never coincide. Chen and Medioni (1992) penalise only the distance to the tangent plane at the target.

The 6 \times 6 matrix \sum_i a_i a_i^\top shows when registration is ill-posed. On a single plane all n_i are equal: translations within the plane and rotations about its normal leave the cost unchanged, and the matrix has rank 3. In a straight corridor the translation along the axis is unobservable.

LiDAR odometry monitors the eigenvalues of this matrix to detect such degenerate places.

import numpy as np
from scipy.spatial import cKDTree
from scipy.spatial.transform import Rotation

rng = np.random.default_rng(1)

def surface(n):                       # a bumpy, non-symmetric patch with its normals
    u = rng.uniform(-1, 1, (n, 2))
    z = 0.3 * np.sin(2 * u[:, 0]) * np.cos(3 * u[:, 1]) + 0.2 * u[:, 0] ** 2
    dzdx = 0.6 * np.cos(2 * u[:, 0]) * np.cos(3 * u[:, 1]) + 0.4 * u[:, 0]
    dzdy = -0.9 * np.sin(2 * u[:, 0]) * np.sin(3 * u[:, 1])
    nrm = np.c_[-dzdx, -dzdy, np.ones(n)]
    return np.c_[u, z], nrm / np.linalg.norm(nrm, axis=1, keepdims=True)

def kabsch(P, Q):
    """R, t minimising sum ||R p_i + t - q_i||^2 (Theorem 3.2)."""
    mp, mq = P.mean(0), Q.mean(0)
    U, _, Vt = np.linalg.svd((P - mp).T @ (Q - mq))
    S = np.diag([1, 1, np.sign(np.linalg.det(Vt.T @ U.T))])
    R = Vt.T @ S @ U.T
    return R, mq - R @ mp

Q, NQ = surface(4000)                                   # reference ("model") scan
R0 = Rotation.from_rotvec(np.radians(20) * np.array([0.3, 0.2, 0.93])).as_matrix()
P = (surface(3000)[0] @ R0.T) + [0.15, -0.10, 0.05]     # new scan, other samples
P += 0.005 * rng.standard_normal(P.shape)
tree = cKDTree(Q)

def icp(P, plane, iters=30):
    R, t, rms = np.eye(3), np.zeros(3), []
    for _ in range(iters):
        X = P @ R.T + t
        d, j = tree.query(X)
        rms.append(np.sqrt(np.mean(d ** 2)))
        if plane:                                       # linearised point-to-plane step
            n = NQ[j]
            A = np.c_[np.cross(X, n), n]
            b = np.einsum("ij,ij->i", Q[j] - X, n)
            x = np.linalg.lstsq(A, b, rcond=None)[0]
            dR = Rotation.from_rotvec(x[:3]).as_matrix()
            R, t = dR @ R, dR @ t + x[3:]
        else:                                           # point-to-point step
            dR, dt = kabsch(X, Q[j])
            R, t = dR @ R, dR @ t + dt
    return R, rms

for plane in (False, True):
    R, rms = icp(P, plane)
    ang = np.degrees(Rotation.from_matrix(R @ R0).magnitude())
    print("point-to-plane" if plane else "point-to-point",
          "| RMS at iterations 1, 5, 10, 30:", " ".join(f"{rms[i]:.4f}" for i in (0, 4, 9, 29)),
          f"| residual rotation error: {ang:.3f} deg")

Robust registration

Real scans contain points with no counterpart: parts seen from one station only, moving cars, noise. Their closest-point residuals are large and pull the least-squares solution away.

For the Cauchy function the weight is 1/(1 + r^2/c^2): a pair with residual 3c counts ten times less than a perfect pair. Each robust ICP iteration is therefore a weighted Theorem 3.2 or a weighted Proposition 3.7.

Global registration

When no initial guess is available, correspondences are found by comparing descriptors instead of positions.

Most descriptor matches are wrong, so the number of draws matters.

For p = 0.99 and s = 3 the required number of draws explodes as w decreases:

Since w is unknown, implementations start with a large N and decrease it each time a better consensus raises the estimate of w (the adaptive rule used in Code 5.10).

4. Segmentation of point clouds

A registered cloud is still a single list of points. Segmentation cuts it into meaningful parts: smooth surface patches, objects, or semantic classes such as ground, building, vegetation and car. This chapter goes from geometric criteria to learned classifiers, and ends with neural networks that consume raw points.

Geometric segmentation

Region growing returns smooth patches: planes of a façade, a road surface, the cylinder of a pole. A simpler variant, Euclidean clustering, keeps the connected components of the graph that joins points closer than r. After the ground has been removed it separates isolated objects such as cars and pedestrians.

Both methods decide locally. A global formulation instead balances how well each point fits a label against how coherent the labels of neighbours are.

Classification with hand-crafted features

The classical supervised pipeline computes, for every point, the features of Definition 2.9 at several scales and trains a standard classifier on labelled points.

import numpy as np
from scipy.spatial import cKDTree
from sklearn.ensemble import RandomForestClassifier

def scene(seed):
    """Toy street with labels: ground 0, pole 1, bush 2, facade 3."""
    rng = np.random.default_rng(seed)
    g = np.c_[rng.uniform(-3, 3, (3000, 2)), 0.01 * rng.standard_normal(3000)]
    th, z = rng.uniform(0, 2 * np.pi, 600), rng.uniform(0, 3, 600)
    c = rng.uniform(-2, 2, 2)
    p = np.c_[0.08 * np.cos(th) + c[0], 0.08 * np.sin(th) + c[1], z]
    b = rng.normal(0, 1, (1200, 3)); b /= np.linalg.norm(b, axis=1, keepdims=True)
    b = b * rng.uniform(0, 0.6, (1200, 1)) ** (1 / 3) + [*rng.uniform(-2, 2, 2), 0.7]
    f = np.c_[np.full(1500, 3.0), rng.uniform(-3, 3, 1500), rng.uniform(0, 4, 1500)]
    f[:, 0] += 0.01 * rng.standard_normal(1500)
    y = np.repeat([0, 1, 2, 3], [3000, 600, 1200, 1500])
    return np.r_[g, p, b, f], y

def features(X, radii=(0.25, 0.5, 1.0)):
    """Per radius: linearity, planarity, sphericity, verticality; then the height."""
    tree, F = cKDTree(X), []
    for r in radii:
        out = np.zeros((len(X), 4))
        for i, nb in enumerate(tree.query_ball_point(X, r)):
            Q = X[nb] - X[nb].mean(0)
            lam, vec = np.linalg.eigh(Q.T @ Q / len(nb))
            l3, l2, l1 = np.maximum(lam, 1e-12)
            out[i] = [(l1 - l2) / l1, (l2 - l3) / l1, l3 / l1, 1 - abs(vec[2, 0])]
        F.append(out)
    return np.hstack(F + [X[:, 2:3]])

(Xtr, ytr), (Xte, yte) = scene(0), scene(1)       # train on one scan, test on another
Ftr, Fte = features(Xtr), features(Xte)
for name, cols in [("one scale (0.25 m)", slice(0, 4)),
                   ("three scales", slice(0, 12)),
                   ("three scales + height", slice(0, 13))]:
    clf = RandomForestClassifier(n_estimators=100, random_state=0).fit(Ftr[:, cols], ytr)
    pred = clf.predict(Fte[:, cols])
    iou = [np.sum((pred == c) & (yte == c)) / np.sum((pred == c) | (yte == c)) for c in range(4)]
    print(f"{name:22s} IoU ground/pole/bush/facade = {np.round(iou, 3)}  mIoU = {np.mean(iou):.3f}")

Learning directly on points

Hand-crafted features fix the description in advance. Deep networks learn it, but a cloud is neither a grid (no convolution) nor a sequence (no order). The architecture must be insensitive to the order of the points.

The maximum is symmetric, so f is permutation invariant by construction. The next theorem shows that nothing essential is lost.

Zaheer and coauthors (2017, “Deep Sets”) prove a sum-decomposition \rho(\sum_i \phi(x_i)) of the same kind. The max version has an extra property: robustness to missing points.

import numpy as np, torch, torch.nn as nn
torch.set_num_threads(1); torch.manual_seed(0); rng = np.random.default_rng(0)

def sample(kind, n=128):
    """n noisy points on a sphere, cube, cylinder or cone, randomly turned and scaled."""
    if kind == 0:
        x = rng.normal(size=(n, 3)); x /= np.linalg.norm(x, axis=1, keepdims=True)
    elif kind == 1:
        x = rng.uniform(-1, 1, (n, 3)); f = rng.integers(0, 3, n)
        x[np.arange(n), f] = rng.choice([-1.0, 1.0], n)
    else:
        th, z = rng.uniform(0, 2 * np.pi, n), rng.uniform(-1, 1, n)
        rad = np.ones(n) if kind == 2 else (1 - z) / 2
        x = np.c_[rad * np.cos(th), rad * np.sin(th), z]
    a = rng.uniform(0, 2 * np.pi)                        # random rotation about z
    Rz = np.array([[np.cos(a), -np.sin(a), 0], [np.sin(a), np.cos(a), 0], [0, 0, 1]])
    return (x @ Rz.T) * rng.uniform(0.7, 1.3) + 0.02 * rng.normal(size=(n, 3))

def dataset(m):
    y = rng.integers(0, 4, m)
    X = np.stack([sample(k) for k in y])
    return torch.tensor(X, dtype=torch.float32), torch.tensor(y)

class PointNet(nn.Module):
    """Shared per-point MLP h, symmetric max-pooling, then a classifier gamma."""
    def __init__(self, k=4):
        super().__init__()
        self.h = nn.Sequential(nn.Linear(3, 32), nn.ReLU(), nn.Linear(32, 64),
                               nn.ReLU(), nn.Linear(64, 128))
        self.gamma = nn.Sequential(nn.ReLU(), nn.Linear(128, 32), nn.ReLU(), nn.Linear(32, k))
    def forward(self, x):                                # x: (batch, n_points, 3)
        return self.gamma(self.h(x).max(dim=1).values)

Xtr, ytr = dataset(600); Xte, yte = dataset(300)
net = PointNet(); opt = torch.optim.Adam(net.parameters(), lr=2e-3)
for epoch in range(20):
    perm = torch.randperm(len(Xtr))
    for i in range(0, len(Xtr), 50):
        b = perm[i:i + 50]
        loss = nn.functional.cross_entropy(net(Xtr[b]), ytr[b])
        opt.zero_grad(); loss.backward(); opt.step()
with torch.no_grad():
    acc = (net(Xte).argmax(1) == yte).float().mean().item()
    same = torch.allclose(net(Xte), net(Xte[:, torch.randperm(128)]), atol=1e-5)
print(f"test accuracy (sphere/cube/cylinder/cone): {acc:.3f}")
print("output unchanged when the points are reordered:", same)

Local structure: hierarchies and point convolutions

Grouping with several radii at once (“multi-scale grouping”) makes the network robust to density changes, the learned counterpart of the multi-scale features of Definition 4.5.

KPConv is a genuine convolution: the same filter is applied at every location, and it depends on relative positions only, so it is translation equivariant.

5. Modelling by primitives

Man-made scenes are mostly planes, cylinders, spheres and cones: walls, pipes, tanks, columns. Describing a scan by a few such primitives compresses it by orders of magnitude and turns points into an editable model. This chapter fits primitives, detects them in clutter with RANSAC and the Hough transform, and chooses between them.

The geometric distance makes sphere and cylinder fitting non-linear. A change of unknowns makes the sphere linear.

Detection with RANSAC

On large scenes a single primitive holds a tiny fraction of the points, and Theorem 3.12 becomes prohibitive. Schnabel, Wahl and Klein (2007) draw the samples locally.

Take 10^6 points, a plane of 1{,}000 points and samples of size 3. Uniform sampling needs about 4.6 \times 10^9 draws for p = 0.99. With cells in which the plane holds half of the points, q = 10^{-3} \times 0.25 and N \approx \frac{4.6}{2.5 \times 10^{-4}} \approx 18{,}400 \text{ draws}.

Efficient RANSAC also estimates scores on random subsets of the cloud, and keeps only the largest connected component of the inliers, so that two separate walls in the same plane become two primitives.

The Hough transform

RANSAC samples the data and tests parameters. The Hough transform does the reverse: every point votes for all the primitives that pass through it, in a discretised parameter space.

With 1° cells, a 10 m façade spreads its votes over 8.7 cm of \rho: finer \rho-cells split the peak, coarser ones merge nearby planes. Cells of 1° also make 360 \times 90 = 32{,}400 angle cells, hence 3 \times 10^{10} votes for 10^6 points. Randomised variants and better accumulator designs (Borrmann et al., 2011) reduce both problems.

Choosing between primitives

import numpy as np
rng = np.random.default_rng(0)

# A room corner: floor z=0, two walls x=0 and y=0, plus 30 % clutter
def patch(axis, n):
    P = rng.uniform(0, 4, (n, 3)); P[:, axis] = 0.01 * rng.standard_normal(n); return P
X = np.r_[patch(2, 2000), patch(0, 1500), patch(1, 1000), rng.uniform(0, 4, (1900, 3))]

def ransac_plane(X, eps=0.03, p=0.99, max_iter=10_000):
    """Plane with the largest consensus set; iteration count updated adaptively."""
    best, N, k = np.zeros(len(X), bool), max_iter, 0
    while k < N:
        a, b, c = X[rng.choice(len(X), 3, replace=False)]
        n = np.cross(b - a, c - a)
        if np.linalg.norm(n) < 1e-9:
            k += 1; continue
        n /= np.linalg.norm(n)
        inl = np.abs((X - a) @ n) < eps
        if inl.sum() > best.sum():
            best = inl
            w = inl.mean()                                    # current inlier ratio
            N = min(max_iter, int(np.ceil(np.log(1 - p) / np.log(1 - w ** 3))))
        k += 1
    Q = X[best] - X[best].mean(0)                             # least-squares refit
    n = np.linalg.svd(Q, full_matrices=False)[2][-1]
    return n, best, k

rest = X.copy()
for found in range(3):                                        # sequential RANSAC
    n, inl, k = ransac_plane(rest)
    print(f"plane {found + 1}: normal {np.round(np.abs(n), 3)}, "
          f"{inl.sum()} inliers, {k} iterations")
    rest = rest[~inl]
print("points left:", len(rest))

# Algebraic sphere fit: |x|^2 = 2 c.x + (r^2 - |c|^2) is linear in (c, d)
S = rng.normal(size=(500, 3)); S /= np.linalg.norm(S, axis=1, keepdims=True)
S = 0.5 * S[S[:, 2] > 0] + [1, 2, 3]                       # upper half only
S += 0.005 * rng.standard_normal(S.shape)
A = np.c_[2 * S, np.ones(len(S))]
sol = np.linalg.lstsq(A, (S ** 2).sum(1), rcond=None)[0]
c = sol[:3]; r = np.sqrt(sol[3] + c @ c)
print("sphere centre", np.round(c, 3), "radius", round(r, 3))

6. Reconstruction of curves and surfaces

Points are samples; most applications need the continuous object they come from: a mesh to render, measure or print, or a function whose zero set is the surface. This chapter asks when samples determine a curve or a surface, and builds reconstructions from Delaunay triangulations and from implicit functions, up to Poisson reconstruction.

Sampling conditions and curve reconstruction

The Voronoi cell of a sample is the set of points closer to it than to any other sample. The Delaunay triangulation is the dual structure: it joins two samples when their cells touch, and no sample lies inside the circumscribed circle of any of its triangles.

import numpy as np
from scipy.spatial import Delaunay, Voronoi

rng = np.random.default_rng(0)
# 200 random samples of the flower r = 1 + 0.3 cos(5 theta)
th = rng.uniform(0, 2 * np.pi, 200)               # random positions on the curve
S = np.c_[(1 + 0.3 * np.cos(5 * th)) * np.cos(th), (1 + 0.3 * np.cos(5 * th)) * np.sin(th)]

def crust(S):
    """2D crust: Delaunay edges between samples that survive adding the Voronoi vertices."""
    V = Voronoi(S).vertices
    T = Delaunay(np.r_[S, V]).simplices
    E = set()
    for tri in T:
        for a, b in [(tri[0], tri[1]), (tri[1], tri[2]), (tri[0], tri[2])]:
            if a < len(S) and b < len(S):
                E.add((min(a, b), max(a, b)))
    return E

E = crust(S)
deg = np.bincount(np.array(list(E)).ravel(), minlength=len(S))
order = np.argsort(th)                           # true cyclic order along the curve
pos = np.empty(len(S), int); pos[order] = np.arange(len(S))
good = sum((pos[a] - pos[b]) % len(S) in (1, len(S) - 1) for a, b in E)
print(f"{len(E)} crust edges, all degrees 2: {np.all(deg == 2)}, "
      f"edges joining consecutive samples: {good}")

# The naive alternative: join each sample to its 2 nearest neighbours
D = np.linalg.norm(S[:, None] - S[None], axis=-1); np.fill_diagonal(D, np.inf)
nn2 = np.argsort(D, axis=1)[:, :2]
bad = sum((pos[i] - pos[j]) % len(S) not in (1, len(S) - 1) for i in range(len(S)) for j in nn2[i])
print("2-nearest-neighbour graph: wrong links =", bad)

Meshes

A surface reconstruction is usually delivered as a triangle mesh.

Implicit surfaces

An implicit method computes a function f : \mathbb{R}^3 \to \mathbb{R} that is negative inside and positive outside, and extracts the surface as its zero set. Noise is averaged and holes are filled by the smoothness of f.

Poisson surface reconstruction

Kazhdan, Bolitho and Hoppe (2006) reconstruct the indicator function \chi of the solid (1 inside, 0 outside) from oriented samples. The key observation is that the gradient of a smoothed indicator is a smoothed field of normals.

The field V is in general not a gradient, so we look for the function whose gradient is closest to it.

In practice the equation is discretised on an adaptive octree with B-spline bases, which makes the linear system sparse and its size proportional to the surface, not the volume. The iso-value \gamma is the average of \chi at the samples, which absorbs any global scaling of V.

Screened Poisson reconstruction (Kazhdan & Hoppe, 2013) adds a term pulling \chi towards \gamma at the samples, which reduces the over-smoothing of thin details.

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
import contourpy

rng = np.random.default_rng(0)
# Oriented samples of the "flower" r = 1 + 0.3 cos(5 theta), with noise
th = np.sort(rng.uniform(0, 2 * np.pi, 400))
r, dr = 1 + 0.3 * np.cos(5 * th), -1.5 * np.sin(5 * th)
P = np.c_[r * np.cos(th), r * np.sin(th)]
T = np.c_[dr * np.cos(th) - r * np.sin(th), dr * np.sin(th) + r * np.cos(th)]
Nout = np.c_[T[:, 1], -T[:, 0]] / np.linalg.norm(T, axis=1, keepdims=True)
P += 0.01 * rng.standard_normal(P.shape)
ds = 0.5 * (np.linalg.norm(np.roll(P, -1, 0) - P, axis=1)
            + np.linalg.norm(np.roll(P, 1, 0) - P, axis=1))   # arc length per sample

# 1. Splat the inward normals on a grid: V ~ grad(indicator)
m, L = 161, 2.0
g = np.linspace(-L, L, m); h = g[1] - g[0]
XX, YY = np.meshgrid(g, g, indexing="ij")
sig = 1.5 * h
V = np.zeros((m, m, 2))
for p, n, w in zip(P, -Nout, ds):
    G = np.exp(-((XX - p[0]) ** 2 + (YY - p[1]) ** 2) / (2 * sig ** 2)) / (2 * np.pi * sig ** 2)
    V += w * G[..., None] * n

# 2. Solve the Poisson equation  Laplacian(chi) = div V, chi = 0 on the border
div = np.gradient(V[..., 0], h, axis=0) + np.gradient(V[..., 1], h, axis=1)
I = sp.identity(m - 2)
D2 = sp.diags([1, -2, 1], [-1, 0, 1], shape=(m - 2, m - 2), dtype=float) / h ** 2
Lap = sp.kron(D2, I) + sp.kron(I, D2)
chi = np.zeros((m, m))
chi[1:-1, 1:-1] = spsolve(Lap.tocsc(), div[1:-1, 1:-1].ravel()).reshape(m - 2, m - 2)

# 3. Iso-value = mean of chi at the samples, then extract the curve
ij = np.clip(np.round((P + L) / h).astype(int), 0, m - 1)
iso = chi[ij[:, 0], ij[:, 1]].mean()
line = max(contourpy.contour_generator(XX, YY, chi).lines(iso), key=len)
x, y = line[:, 0], line[:, 1]
area = 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
print(f"chi inside ~ {chi[m // 2, m // 2]:.3f}, outside ~ {chi[5, 5]:.3f}, iso = {iso:.3f}")
print(f"area of the reconstruction {area:.4f} vs exact {np.pi * 1.045:.4f}")
rho = np.hypot(x, y); phi = np.arctan2(y, x)
print(f"max radial error: {np.max(np.abs(rho - (1 + 0.3 * np.cos(5 * phi)))):.4f}")

7. Rendering point clouds and meshes

Every step of the pipeline is checked by looking at the data, and many applications end with an image: inspection, virtual visits, simulation. This chapter explains how a mesh or a cloud becomes pixels, why a naive point display fails, and how splatting and level-of-detail hierarchies fix it.

Meshes: rasterisation and shading

A triangle is rasterised by testing pixel centres against its three edges. Inside the triangle, attributes such as depth, colour and normal are interpolated from the vertices.

Shading needs oriented normals: this is why the orientation step of ch. 2 matters even when no surface is reconstructed.

Rendering points directly

A cloud can be rendered without a mesh by giving each point an extent on the surface.

import numpy as np
from scipy.spatial import cKDTree
rng = np.random.default_rng(0)

# 3000 surfels on the unit sphere (position = normal), seen by a pinhole camera
N = 3000
X = rng.normal(size=(N, 3)); X /= np.linalg.norm(X, axis=1, keepdims=True)
W, f, cam = 200, 300.0, np.array([0.0, 0.0, -4.0])   # image size, focal (px), centre
Z = X[:, 2] - cam[2]                                  # depth along the optical axis
u, v = f * X[:, 0] / Z + W / 2, f * X[:, 1] / Z + W / 2
back = ((X - cam) * X).sum(1) > 0                     # normal points away from the camera

def splat(radius_px):
    """Z-buffer: each point covers a square of half-width radius_px (in pixels)."""
    depth, who = np.full((W, W), np.inf), np.full((W, W), -1)
    for k in range(N):
        r = int(np.ceil(radius_px[k]))
        i0, i1 = max(int(v[k]) - r, 0), min(int(v[k]) + r + 1, W)
        j0, j1 = max(int(u[k]) - r, 0), min(int(u[k]) + r + 1, W)
        closer = Z[k] < depth[i0:i1, j0:j1]
        depth[i0:i1, j0:j1][closer] = Z[k]
        who[i0:i1, j0:j1][closer] = k
    return who

# True silhouette: pixels whose viewing ray meets the sphere
jj, ii = np.meshgrid(np.arange(W) + 0.5, np.arange(W) + 0.5)
d = np.stack([(jj - W / 2) / f, (ii - W / 2) / f, np.ones_like(jj)], -1)
d /= np.linalg.norm(d, axis=-1, keepdims=True)
inside = (d @ cam) ** 2 - (cam @ cam - 1) > 0

s = np.sqrt(4 * np.pi / N)                            # mean spacing if samples were regular
d4 = cKDTree(X).query(X, k=5)[0][:, -1]               # distance to the 4th neighbour
for name, rad in [("1-pixel points", np.zeros(N)),
                  ("fixed r = f s / 2z", f * s / (2 * Z)),
                  ("adaptive r = f d4 / z", f * d4 / Z)]:
    who = splat(rad)
    holes = np.mean(who[inside] < 0)
    seen = who[inside & (who >= 0)]
    print(f"{name:22s} holes {100 * holes:5.1f} % | pixels showing a back-facing "
          f"point {100 * np.mean(back[seen]):5.1f} %")

Square splats of constant colour give a blocky image and aliasing at a distance. Gaussian kernels give smooth blending, and their projection has a closed form.

Level of detail

Large clouds hold billions of points, far more than the pixels of a screen.

The number of splats drawn depends on the screen area covered and on \tau, not on the size of the cloud. Lowering \tau while the camera is still refines the image progressively. Octree-based viewers stream the nodes from disk in the same way.

8. Research directions and open problems

This chapter is for self-study. It points to active research lines that extend each chapter of the course, with the questions that remain open. It is also a reservoir of papers for the course project.

Learned and certifiable registration, LiDAR odometry

Learning has entered registration at two levels. Learned descriptors replace FPFH in the pipeline of Definition 3.11. End-to-end networks go further: Deep Closest Point (Wang & Solomon, 2019) computes soft correspondences with attention and back-propagates through the SVD of Theorem 3.2.

Predator (Huang et al., 2021) focuses on scans with low overlap, and GeoTransformer (Qin et al., 2022) encodes pairwise distances and angles, which are invariant to rigid motions.

At the other end, TEASER (Yang, Shi & Carlone, 2021) minimises a truncated least-squares cost, robust to a large majority of outliers. A semidefinite relaxation can then certify that the solution found is globally optimal.

For mobile robots and mapping vehicles, registration becomes LiDAR odometry. LOAM (Zhang & Singh, 2014) registers edge and plane features. IMLS-SLAM (Deschaud, 2018) registers each scan to an implicit MLS surface of the map (Definition 6.9).

CT-ICP (Dellenbach, Deschaud, Jacquet & Goulette, 2022) estimates a pose at the beginning and at the end of each scan, which corrects the motion distortion of Note 1.10. KISS-ICP (Vizzo et al., 2023) shows that a carefully tuned point-to-point ICP with an adaptive threshold is competitive with far more complex systems.

Open questions include:

• degenerate places (tunnels, corridors, open fields), detected through the rank of the matrix of Proposition 3.7 but not yet handled reliably;
• moving objects, which violate the rigidity assumption;
• learned descriptors that transfer between sensors of different resolutions (Proposition 1.6);
• certifiable methods fast enough for real time.

Large-scale semantic segmentation

Outdoor scans contain billions of points. Superpoint graphs (Landrieu & Simonovsky, 2018) first partition the cloud into geometrically homogeneous parts, then classify the parts with a graph network. RandLA-Net (Hu et al., 2020) shows that random subsampling, combined with local feature aggregation, is enough to process large scenes efficiently.

Transformer architectures now lead most benchmarks; Point Transformer V3 (Wu et al., 2024) trades exact neighbourhoods for serialised point orderings to scale up.

Annotation is the bottleneck: labelling a scan point by point takes hours. Self-supervised pre-training on unlabelled clouds, by contrastive learning (PointContrast, Xie et al., 2020) or masked auto-encoding (Point-MAE, Pang et al., 2022), reduces the labels needed.

Domain shift between sensors, cities and seasons, long-tailed classes and calibrated uncertainty remain largely open. Work by Thomas, Deschaud, Goulette and coauthors (KPConv, 2019) and the Paris-Lille-3D data set (Roynard, Deschaud & Goulette, 2018) are reference points for this line.

Neural implicit surfaces

Occupancy Networks (Mescheder et al., 2019) and DeepSDF (Park et al., 2019) represent a shape by a neural network f_\theta(x), an occupancy probability or a signed distance, learned over collections of shapes. The surface is extracted by marching cubes (Definition 6.8).

A true signed distance satisfies the eikonal equation \|\nabla f\| = 1. Gropp and coauthors (2020) fit a network to a single raw cloud by asking f to vanish at the samples while penalising the deviation from this equation: \min_\theta \ \sum_i |f_\theta(x_i)| + \lambda\, \mathbb{E}_x \big( \|\nabla f_\theta(x)\| - 1 \big)^2 .

Shape As Points (Peng et al., 2021) makes Poisson reconstruction (Theorem 6.11) differentiable by solving it spectrally, so that the oriented points themselves can be learned.

Open questions: theoretical guarantees comparable to those of \varepsilon-samples, sharp features, and very large scenes.

Neural rendering and Gaussian splatting

NeRF (Mildenhall et al., 2020) represents a scene by a network mapping a position and a direction to a density and a colour, and renders it by integrating along rays. It produces photorealistic views but is slow to train and to render.

3D Gaussian splatting (Kerbl, Kopanas, Leimkühler & Drettakis, 2023) returns to points. The scene is a cloud of anisotropic 3D Gaussians, initialised from the sparse cloud of structure-from-motion and projected with the EWA formula of Proposition 7.8. Each Gaussian i has a projected opacity \alpha_i and a colour c_i.

Gaussian splatting renders in real time at high quality. Extracting accurate surfaces from it, controlling its memory footprint, and relighting are active questions; the point-based graphics of ch. 7 is its direct ancestor.

Generating and comparing point clouds

Generative models of point clouds, based on normalising flows (PointFlow, Yang et al., 2019) or diffusion (Luo & Hu, 2021), need a distance between clouds to train and to evaluate. Optimal transport gives the earth mover’s distance; the cheaper Chamfer distance is the most common choice.

Exercises

1. ★ Depth error of a stereo rig. A stereo rig has f = 1000 px and baseline b = 0.3 m, and matches disparities to 0.25 px. Compute the depth error at 10 m and at 20 m. Which baseline gives 5 cm at 20 m?

   Solution: by Proposition 1.2, \delta z = z^2 \delta d / (f b), which gives 100 \times 0.25 / 300 \approx 8.3 cm at 10 m and four times more, 33 cm, at 20 m.

   Solving for the baseline: b = \frac{z^2\, \delta d}{f\, \delta z} = \frac{400 \times 0.25}{1000 \times 0.05} = 2\ \text{m}.

2. ★ Phase-shift ambiguity. A phase-shift scanner modulates at 5 MHz. What is its ambiguity interval? It adds a second frequency of 4.5 MHz: over which range can the two phases together identify the distance?

   Solution: the interval is c / (2 f_m) = 3 \times 10^8 / 10^7 = 30 m. The pair of phases repeats only when both do, which happens with the beat frequency f_1 - f_2 = 0.5 MHz: the combined interval is \frac{c}{2 (f_1 - f_2)} = 300\ \text{m}.

3. ★ Farthest point sampling by hand. On the line, take the points 0, 1, 2, 3, 10. Run FPS from 0 with k = 3 and compute the covering radius. Compare with the best possible radius with 3 centres anywhere on the line.

   Solution: FPS picks 0, then 10 (the farthest), then 3 (at distance 3 from the set \{0, 10\}). The points 1 and 2 are then at distance 1, so \rho = 1.

   Centres at 0.5, 2.5 and 10 give \rho^\star = 0.5. The ratio is exactly 2: the bound of Proposition 2.4 is tight.

4. ★ A normal by hand. Compute the PCA normal and the surface variation of the four points (1, 0, 0), (-1, 0, 0), (0, 1, 0.1), (0, -1, -0.1).

   Solution: the centroid is 0. The covariance has C_{xx} = C_{yy} = 0.5, C_{zz} = 0.005, C_{yz} = 0.05, and 0 elsewhere.

   The (y, z) block has trace 0.505 and determinant 0, so the eigenvalues are 0.5, 0.505, 0. The eigenvector of 0 solves 0.5 y + 0.05 z = 0: n = \frac{(0, -0.1, 1)}{\sqrt{1.01}}, \qquad \sigma = 0 . The four points lie exactly on the plane z = 0.1\, y.

5. ★ RANSAC budgets. With p = 0.99, how many iterations are needed to detect a plane holding 20 % of the points with 3-point samples? A cylinder holding 20 % with samples of 2 oriented points? A cylinder holding 5 %?

   Solution: Theorem 3.12 gives N = 574 for the plane (w^3 = 0.008), N = 113 for the cylinder (w^2 = 0.04), and N = 1{,}840 for the cylinder at 5 %. Smaller minimal samples, made possible by normals, pay off quickly.

6. ★ Pooling and invariance. For a shared map h, compare the pooled features \max_i h(x_i), \sum_i h(x_i) and \frac1n \sum_i h(x_i). Which are permutation invariant? Which are unchanged when every point is duplicated?

   Solution: all three are permutation invariant, since max, sum and mean are symmetric. Duplicating every point leaves the max and the mean unchanged and doubles the sum, so max and mean are insensitive to a uniform change of density.

The mean still shifts when the density varies within the cloud, since densely sampled parts weigh more. The max depends only on which values are reached, not on how often, which suits the uneven clouds of ch. 1.

7. ★ Euler characteristic. A closed orientable triangle mesh has V = 10{,}000 and F = 20{,}000. What is its genus? Same question with V = 5{,}000 and F = 9{,}996.

   Solution: E = 3F/2 = 30{,}000, so V - E + F = 0 = 2 - 2g and g = 1, a torus. In the second case E = 14{,}994 and V - E + F = 2, so g = 0, a sphere (Proposition 6.5).

8. ★ Perspective-correct interpolation. A segment goes from a point A at depth 1 (attribute a = 0) to a point B at depth 3 (a = 1). Which value of a appears at the midpoint of its image?

   Solution: interpolate 1/Z and a/Z on the screen (Proposition 7.4): 1/Z = (1 + 1/3)/2 = 2/3 and a/Z = (0 + 1/3)/2 = 1/6, so a = \frac{1/6}{2/3} = 0.25 . The image midpoint shows the point a quarter of the way along the segment, at depth 1.5, not the true midpoint.

9. ★★ Optimal neighbourhood radius. Using Proposition 2.8, compute the optimal radius for a scanner noise of 5 mm on a surface with radius of curvature 2 m, and the corresponding value of \lambda_3 / \lambda_1. Why does a real implementation also need a lower bound on the number of points?

   Solution: with \kappa = 0.5 m⁻¹ and \sigma = 0.005 m, r^\star = 48^{1/4} \sqrt{0.01} \approx 0.26 m.

   At the optimum, ar^2 = b/r^2, so the ratio is 2\sqrt{ab} = 2 \kappa \sigma / \sqrt{3} \approx 0.0029. The proposition assumes uniform dense sampling; with few points the eigenvalues are themselves noisy, so the radius must also contain enough points (Mitra & Nguyen, 2003).

10. ★★ Statistical outlier removal (Python). Generate a 2 m × 2 m noisy plane of 4,750 points and 250 uniform outliers in a box around it. Apply the filter of Definition 2.12 with 8 neighbours and \alpha \in \{1, 2, 3\}, and report precision and recall.

    Solution:

    import numpy as np
    from scipy.spatial import cKDTree
    
    rng = np.random.default_rng(0)
    plane = np.c_[rng.uniform(0, 2, (4750, 2)), 0.005 * rng.standard_normal(4750)]
    noise = rng.uniform([0, 0, -0.5], [2, 2, 0.5], (250, 3))       # 5 % outliers
    X = np.r_[plane, noise]
    is_out = np.r_[np.zeros(4750, bool), np.ones(250, bool)]
    
    d, _ = cKDTree(X).query(X, k=9)                  # itself + 8 neighbours
    dbar = d[:, 1:].mean(1)                          # mean distance to the 8 neighbours
    mu, s = dbar.mean(), dbar.std()
    for alpha in (1.0, 2.0, 3.0):
        removed = dbar > mu + alpha * s
        prec = (removed & is_out).sum() / max(removed.sum(), 1)
        rec = (removed & is_out).sum() / is_out.sum()
        print(f"alpha = {alpha}: removed {removed.sum():4d}, precision {prec:.3f}, recall {rec:.3f}")

    Precision is 1.000 for the three values: no plane point is removed. Recall falls from 0.844 (\alpha = 1) to 0.728 and 0.616. The outliers that survive lie close to the plane, where they are geometrically indistinguishable from noise; projection smoothing then absorbs them.

11. ★★ Orienting normals by a spanning tree (Python). Implement Proposition 2.6: PCA normals with random signs, a minimum spanning tree of the 10-nearest-neighbour graph with weights 1 - |n_i^\top n_j|, and propagation from the highest point. Test on a sphere and on a thin plate whose two faces are 2 cm apart.

    Solution:

    import numpy as np
    from scipy.spatial import cKDTree
    from scipy.sparse import coo_matrix
    from scipy.sparse.csgraph import minimum_spanning_tree, breadth_first_order
    
    rng = np.random.default_rng(0)
    
    def orient(X, k=10):
        """PCA normals with random signs, then Hoppe-style propagation along an MST."""
        _, idx = cKDTree(X).query(X, k=k)
        N = np.zeros_like(X)
        for i, nb in enumerate(idx):
            Q = X[nb] - X[nb].mean(0)
            N[i] = np.linalg.eigh(Q.T @ Q)[1][:, 0]
        N *= rng.choice([-1.0, 1.0], (len(X), 1))                  # unknown signs
        rows = np.repeat(np.arange(len(X)), k - 1); cols = idx[:, 1:].ravel()
        w = 1 - np.abs(np.einsum("ij,ij->i", N[rows], N[cols])) + 1e-9
        T = minimum_spanning_tree(coo_matrix((w, (rows, cols)), shape=(len(X),) * 2))
        seed = np.argmax(X[:, 2])                                  # highest point: normal up
        if N[seed, 2] < 0:
            N[seed] *= -1
        order, pred = breadth_first_order(T, seed, directed=False)
        for v in order[1:]:
            if N[v] @ N[pred[v]] < 0:
                N[v] *= -1
        return N
    
    S = rng.normal(size=(3000, 3)); S /= np.linalg.norm(S, axis=1, keepdims=True)
    print("sphere: outward normals", np.mean(np.einsum("ij,ij->i", orient(S), S) > 0))
    
    # A thin plate: two parallel faces 2 cm apart, sampled every ~3 cm
    u = rng.uniform(0, 1, (1500, 2))
    P = np.r_[np.c_[u, np.zeros(1500)], np.c_[rng.uniform(0, 1, (1500, 2)), np.full(1500, 0.02)]]
    out = np.r_[-np.ones(1500), np.ones(1500)]                     # true outward n_z
    Np = orient(P)
    print("thin plate: correct normals", np.mean(np.sign(Np[:, 2]) == out))

    On the sphere, 100 % of the normals point outwards. On the plate only 49.6 % are correct: the neighbourhood graph joins the two faces, so propagation gives both the same sign and one face is entirely wrong. The fix is a neighbourhood smaller than the thickness, or the sensor position.

12. ★★ The best rotation in the plane. Show that for 2D point sets the minimum of E in Theorem 3.2 is \sum_i \|p_i'\|^2 + \sum_i \|q_i'\|^2 - 2 \sqrt{(H_{xx} + H_{yy})^2 + (H_{xy} - H_{yx})^2} , and evaluate it on Example 3.3.

    Solution: by the proof of Theorem 3.2 we maximise A \cos\theta + B \sin\theta with A = H_{xx} + H_{yy} and B = H_{xy} - H_{yx}. Its maximum is \sqrt{A^2 + B^2} by the Cauchy–Schwarz inequality, at \theta^\star = \operatorname{atan2}(B, A).

    In Example 3.3, both sums of squared norms equal 4/3 and \sqrt{A^2 + B^2} = 4/3, so the minimum is 4/3 + 4/3 - 8/3 = 0: the motion is exact.

13. ★★ Point-to-plane ICP on a plane. Suppose all target normals are equal to n and the source points span a plane orthogonal to n. Show that the matrix \sum_i a_i a_i^\top of Proposition 3.7 has rank exactly 3, and describe its kernel.

    Solution: a_i^\top \xi = p_i^\top (n \times \omega) + n^\top t.

    It vanishes for every i when \omega is parallel to n and t is orthogonal to n: a 3-dimensional kernel (rotation about n, translations in the plane).

    Conversely, if a_i^\top \xi = 0 for all i, the affine function p \mapsto p^\top (n \times \omega) + n^\top t vanishes on the plane.

    Since n \times \omega is orthogonal to n, this forces it to be 0, then n^\top t = 0. The kernel is exactly the space above, so the rank is 6 - 3 = 3.

14. ★★ Weights of robust ICP. Write the Huber and Cauchy functions as \rho(r) = \phi(r^2), check that \phi is concave, and compute the IRLS weights of Proposition 3.10. How much does a residual of 3c count relative to a perfect match?

    Solution: for Huber, \phi(s) = s/2 for s \le c^2 and c\sqrt{s} - c^2/2 beyond; for Cauchy, \phi(s) = \frac{c^2}{2} \log(1 + s / c^2). Both are concave and non-decreasing. Up to a factor 2 the weights are w_{\text{Huber}}(r) = \min\Big(1, \frac{c}{|r|}\Big), \qquad w_{\text{Cauchy}}(r) = \frac{1}{1 + r^2/c^2} . At r = 3c: one third for Huber and one tenth for Cauchy.

15. ★★ A graph cut on a chain. Four points in a row have data costs (U_i(0), U_i(1)) equal to (0, 2), (0.5, 1), (1.5, 0.4) and (2, 0), and neighbours are consecutive points. Find the labelling of minimal Potts energy (Definition 4.3) for \beta = 0, 1, 2 and 3.

    Solution: a labelling with two or more changes is never optimal here, so compare the monotone ones. With \beta = 0, choose the cheaper label pointwise: 0011, energy 0.9.

    For \beta = 1, 0011 costs 0.9 + 1 = 1.9, against 2.4 for 0111, 3 for 0001 and 3.4 for 1111. For \beta = 2, 0011 still wins with 2.9.

    For \beta = 3, 0011 costs 3.9 and the constant labelling 1111, at 3.4, wins: a strong prior erases the small region.

16. ★★ The budget of a Hough transform. For a 20 m façade and angle cells of 0.5°, over how many centimetres of \rho do its votes spread (Proposition 5.7)? How many votes do 10^7 points cast over a hemisphere of normals?

    Solution: 20 \times 0.0087 / 2 \approx 8.7 cm, plus the noise. The hemisphere has 720 \times 180 = 129{,}600 angle cells, hence about 1.3 \times 10^{12} votes. This is why practical detectors use RANSAC, randomised Hough transforms or votes restricted by normals: a point with a known normal votes for a single angle cell.

17. ★★ Screened Poisson. Add to the functional of Theorem 6.11 the term \alpha \int_\Omega w\, (\chi - \gamma)^2, where w \ge 0 concentrates near the samples. Derive the new Euler–Lagrange equation and explain its effect.

    Solution: the first variation gains 2\alpha \int_\Omega w (\chi - \gamma) \varphi, so the equation becomes \Delta \chi - \alpha\, w\, (\chi - \gamma) = \nabla \cdot V . Near the samples, \chi is pulled towards the iso-value \gamma, so the extracted surface passes closer to the data. This is the screened Poisson equation of Kazhdan and Hoppe (2013), which preserves thin details that plain Poisson smooths out.

18. ★★★ The basin of convergence of ICP. Using Code 3.8, measure the success rate of point-to-point and point-to-plane ICP (final rotation error below 1°) as a function of the initial rotation angle, from 10° to 90°, and of the overlap, obtained by cropping the source. Propose and test a strategy that enlarges the basin.

    Hints and expected result: run 20 random axes per angle. With 30 iterations and full overlap, both variants succeed up to about 50°. At 70° point-to-plane still succeeds in most runs while point-to-point fails in about half, and at 90° both fail in most runs. Reducing the overlap makes the drop earlier unless pairs are trimmed (Definition 3.9).

    Multi-start over a grid of initial rotations, coarse-to-fine subsampling (Definition 2.3), or a global initialisation by FPFH and RANSAC (Definition 3.11) all enlarge the basin. Report a curve of success rate against angle for each strategy.

19. ★★★ From EWA to Gaussian splatting. A 3D Gaussian of covariance \Sigma is seen by a camera of rotation W, with J the Jacobian of the perspective map at its centre. Show that its screen footprint has covariance J W \Sigma W^\top J^\top, the upper 2 \times 2 block being kept. Then write a small 2D experiment that fits a set of 2D Gaussians to a target image by gradient descent on Definition 8.1.

    Hints and expected result: the first part is Proposition 7.8 applied twice: the rotation W maps \Sigma to W \Sigma W^\top, and the local affine approximation of the projection maps it again by J. For the experiment, in PyTorch, parametrise each covariance as L L^\top with L lower triangular, so that it stays positive semi-definite.

    Sort the Gaussians by a learned depth, composite them, and minimise the squared error to a synthetic image of a few disks. With about 50 Gaussians the error should fall steadily and the disks appear; Kerbl et al. (2023) add densification and pruning of Gaussians, which you can imitate.

Projects and assessment

The official assessment (master-mva.com) is a project based on a research paper, in four steps: choose a paper from the instructors’ proposals, implement and test its algorithm(s), analyse the algorithm and its results and propose variants when relevant, and hand in a document explaining the method and the results together with the code. A student who fails may take a new paper from a new list.

The Université Paris-Saclay course detail describes the assessment as partial continuous assessment; no weights are displayed. The solo subjects below follow the four official steps on named papers. The group subjects train the same skills at a larger scale. Every subject is evaluated with the rubric at the end of this section, on its stated success criterion; the written exam is a self-assessment proposal, not part of the official assessment.

Solo projects

1. [Solo · ≈ 20 h] Paper study: Poisson surface reconstruction. Reproduce Kazhdan, Bolitho and Hoppe (2006) on synthetic shapes.
   • Question: how does accuracy depend on noise, sampling density and grid resolution, compared with the signed distance of Hoppe et al. (1992)?
   • Work: extend Code 6.12 to 3D on a regular grid of 64^3 to 128^3 cells; extract the surface with marching cubes (Definition 6.8, implemented or taken from a library); test on a sphere, a torus and a union of boxes with 0–2 % noise.
   • Success criterion: square root of the Chamfer distance (Definition 8.2) to dense samples of the true surface below 1 % of the bounding-box diagonal at 1 % noise, and a curve of error against noise for both methods.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 2 and 6.
2. [Solo · ≈ 20 h] Paper study: efficient RANSAC. Reimplement Schnabel, Wahl and Klein (2007).
   • Work: octree-localised sampling (Proposition 5.5), planes, spheres and cylinders from oriented points (Propositions 5.2–5.3), inlier connectivity, and plain sequential RANSAC as a baseline.
   • Data: synthetic scenes of 10^5 to 10^6 points (a room with pipes, a tank, columns) with clutter, where the true primitives are known.
   • Success criterion: every primitive detected with at least 90 % of its points, and a speed-up of at least 10 over the baseline at 10^6 points.
   • Deliverables: code, report, talk. Chapters 2 and 5.
3. [Solo · ≈ 20 h] Paper study: Generalized-ICP. Reimplement Segal, Haehnel and Thrun (2009) and compare it with point-to-point and point-to-plane ICP.
   • Work: local covariances from ch. 2, the GICP cost minimised by Gauss–Newton on rigid motions, and two scans simulated from different poses with Code 1.7.
   • Success criterion: success-rate curves against the initial rotation (as in Exercise 18) and final errors for the three variants, with the paper’s qualitative ordering confirmed or explained.
   • Deliverables: code, report, talk. Chapters 1, 2 and 3.

Group projects

1. [Group of 2–3 · ≈ 30 h each] LiDAR odometry on simulated drives. Build a scan-to-map odometry in the spirit of IMLS-SLAM (Deschaud, 2018) and CT-ICP (Dellenbach et al., 2022).
   • Work: a simulator of a vehicle driving a synthetic street, extending Code 1.7 with motion distortion; scan-to-map point-to-plane ICP on a voxel map; per-point motion compensation.
   • Split: one member writes the simulator and the ground truth, one the registration and the map, one the motion compensation and the evaluation.
   • Success criterion: drift below 1 % of the distance travelled on a 200 m loop with 2 cm noise; an ablation without motion compensation.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 1, 2 and 3.
2. [Group of 2–3 · ≈ 30 h each] Semantic segmentation of an urban scan. Compare multi-scale features with a random forest (Weinmann et al., 2015; Hackel, Wegner & Schindler, 2016) with a point network (PointNet++ or KPConv, Thomas et al., 2019).
   • Data: a public urban LiDAR data set such as Paris-Lille-3D (Roynard, Deschaud & Goulette, 2018), with training and test areas kept apart; the synthetic scenes of Code 4.6 for debugging.
   • Split: features and forest; network (free GPU tier); graph-cut regularisation (Theorem 4.4) and evaluation.
   • Success criterion: per-class IoU and mIoU on the held-out area for the three pipelines, with an analysis of the classes where each one fails.
   • Deliverables: code, report, talk. Chapters 2 and 4.
3. [Group of 2–3 · ≈ 30 h each] A point-based renderer with level of detail. Implement the hierarchy of QSplat (Rusinkiewicz & Levoy, 2000) with Gaussian splats (Proposition 7.8) and normalised blending.
   • Data: synthetic clouds of 10^6 to 10^7 points sampled on analytic surfaces with known normals, or a scan from the Stanford 3D Scanning Repository.
   • Split: hierarchy and traversal; splatting, blending and shading; evaluation and an interactive viewer.
   • Success criterion: images without holes or show-through (as measured in Code 7.7), and curves of image error against a full-resolution reference and of the number of splats against the threshold \tau.
   • Deliverables: code, report, talk with a live demonstration. Chapters 2 and 7.

Written exam (proposal)

Format: 2 hours, one handwritten A4 sheet allowed, no electronic devices, four exercises.

1. Question (4 points). Compare triangulation and time-of-flight sensors: give the range formula of each and the depth error of triangulation, and say which you would use to scan a 50 m façade, and why. Expected: z = fb/d and \delta z \approx z^2 \delta d / (fb); pulsed r = c\tau/2 and phase shift with its ambiguity interval; for 50 m the quadratic growth rules out short-baseline triangulation, so a pulsed or multi-frequency phase-shift scanner is the answer.

2. Question (6 points). (a) Align the points (0, 0), (2, 0), (0, 1) to the points (3, 1), (3, −1), (4, 1) in closed form. (b) A global registration finds 40 % of correct matches: how many 3-match RANSAC draws give a 99 % success probability? Expected: centroids (2/3, 1/3) and (10/3, 1/3); H_{xx} + H_{yy} = 0 and H_{xy} - H_{yx} = -10/3, so \theta^\star = -90° and t^\star = (3, 1); (b) N = 70 by Theorem 3.12.

3. Question (5 points). Prove that the costs of ICP are non-increasing, and explain why ICP may nevertheless return a wrong alignment. Expected: the two inequalities of Theorem 3.5, one for matching and one for alignment; the limit is only a local minimum, whose basin depends on the initial rotation, the overlap and the symmetries of the scene; remedies are a global initialisation and robust weights.

4. Question (5 points). State the variational problem of Poisson reconstruction, derive its Euler–Lagrange equation, and explain why normals must be consistently oriented and how the iso-value is chosen. Expected: minimise the integral of \|\nabla\chi - V\|^2, first variation and integration by parts giving \Delta\chi = \nabla\cdot V; a flipped normal reverses its contribution to V and creates spurious surface sheets; the iso-value is the mean of \chi at the samples, which absorbs the scale of V.

Grading rubric (proposal)

References

Official references (master-mva.com)

• J. Gallice (ed.), Images de profondeur, Hermès Science, 2001.
• F. Goulette, Modélisation 3D automatique, Presses de l’École des Mines, 1998.
• J.-E. Deschaud, Traitement de nuages de points denses et modélisation 3D d’environnements par système mobile LiDAR/Caméra, PhD thesis, MINES ParisTech, 2010.
• M. Gross and H. Pfister (eds.), Point-Based Graphics, Morgan Kaufmann, 2007.

Textbooks and surveys

• M. Botsch, L. Kobbelt, M. Pauly, P. Alliez and B. Lévy, Polygon Mesh Processing, A K Peters, 2010.
• T. K. Dey, Curve and Surface Reconstruction: Algorithms with Mathematical Analysis, Cambridge University Press, 2007.
• R. Hartley and A. Zisserman, Multiple View Geometry in Computer Vision, 2nd ed., Cambridge University Press, 2004.
• G. Vosselman and H.-G. Maas (eds.), Airborne and Terrestrial Laser Scanning, Whittles Publishing, 2010.
• T. Akenine-Möller, E. Haines, N. Hoffman, A. Pesce, M. Iwanicki and S. Hillaire, Real-Time Rendering, 4th ed., CRC Press, 2018.
• M. Berger, A. Tagliasacchi, L. M. Seversky, P. Alliez, G. Guennebaud, J. A. Levine, A. Sharf and C. T. Silva, “A survey of surface reconstruction from point clouds”, Computer Graphics Forum 36(1), 2017.
• F. Pomerleau, F. Colas and R. Siegwart, “A review of point cloud registration algorithms for mobile robotics”, Foundations and Trends in Robotics 4(1), 2015.
• Y. Guo, H. Wang, Q. Hu, H. Liu, L. Liu and M. Bennamoun, “Deep learning for 3D point clouds: a survey”, IEEE Transactions on Pattern Analysis and Machine Intelligence 43(12), 2021.

Sensors, data structures and operators (ch. 1–2)

• J. L. Bentley, “Multidimensional binary search trees used for associative searching”, Communications of the ACM 18(9), 1975.
• J. H. Friedman, J. L. Bentley and R. A. Finkel, “An algorithm for finding best matches in logarithmic expected time”, ACM Transactions on Mathematical Software 3(3), 1977.
• T. F. Gonzalez, “Clustering to minimize the maximum intercluster distance”, Theoretical Computer Science 38, 1985.
• H. Hoppe, T. DeRose, T. Duchamp, J. McDonald and W. Stuetzle, “Surface reconstruction from unorganized points”, SIGGRAPH, 1992.
• M. Pauly, M. Gross and L. Kobbelt, “Efficient simplification of point-sampled surfaces”, IEEE Visualization, 2002.
• N. J. Mitra and A. Nguyen, “Estimating surface normals in noisy point cloud data”, ACM Symposium on Computational Geometry, 2003.
• J. Demantké, C. Mallet, N. David and B. Vallet, “Dimensionality based scale selection in 3D lidar point clouds”, ISPRS Archives XXXVIII-5/W12 (Laser Scanning), 2011.
• R. B. Rusu, N. Blodow and M. Beetz, “Fast point feature histograms (FPFH) for 3D registration”, IEEE ICRA, 2009.
• M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin and C. T. Silva, “Computing and rendering point set surfaces”, IEEE Transactions on Visualization and Computer Graphics 9(1), 2003.
• R. B. Rusu and S. Cousins, “3D is here: Point Cloud Library (PCL)”, IEEE ICRA, 2011.
• Q.-Y. Zhou, J. Park and V. Koltun, “Open3D: a modern library for 3D data processing”, arXiv:1801.09847, 2018.

Registration (ch. 3)

• K. S. Arun, T. S. Huang and S. D. Blostein, “Least-squares fitting of two 3-D point sets”, IEEE Transactions on Pattern Analysis and Machine Intelligence 9(5), 1987.
• B. K. P. Horn, “Closed-form solution of absolute orientation using unit quaternions”, Journal of the Optical Society of America A 4(4), 1987.
• S. Umeyama, “Least-squares estimation of transformation parameters between two point patterns”, IEEE Transactions on Pattern Analysis and Machine Intelligence 13(4), 1991.
• P. J. Besl and N. D. McKay, “A method for registration of 3-D shapes”, IEEE Transactions on Pattern Analysis and Machine Intelligence 14(2), 1992.
• Y. Chen and G. Medioni, “Object modelling by registration of multiple range images”, Image and Vision Computing 10(3), 1992.
• Z. Zhang, “Iterative point matching for registration of free-form curves and surfaces”, International Journal of Computer Vision 13(2), 1994.
• S. Rusinkiewicz and M. Levoy, “Efficient variants of the ICP algorithm”, 3-D Digital Imaging and Modeling (3DIM), 2001.
• D. Chetverikov, D. Svirko, D. Stepanov and P. Krsek, “The trimmed iterative closest point algorithm”, ICPR, 2002.
• K.-L. Low, “Linear least-squares optimization for point-to-plane ICP surface registration”, Technical Report TR04-004, University of North Carolina at Chapel Hill, 2004.
• P. Biber and W. Straßer, “The normal distributions transform: a new approach to laser scan matching”, IEEE/RSJ IROS, 2003.
• A. Segal, D. Haehnel and S. Thrun, “Generalized-ICP”, Robotics: Science and Systems, 2009.
• D. Aiger, N. J. Mitra and D. Cohen-Or, “4-points congruent sets for robust pairwise surface registration”, ACM Transactions on Graphics 27(3), 2008.
• H. Yang, J. Shi and L. Carlone, “TEASER: fast and certifiable point cloud registration”, IEEE Transactions on Robotics 37(2), 2021.
• P. J. Huber, “Robust estimation of a location parameter”, Annals of Mathematical Statistics 35(1), 1964.
• M. A. Fischler and R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography”, Communications of the ACM 24(6), 1981.

Segmentation and learning (ch. 4)

• T. Rabbani, F. A. van den Heuvel and G. Vosselman, “Segmentation of point clouds using smoothness constraint”, ISPRS Archives 36(5), 2006.
• D. M. Greig, B. T. Porteous and A. H. Seheult, “Exact maximum a posteriori estimation for binary images”, Journal of the Royal Statistical Society B 51(2), 1989.
• Y. Boykov, O. Veksler and R. Zabih, “Fast approximate energy minimization via graph cuts”, IEEE Transactions on Pattern Analysis and Machine Intelligence 23(11), 2001.
• L. Breiman, “Random forests”, Machine Learning 45(1), 2001.
• M. Weinmann, B. Jutzi, S. Hinz and C. Mallet, “Semantic point cloud interpretation based on optimal neighborhoods, relevant features and efficient classifiers”, ISPRS Journal of Photogrammetry and Remote Sensing 105, 2015.
• T. Hackel, J. D. Wegner and K. Schindler, “Fast semantic segmentation of 3D point clouds with strongly varying density”, ISPRS Annals III-3, 2016.
• C. R. Qi, H. Su, K. Mo and L. J. Guibas, “PointNet: deep learning on point sets for 3D classification and segmentation”, CVPR, 2017.
• C. R. Qi, L. Yi, H. Su and L. J. Guibas, “PointNet++: deep hierarchical feature learning on point sets in a metric space”, NeurIPS, 2017.
• M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Póczos, R. Salakhutdinov and A. Smola, “Deep sets”, NeurIPS, 2017.
• H. Thomas, C. R. Qi, J.-E. Deschaud, B. Marcotegui, F. Goulette and L. J. Guibas, “KPConv: flexible and deformable convolution for point clouds”, ICCV, 2019.
• Y. Wang, Y. Sun, Z. Liu, S. E. Sarma, M. M. Bronstein and J. M. Solomon, “Dynamic graph CNN for learning on point clouds”, ACM Transactions on Graphics 38(5), 2019.
• C. Choy, J. Gwak and S. Savarese, “4D spatio-temporal ConvNets: Minkowski convolutional neural networks”, CVPR, 2019.
• H. Zhao, L. Jiang, J. Jia, P. Torr and V. Koltun, “Point transformer”, ICCV, 2021.
• Z. Wu, S. Song, A. Khosla, F. Yu, L. Zhang, X. Tang and J. Xiao, “3D ShapeNets: a deep representation for volumetric shapes”, CVPR, 2015.
• I. Armeni, O. Sener, A. R. Zamir, H. Jiang, I. Brilakis, M. Fischer and S. Savarese, “3D semantic parsing of large-scale indoor spaces”, CVPR, 2016.
• T. Hackel, N. Savinov, L. Ladicky, J. D. Wegner, K. Schindler and M. Pollefeys, “Semantic3D.net: a new large-scale point cloud classification benchmark”, ISPRS Annals IV-1/W1, 2017.
• X. Roynard, J.-E. Deschaud and F. Goulette, “Paris-Lille-3D: a large and high-quality ground-truth urban point cloud dataset for automatic segmentation and classification”, International Journal of Robotics Research 37(6), 2018.
• J. Behley, M. Garbade, A. Milioto, J. Quenzel, S. Behnke, C. Stachniss and J. Gall, “SemanticKITTI: a dataset for semantic scene understanding of LiDAR sequences”, ICCV, 2019.

Primitives (ch. 5)

• R. Schnabel, R. Wahl and R. Klein, “Efficient RANSAC for point-cloud shape detection”, Computer Graphics Forum 26(2), 2007.
• D. Borrmann, J. Elseberg, K. Lingemann and A. Nüchter, “The 3D Hough transform for plane detection in point clouds: a review and a new accumulator design”, 3D Research 2(2), 2011.
• F. Lafarge and C. Mallet, “Creating large-scale city models from 3D-point clouds: a robust approach with hybrid representation”, International Journal of Computer Vision 99(1), 2012.
• L. Nan and P. Wonka, “PolyFit: polygonal surface reconstruction from point clouds”, ICCV, 2017.

Reconstruction (ch. 6)

• N. Amenta, M. Bern and D. Eppstein, “The crust and the β-skeleton: combinatorial curve reconstruction”, Graphical Models and Image Processing 60(2), 1998.
• N. Amenta and M. Bern, “Surface reconstruction by Voronoi filtering”, Discrete & Computational Geometry 22(4), 1999.
• N. Amenta, S. Choi and R. K. Kolluri, “The power crust”, ACM Symposium on Solid Modeling and Applications, 2001.
• N. Amenta, S. Choi, T. K. Dey and N. Leekha, “A simple algorithm for homeomorphic surface reconstruction”, International Journal of Computational Geometry & Applications 12(1–2), 2002.
• H. Edelsbrunner and E. P. Mücke, “Three-dimensional alpha shapes”, ACM Transactions on Graphics 13(1), 1994.
• F. Bernardini, J. Mittleman, H. Rushmeier, C. Silva and G. Taubin, “The ball-pivoting algorithm for surface reconstruction”, IEEE Transactions on Visualization and Computer Graphics 5(4), 1999.
• W. E. Lorensen and H. E. Cline, “Marching cubes: a high resolution 3D surface construction algorithm”, SIGGRAPH, 1987.
• G. M. Nielson and B. Hamann, “The asymptotic decider: resolving the ambiguity in marching cubes”, IEEE Visualization, 1991.
• J. C. Carr, R. K. Beatson, J. B. Cherrie, T. J. Mitchell, W. R. Fright, B. C. McCallum and T. R. Evans, “Reconstruction and representation of 3D objects with radial basis functions”, SIGGRAPH, 2001.
• R. Kolluri, “Provably good moving least squares”, ACM Transactions on Algorithms 4(2), 2008.
• M. Kazhdan, M. Bolitho and H. Hoppe, “Poisson surface reconstruction”, Eurographics Symposium on Geometry Processing, 2006.
• M. Kazhdan and H. Hoppe, “Screened Poisson surface reconstruction”, ACM Transactions on Graphics 32(3), 2013.

Rendering (ch. 7)

• E. Catmull, A Subdivision Algorithm for Computer Display of Curved Surfaces, PhD thesis, University of Utah, 1974.
• B. T. Phong, “Illumination for computer generated pictures”, Communications of the ACM 18(6), 1975.
• J. F. Blinn, “Models of light reflection for computer synthesized pictures”, SIGGRAPH, 1977.
• H. Pfister, M. Zwicker, J. van Baar and M. Gross, “Surfels: surface elements as rendering primitives”, SIGGRAPH, 2000.
• S. Rusinkiewicz and M. Levoy, “QSplat: a multiresolution point rendering system for large meshes”, SIGGRAPH, 2000.
• M. Zwicker, H. Pfister, J. van Baar and M. Gross, “Surface splatting”, SIGGRAPH, 2001.
• M. Botsch, A. Hornung, M. Zwicker and L. Kobbelt, “High-quality surface splatting on today’s GPUs”, Eurographics Symposium on Point-Based Graphics, 2005.
• T. Boubekeur and M. Alexa, “Phong tessellation”, ACM Transactions on Graphics 27(5), 2008.

Research directions (ch. 8)

• J. Zhang and S. Singh, “LOAM: lidar odometry and mapping in real-time”, Robotics: Science and Systems, 2014.
• J.-E. Deschaud, “IMLS-SLAM: scan-to-model matching based on 3D data”, IEEE ICRA, 2018.
• P. Dellenbach, J.-E. Deschaud, B. Jacquet and F. Goulette, “CT-ICP: real-time elastic LiDAR odometry with loop closure”, IEEE ICRA, 2022.
• I. Vizzo, T. Guadagnino, B. Mersch, L. Wiesmann, J. Behley and C. Stachniss, “KISS-ICP: in defense of point-to-point ICP — simple, accurate, and robust registration if done the right way”, IEEE Robotics and Automation Letters 8(2), 2023.
• Y. Wang and J. M. Solomon, “Deep closest point: learning representations for point cloud registration”, ICCV, 2019.
• S. Huang, Z. Gojcic, M. Usvyatsov, A. Wieser and K. Schindler, “PREDATOR: registration of 3D point clouds with low overlap”, CVPR, 2021.
• Z. Qin, H. Yu, C. Wang, Y. Guo, Y. Peng and K. Xu, “Geometric transformer for fast and robust point cloud registration”, CVPR, 2022.
• L. Landrieu and M. Simonovsky, “Large-scale point cloud semantic segmentation with superpoint graphs”, CVPR, 2018.
• Q. Hu, B. Yang, L. Xie, S. Rosa, Y. Guo, Z. Wang, N. Trigoni and A. Markham, “RandLA-Net: efficient semantic segmentation of large-scale point clouds”, CVPR, 2020.
• X. Wu, L. Jiang, P.-S. Wang, Z. Liu, X. Liu, Y. Qiao, W. Ouyang, T. He and H. Zhao, “Point Transformer V3: simpler, faster, stronger”, CVPR, 2024.
• S. Xie, J. Gu, D. Guo, C. R. Qi, L. Guibas and O. Litany, “PointContrast: unsupervised pre-training for 3D point cloud understanding”, ECCV, 2020.
• Y. Pang, W. Wang, F. E. H. Tay, W. Liu, Y. Tian and L. Yuan, “Masked autoencoders for point cloud self-supervised learning”, ECCV, 2022.
• L. Mescheder, M. Oechsle, M. Niemeyer, S. Nowozin and A. Geiger, “Occupancy networks: learning 3D reconstruction in function space”, CVPR, 2019.
• J. J. Park, P. Florence, J. Straub, R. Newcombe and S. Lovegrove, “DeepSDF: learning continuous signed distance functions for shape representation”, CVPR, 2019.
• A. Gropp, L. Yariv, N. Haim, M. Atzmon and Y. Lipman, “Implicit geometric regularization for learning shapes”, ICML, 2020.
• S. Peng, C. M. Jiang, Y. Liao, M. Niemeyer, M. Pollefeys and A. Geiger, “Shape as points: a differentiable Poisson solver”, NeurIPS, 2021.
• B. Mildenhall, P. P. Srinivasan, M. Tancik, J. T. Barron, R. Ramamoorthi and R. Ng, “NeRF: representing scenes as neural radiance fields for view synthesis”, ECCV, 2020.
• B. Kerbl, G. Kopanas, T. Leimkühler and G. Drettakis, “3D Gaussian splatting for real-time radiance field rendering”, ACM Transactions on Graphics 42(4), 2023.
• G. Yang, X. Huang, Z. Hao, M.-Y. Liu, S. Belongie and B. Hariharan, “PointFlow: 3D point cloud generation with continuous normalizing flows”, ICCV, 2019.
• S. Luo and W. Hu, “Diffusion probabilistic models for 3D point cloud generation”, CVPR, 2021.