bestfit_plane_numpy returns wrong normals on near-collinear point clouds (covariance step squares the conditioning)

[mengxihex]

Describe the bug

bestfit_plane_numpy (via pca_numpy) fits by building the covariance matrix and SVD-ing it:

# src/compas/geometry/pca_numpy.py
C = Y.T.dot(Y) / (n - 1)
u, s, vT = svd(C, full_matrices=False)

Forming Y.T @ Y squares the condition number: a point cloud whose second in-plane extent is ~1e-8 of its first (a thin sliver — still a perfectly well-defined plane) has a covariance eigenvalue ratio of ~1e-16, at machine epsilon, so the smallest principal direction is lost to rounding and the returned plane normal tips into the sliver plane. An SVD of the centered data matrix Y itself (whose right singular vectors are mathematically the same) only has to resolve a ~1e-8 singular-value ratio and recovers the normal to full precision.

Slivers like this are routine in AEC data — footprints of thin panels, wall edges, extracted face outlines — which is where we hit it (migrating a timber-engineering geometry engine onto COMPAS; our legacy data-matrix fit and bestfit_plane_numpy disagreed by up to 1.31 per normal component on real inputs).

To Reproduce

import numpy as np
from compas.geometry import bestfit_plane_numpy

rng = np.random.default_rng(7)

# A thin sliver in a known plane: long axis ~1.0, short in-plane axis ~1e-8.
n = 200
t = rng.uniform(-1.0, 1.0, n)          # long in-plane direction
s = rng.uniform(-1e-8, 1e-8, n)        # short in-plane direction (near-collinear)
pts_local = np.column_stack([t, s, np.zeros(n)])  # exactly planar, normal = +Z

# Rotate into a generic orientation so nothing is axis-aligned.
a, b = 0.6, -1.1
Rx = np.array([[1, 0, 0], [0, np.cos(a), -np.sin(a)], [0, np.sin(a), np.cos(a)]])
Rz = np.array([[np.cos(b), -np.sin(b), 0], [np.sin(b), np.cos(b), 0], [0, 0, 1]])
R = Rz @ Rx
pts = pts_local @ R.T + np.array([100.0, -40.0, 7.0])
true_normal = R @ np.array([0.0, 0.0, 1.0])

_, normal = bestfit_plane_numpy(pts)
normal = np.asarray(normal) / np.linalg.norm(normal)

Y = pts - pts.mean(axis=0)
normal_svd = np.linalg.svd(Y, full_matrices=False)[2][2]

for name, nrm in (("bestfit_plane_numpy", normal), ("data-matrix SVD", normal_svd)):
    ang = np.degrees(np.arccos(np.clip(abs(nrm @ true_normal), -1, 1)))
    print(f"{name:20s}: angle off = {ang:.6f} deg")

Output (compas 2.15.1, numpy 2.4.6, CPython 3.13, macOS arm64):

bestfit_plane_numpy : angle off = 11.638995 deg
data-matrix SVD     : angle off = 0.000001 deg

Shrinking the short axis to 1e-9 makes it worse (37.28 deg vs 0.000004 deg). The points are exactly planar in both cases — the fit problem is perfectly posed; only the covariance construction loses it.

Confirmed present on current main as well as 2.15.1 (the C = Y.T.dot(Y) / (n - 1) construction is unchanged).

Expected behavior

The normal of an exactly-planar (if skinny) point cloud should be recovered — the direct SVD of the centered coordinates does so to machine precision.

Suggested fix

In pca_numpy, run the SVD on the centered data matrix directly:

u, s, vT = svd(Y, full_matrices=False)   # right singular vectors == eigenvectors of C
eigenvalues = (s ** 2) / (n - 1)          # if the current return contract needs them

Right singular vectors of Y are identical (in exact arithmetic) to the eigenvectors of C, so callers see the same results on well-conditioned inputs — only the ill-conditioned ones improve. Cost-wise both routes are O(n) once the Y.T @ Y product is counted; for typical point-cloud sizes the difference is noise.

Happy to open a PR with the change + tests along these lines if you agree with the direction.