Segmenting a point cloud with a 2004 image algorithm and a statistical distance

This is the third post built on the same octree extension, and it is really one argument spread across three. The first post built a chunked octree that finds a point’s voxel with a handful of bit operations. The second put a running Gaussian, mean and covariance, in every voxel that structure indexes, updated one point at a time and never revisited. This one takes those Gaussians and asks a question about pairs of them: are two neighbouring voxels similar enough to be part of the same surface? Answering that with a real statistical distance rather than “how far apart are their centres” is the whole point of what follows, and the reason is concrete rather than abstract: a voxel sitting on a floor right next to a corner and a voxel sitting on a wall right next to the same corner can be centimetres apart, closer to each other than either is to a voxel in the middle of its own plane, even though the two surfaces face completely different directions.

The merging itself is not new. Felzenszwalb and Huttenlocher’s 2004 graph-based segmentation, already covered in full in post 69, sorts every edge by weight and joins two components whenever the edge is cheap relative to how much variation each component has already absorbed. That post explains the algorithm over pixels; this one runs the identical merge rule over voxels, with the edge weight replaced by a statistical distance between two Gaussians. I am not re-deriving the merge rule here, it is exactly post 69’s, read there for the threshold formula and why one ascending pass is enough. What is new is the distance function, and what changes when the domain is 3D rather than 2D.

The distance the repo actually computes

My octree extension’s segmentation step weights an edge between two voxels by a Bhattacharyya-style distance between their Gaussians, added in a commit titled “Added HGellinger distance”, my own typo for Hellinger, on 2020-09-27. Trimmed to the arithmetic:

// worldmap.cpp:459-484 (WorldMap::calculateAndSortScores, trimmed)
torch::Tensor meanA = m_loadedChunks[edge.fromChnk]->getFeatMeanAt(edge.fromIdx);
torch::Tensor meanB = m_loadedChunks[edge.toChnk]->getFeatMeanAt(edge.toIdx);
torch::Tensor stdA = m_loadedChunks[edge.fromChnk]->getFeatCovAt(edge.fromIdx);
torch::Tensor stdB = m_loadedChunks[edge.toChnk]->getFeatCovAt(edge.toIdx);
// For diagonal matrix |A| = product(trace(A))
torch::Tensor halfAB = 0.5*(stdA+stdB);
torch::Tensor u = meanA - meanB;
float scale = ( stdA.prod().pow(0.25).item<float>() * stdB.prod().pow(0.25).item<float>() ) /
              ( halfAB.prod().pow(0.5).item<float>() );
float pow = -0.125*(u.square()/halfAB).sum().item<float>();
// TODO: Ensure it is 1.0- and not just exp
edge.score = 1.0 - scale*exp(pow);

That comment, “TODO: Ensure it is 1.0- and not just exp”, is five years old and it turns out to be answerable now: scale*exp(pow) is the Bhattacharyya coefficient $BC$ between two Gaussians, a number between 0 and 1 that says how much two distributions overlap, and 1.0 - scale*exp(pow) is exactly the squared Hellinger distance, $H^2 = 1 - BC$, which is a real distance (symmetric, zero only when the two Gaussians coincide, bounded by 1 when they share nothing). So the TODO’s own question was already correctly answered in the code that shipped; I just never went back and named the thing I’d built. Working through why: for two univariate Gaussians with variances $\sigma_1^2, \sigma_2^2$ and means $\mu_1, \mu_2$, the closed form for that coefficient is

$$BC = \sqrt{\frac{2\sigma_1\sigma_2}{\sigma_1^2+\sigma_2^2}} \exp\!\left(-\frac{(\mu_1-\mu_2)^2}{4(\sigma_1^2+\sigma_2^2)}\right)$$

and multiplying that formula independently across every feature dimension, which is exactly what treating the covariance as diagonal means, turns the square root and the exponent’s sum into the .prod() and .sum() calls above. getFeatCovAt’s own comment says as much:

