[geometry] Improve CalcOrientationByPca() for covariance matrix with repeated eigenvalues.

[DamrongGuoy]

Introduction

In PR #13808, we use Principal Component Analysis (PCA) to calculate initial orientation of an oriented bounding box (OBB) of a point set. However, when the covariance matrix has repeated eigenvalues, there are infinite solutions of PCA, and ObbMaker::CalcOrientationByPca() does not always pick a good solution. These two examples illustrate the problem:

• The set of vertices of a regular octahedron has its covariance matrix with eigenvalue of multiplicity 3. Any three mutually orthonormal vectors are valid solutions to PCA. However, this picture shows that one choice of the basis gives a much better OBB than another choice of basis.

• Four vertices of a square embedded in ℝ³. Their covariance matrix has one eigenvalue of multiplicity 2 (and the one remaining eigenvalue is zero). Any two mutually orthonormal vectors spanning the plane of the square are valid for the first two principal components. However, only the basis that is aligned with the square is optimal.

Analysis

For an eigenvalue with multiplicity m, SelfAdjointEigenSolver will pick an orthonormal basis to span the m-dimensional subspace arbitrarily. For m=2, it can rotate two basis vectors around a plane arbitrarily keeping the remaining basis vector fixed. For m=3, it can rotate the three basis vectors arbitrarily to anywhere in 3D rotation group (SO(3)). The arbitrary choice is deterministic and depends on the actual Eigen implementation that we have no control.

Let a ≥ b ≥ c be the three possibly repeated eigenvalues of the covariance matrix of an input point set. For a co-planar point set, there is only one case to consider:

• a = b > 0 with c = 0. For example, any uniform sampling of a circle (here uniformity means every pair of two consecutive points have the same distance) including four vertices of a square or three vertices of an equilateral triangle.

For a true 3D point set (its convex hull has non-zero volume), there are three cases of repeated eigenvalues to consider:

1. a = b = c. For example, a regular octahedron. Intuitively the octahedron has no preferred direction under PCA.
2. a > b = c. For example, stretch a regular octahedron to be longer along a diagonal (direction of a).
3. a = b > c. For example, shrink a regular octahedron to be shorter along a diagonal (direction of c).

Notice that for a collinear point set, a > b = c = 0. SelfAdjointEigenSolver can rotate the two basis vectors that are perpendicular to the first principal component arbitrarily. Luckily it does not concern us because all OBBs of those bases have the same volume.

Available unit tests

We used six vertices of a standard octahedron in the unit tests in geometry/proximity/test/obb_test.cc:

• TEST_F(ObbMakerTestOctahedron, CalcOrientationByPca)
• TEST_F(ObbMakerTestOctahedron, ObbMakerCompute)

The first test showed that the returned principal components were not stable under rigid transform. Due to the triple multiplicity, a small numerical perturbation (due to rigid transform) triggered an arbitrary (but deterministic) basis in SO(3).

The second test showed that two arbitrary bases from PCA (due to rigid transform) gave very different oriented bounding boxes with volume difference about 67%.

[SeanCurtis-TRI]

Several thoughts in passing:

1. Ultimately, this is in service to producing tighter-fitting Obbs. Currently, we use PCA to create an initial guess and then apply a refinement step to improve that guess. There are two levers to pull on: a better guess or better refinement. This focuses solely on the first part (#14081 focuses on the second). It may be that better refinement renders the better guess moot.
2. While the issue description provides a clear characterization of the problem, it doesn't suggest any possible solutions.
3. One possible solution that has been mooted is to supplant PCA analysis with inertia tensor calculation. It works as follows:
   • The approach
     • Compute an inertia tensor for each element (triangle or tet).
     • As we partition elements to build the hierarchy, compute the combined inertia tensor for all elements included in the domain.
     • Use the principal axes of rotation as the initial guess for the basis for the Obb.
     • Refine from that.
   • Issues
     • For symmetric elements, the inertia tensor will suffer from the same issues as PCA in terms of being arbitrarily rotated. One can only hope that in aggregate the quality of the guess will be higher quality.
     • This is definitely more expensive than doing the PCA.
     • It isn't clear that this provides a strictly better initial guess.