RasmahRasmah

From points to mesh

Overview

A 3D scanner does not produce a mesh. It produces a point cloud — a million or so positions on the object's surface, with no triangles, no connectivity, no idea of inside versus outside. Reconstructing a usable mesh from those points is the inverse problem of meshing.

This chapter covers the two halves of that problem: estimating the missing surface information (normals), and then reconstructing a surface or volume from the points — either directly, or by first recognizing that the cloud is really just a sphere, a plane, or a box that someone scanned.

The point cloud

Rasmah.PointCloud — Type
PointCloud{P,N}

A 3×n point set with optional per-point normals.

Fields

  • points::P: a 3×n point matrix.
  • normals::N: a 3×n normal matrix, or nothing.
source

As a running example, take a cloud of points sampled on a unit sphere:

pts = surface_mesh(sphere(1.0), n=16).vertices
size(pts, 2)
162

Estimating normals

A point cloud has no surface normal at each point — you have to infer one from the point's neighbours. estimate_normals fits a local plane by PCA and takes its normal; orient_normals then flips them so they all point the same (outward) way:

Rasmah.estimate_normals — Function
estimate_normals(points::AbstractMatrix; k = 10) -> Matrix{Float64}

Estimate per-point unit normals of a 3×n point cloud by PCA of the local k-neighborhood (the eigenvector of the smallest covariance eigenvalue). The sign is arbitrary per point (see orient_normals).

source
Rasmah.orient_normals — Function
orient_normals(points::AbstractMatrix, normals::AbstractMatrix;
               k = 10, outward = true) -> Matrix{Float64}

Orient a set of per-point normals consistently by propagating across a minimum spanning tree of the k-neighbor graph. When outward=true, each component is additionally flipped so the majority of normals point away from its centroid.

source
normals = orient_normals(pts, estimate_normals(pts))
size(normals, 2)
162

Ball-pivoting

The simplest reconstruction method imagines a ball rolling over the point cloud, pivoting around each edge to find the next triangle:

Rasmah.ball_pivot — Function
ball_pivot(points; ρ = nothing) -> TriangleMesh

Ball-pivoting surface reconstruction: roll a virtual ball of radius ρ over the point cloud to build a triangle mesh that interpolates the input points. ρ = nothing auto-selects a radius from alpha_shape_radius.

source
vtk(ball_pivot(pts))

Screened Poisson reconstruction

A more robust method solves a Poisson equation whose gradient best matches the estimated normals, then extracts the surface where the implied indicator changes sign — giving a smooth, watertight mesh even from noisy input:

Rasmah.poisson_reconstruction — Function
poisson_reconstruction(points; normals=nothing, n=64, k=10, margin=0.05, screening=0.0, iso=nothing, maxit=2000) -> TriangleMesh

Screened Poisson surface reconstruction (Kazhdan & Hoppe 2013) from an oriented point cloud points (3×N). If normals are not given they are estimated with k-NN and oriented. Solves (Δ − screening·I)χ = ∇·V by conjugate gradients on an n³ grid and extracts the χ = iso iso-surface by marching cubes. Returns a watertight, outward-oriented TriangleMesh.

source
vtk(poisson_reconstruction(pts; normals=normals, n=32))

Alpha shapes

An alpha shape is the set of points plus the triangles whose circumradius is at most $\alpha$. A small $\alpha$ keeps only the tight triangles (the surface); a huge $\alpha$ fills in everything (the convex hull):

Rasmah.alpha_shape — Function
alpha_shape(points; α = nothing) -> TriangleMesh

Boundary surface of the 3D alpha solid: keep the Delaunay tetrahedra whose circumsphere radius is < α and return the faces bounding exactly one kept tetrahedron, oriented outward. α = nothing auto-selects via alpha_shape_radius.

source
vtk(alpha_shape(pts; α=1.0))

Delaunay tetrahedralization

For a volume reconstruction, the Delaunay tetrahedralization connects all the points into a tetrahedral mesh whose circumspheres contain no other point:

Rasmah.delaunay_tetrahedralization — Function
delaunay_tetrahedralization(points::AbstractMatrix; weights = nothing) -> TetMesh

Delaunay tetrahedralization of a 3×N matrix of points (tets positively oriented). Passing weights (length-N, or nothing) yields a regular (weighted / power) triangulation.

Returns

A TetMesh of the Delaunay tetrahedralization.

source
Rasmah.delaunay_tetmesh — Function
delaunay_tetmesh(g::ImplicitGeometry, n = 32; conform = true) -> TetMesh

Delaunay volume mesh of an implicit geometry: sample the surface and a coarse interior grid, tetrahedralize, keep the inside tets, and (optionally) snap the boundary onto the surface. n controls the surface/interior sampling resolution.

source
vtk(delaunay_tetrahedralization(pts))

Reverse engineering with RANSAC

Sometimes the fastest path is to notice the cloud is not a million points at all — it is a sphere (or a plane, or a cylinder) plus noise. RANSAC guesses a primitive from a few random points, counts how many points agree, and keeps the best fit:

Rasmah.detect_primitives — Function
detect_primitives(points; normals=nothing, max_primitives=4, min_support=10, distance_threshold=0.05, angle_threshold=deg2rad(5), max_iterations=500, seed=nothing)

Greedily detect up to max_primitives analytic primitives (planes, and — when normals are given — spheres, cylinders, cones) in a point cloud: run all detectors, keep the best-supported one, remove its inliers, and repeat on the remainder. Returns a (possibly heterogeneous) vector of PlaneFit/SphereFit/ CylinderFit/ConeFit.

source
fit = ransac_sphere(pts; seed=1)
SphereFit([5.178846433198768e-17, -1.1846767626105591e-16, 3.200703118233001e-17], 0.9999999999999999, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  153, 154, 155, 156, 157, 158, 159, 160, 161, 162], 162)

The fit can be turned straight back into a Feature node — recovering the parametric shape that was scanned, not just a triangle soup:

volume(to_feature(fit))
4.188790204786389 m^3
volume(sphere(1.0))
4.1887902047863905 m^3

The reconstructed sphere's analytic volume matches the original exactly ($4\pi/3 \approx 4.1888$).

Next steps

This completes the meshing story. Meshes feed directly into simulation: see the FEM & simulation chapters.