// worldmap.h:146-148 (Chunk::getFeatCovAt)
const torch::Tensor getFeatCovAt(int i) const {
    // Assumes features take form of diagonal matrix (No dependance between variables and position)
    return m_obs_feat_m2.index({i, torch::indexing::Slice(), SPATIAL_DIMS})/m_obs_count[i];
}

That assumption is fine for the thing it was built for: the “features” are per-point colour or a surface normal, channels that genuinely don’t covary with each other by construction of how they’re measured. It is not fine for the thing this post wants to distinguish, a voxel’s own positional spread, which is exactly correlated whenever the surface it sits on isn’t lined up with the axes. A tilted plane’s covariance has real off-diagonal terms; treating it as diagonal throws away the one signal, orientation, that would tell a floor from a wall.

What changes in 3D: the full covariance, not just its diagonal

Generalising the formula above from “elementwise product and sum over independent channels” to “determinant and matrix inverse over a real covariance matrix” is what a genuinely oriented Gaussian needs. For two multivariate Gaussians with covariances $\Sigma_1, \Sigma_2$ and means $\mu_1, \mu_2$, letting $\Sigma = \tfrac12(\Sigma_1+\Sigma_2)$:

$$BC = \frac{|\Sigma_1|^{1/4}|\Sigma_2|^{1/4}}{|\Sigma|^{1/2}} \exp\!\left(-\frac18(\mu_1-\mu_2)^\top \Sigma^{-1}(\mu_1-\mu_2)\right)$$

Set $\Sigma_1, \Sigma_2$ diagonal and this collapses back to the repo’s own formula exactly, determinant becomes a product, the quadratic form becomes a sum of squares over variance. The determinant and the $\Sigma^{-1}$ applied to a vector are both 3x3, hand-rolled the same way post 65’s getProbabilities port computes its own determinant and Cramer’s-rule inverse, cofactors expanded along one row rather than a general-purpose matrix library:

// src/widgets/point-cloud-segmentation/graph.ts (quadFormSym3, trimmed)
const a11 = yy * zz - yz * yz;
const a12 = xz * yz - xy * zz;
const a13 = xy * yz - yy * xz;
const a22 = xx * zz - xz * xz;
const a23 = xy * xz - xx * yz;
const a33 = xx * yy - xy * xy;
const num = a11*dx*dx + a22*dy*dy + a33*dz*dz + 2*a12*dx*dy + 2*a13*dx*dz + 2*a23*dy*dz;
return num / det;

hellingerDistance in the same file puts the pieces together exactly as the repo’s own formula does, scale * exp(-0.125 * quad) for the Bhattacharyya coefficient, 1 - BC for the distance, just with detSym3 (a full 3x3 determinant) and this cofactor expansion in place of .prod() and elementwise division. A tiny diagonal regulariser (regularise, proportional to the voxel size squared) keeps the determinant off exactly zero for a voxel whose points are almost perfectly coplanar, which this post’s synthetic scene deliberately produces; nothing in the source code needed one, because a feature variance is never that close to degenerate.

A scene built to need this

Showing why the statistical distance earns its keep needs a case where two voxels are close in position but differ in orientation or spread, and the classic one is a plane meeting a plane at an edge. So the widget’s default scene is a small synthetic room corner: a floor, two mutually orthogonal walls meeting it (and each other), and a small tabletop floating clear of all three, each surface sampled with about 1.2 cm of out-of-plane noise standing in for sensor jitter. No real scan is needed to make the point, and it means every number below is exactly reproducible rather than dependent on a particular capture.

With JavaScript on, this becomes a WebGL view of the room-corner scene (or the bundled Stanford Bunny, or your own small PLY/CSV), voxelised and coloured by segment. A voxel-size slider rebuilds the graph; k and minimum-size sliders drive Felzenszwalb & Huttenlocher’s merge rule live; a toggle switches the edge weight between plain centroid distance and the Hellinger distance above, and a readout compares both directly on one floor voxel and one wall voxel right where the two surfaces meet.

The numbers behind the corner

At the scene’s default voxel size (about 15 cm, one-fourteenth of the room’s own extent), 595 of the theoretically addressable cells are occupied by at least 4 points, 2,334 edges connect them. Pick the specific pair the widget’s own readout tracks, one floor voxel and one wall voxel, 18.3 cm apart centre to centre, right where the two surfaces meet:

That last pair of numbers is the whole argument compressed into two figures: centroid distance cannot separate “this pair should merge” from “this pair should not” because both kinds of pair land in overlapping ranges (same-plane 0.54 to 1.61, cross-plane 0.97 to 1.81), but the Hellinger distance keeps 114 cross-plane pairs almost entirely above 0.91 while 1,956 same-plane pairs mostly sit well below that, right up to a maximum of 0.997. Push both metrics through Felzenszwalb & Huttenlocher’s actual merge rule at the same $k = 0.7$ and the difference shows up as contamination, a segment spanning more than one true surface: under centroid distance 15 of the resulting 80 segments already mix two different surfaces, including five of its six largest segments (19, 19, 18, 15 and 14 voxels, each a blend of two different walls or a wall and the floor). Under Hellinger distance only 2 of 59 segments do, and both are small minorities inside an otherwise-correct region, not roughly-even splits: 21 of 24 voxels genuinely wall B, 12 of 21 genuinely floor. Hellinger’s own single largest segment, 30 voxels, is the tabletop, entirely correct, with nothing from any wall or the floor in it at all.

Raise $k$ far enough to make centroid distance consolidate an entire plane into one piece, $k \geq 1.05$ here, and it cannot do that without also fusing the floor, wall A and wall B into a single 530-voxel blob, 89% of the whole scene, because the edge weight that finally lets same-plane voxels merge is the same size as the weight across the corner. The Hellinger distance eventually does the same thing too, past $k \approx 0.86$, collapsing to a 547-voxel blob, 92% of the scene, well before its own theoretical ceiling of 1: the callout below is why.

Drag $k$ yourself, on either metric, and watch this happen: fragments consolidate, then, past a threshold that is far higher and far narrower for the statistical distance, the whole scene collapses into one blob. The widget’s readout reports the live segment count and the same floor/wall comparison at whatever $k$ you land on.

Loading a real scan instead

The bundled Stanford Bunny is offered as a second scene, the same decimated cloud post 64 uses, credited there and here alike. Its curvature changes gradually almost everywhere rather than snapping between orthogonal planes, so there is no sharp orientation change anywhere near as clean as the corner’s for either metric to exploit: raising $k$ on the bunny fragments and then consolidates both metrics without the wide, clean margin the corner shows, and without one metric clearly losing to the other the way centroid distance does at the corner. That is itself the honest finding, not a failure of the widget: the statistical distance’s advantage is largest exactly where two surfaces meet at a sharp, real change of orientation, which a smoothly curved scan mostly doesn’t have.

Why this is plain JavaScript, again

Neither the graph nor the merge loop is a WASM candidate here. A few thousand voxels and a few thousand edges is a sort and a union-find over an array that size, not a per-pixel image operation or dense linear algebra at any real scale; felzenszwalb.rs exists because a camera frame is 50,000 to 150,000 pixels processed 30 times a second, a completely different budget. buildEdges and segmentFH (graph.ts) together take a few milliseconds on this scene, most of it the edge sort, comfortably inside “plain JS is faster to ship” territory, the same call post 65 made for its own 3x3 covariance work.

Where the series lands

Three posts, one structure: a chunked octree that finds any voxel without walking from the root, a running Gaussian in every voxel that never stores a raw point, and now a statistical distance between two of those Gaussians standing in for Felzenszwalb & Huttenlocher’s colour difference. None of the three pieces needed the others to be interesting on their own, but put together they answer a question none of them answers alone: not just “where are the points” or “what does this patch of surface look like”, but “which patches belong to the same surface”, using nothing more than the statistics post 65 already had to compute anyway.