feat(s2): centroids

s2/centroids.ts

@@ -0,0 +1,125 @@
+import { Vector } from '../r3/Vector'
+import { Point } from './Point'
+
+/**
+ * There are several notions of the "centroid" of a triangle. First, there
+ * is the planar centroid, which is simply the centroid of the ordinary
+ * (non-spherical) triangle defined by the three vertices. Second, there is
+ * the surface centroid, which is defined as the intersection of the three
+ * medians of the spherical triangle. It is possible to show that this
+ * point is simply the planar centroid projected to the surface of the
+ * sphere. Finally, there is the true centroid (mass centroid), which is
+ * defined as the surface integral over the spherical triangle of (x,y,z)
+ * divided by the triangle area. This is the point that the triangle would
+ * rotate around if it was spinning in empty space.
+ *
+ * The best centroid for most purposes is the true centroid. Unlike the
+ * planar and surface centroids, the true centroid behaves linearly as
+ * regions are added or subtracted. That is, if you split a triangle into
+ * pieces and compute the average of their centroids (weighted by triangle
+ * area), the result equals the centroid of the original triangle. This is
+ * not true of the other centroids.
+ *
+ * Also note that the surface centroid may be nowhere near the intuitive
+ * "center" of a spherical triangle. For example, consider the triangle
+ * with vertices A=(1,eps,0), B=(0,0,1), C=(-1,eps,0) (a quarter-sphere).
+ * The surface centroid of this triangle is at S=(0, 2*eps, 1), which is
+ * within a distance of 2*eps of the vertex B. Note that the median from A
+ * (the segment connecting A to the midpoint of BC) passes through S, since
+ * this is the shortest path connecting the two endpoints. On the other
+ * hand, the true centroid is at M=(0, 0.5, 0.5), which when projected onto
+ * the surface is a much more reasonable interpretation of the "center" of
+ * this triangle.
+ */
+
+/**
+ * Returns the true centroid of the spherical triangle ABC
+ * multiplied by the signed area of spherical triangle ABC. The reasons for
+ * multiplying by the signed area are (1) this is the quantity that needs to be
+ * summed to compute the centroid of a union or difference of triangles, and
+ * (2) it's actually easier to calculate this way. All points must have unit length.
+ *
+ * Note that the result of this function is defined to be Point(0, 0, 0) if
+ * the triangle is degenerate.
+ */
+export const trueCentroid = (a: Point, b: Point, c: Point): Point => {
+  // Use Distance to get accurate results for small triangles.
+  let ra = 1
+  const sa = b.distance(c)
+  if (sa !== 0) ra = sa / Math.sin(sa)
+
+  let rb = 1
+  const sb = c.distance(a)
+  if (sb !== 0) rb = sb / Math.sin(sb)
+
+  let rc = 1
+  const sc = a.distance(b)
+  if (sc !== 0) rc = sc / Math.sin(sc)
+
+  // Now compute a point M such that:
+  //
+  //  [Ax Ay Az] [Mx]                       [ra]
+  //  [Bx By Bz] [My]  = 0.5 * det(A,B,C) * [rb]
+  //  [Cx Cy Cz] [Mz]                       [rc]
+  //
+  // To improve the numerical stability we subtract the first row (A) from the
+  // other two rows; this reduces the cancellation error when A, B, and C are
+  // very close together. Then we solve it using Cramer's rule.
+  //
+  // The result is the true centroid of the triangle multiplied by the
+  // triangle's area.
+  //
+  // This code still isn't as numerically stable as it could be.
+  // The biggest potential improvement is to compute B-A and C-A more
+  // accurately so that (B-A)x(C-A) is always inside triangle ABC.
+  const x = new Vector(a.x, b.x - a.x, c.x - a.x)
+  const y = new Vector(a.y, b.y - a.y, c.y - a.y)
+  const z = new Vector(a.z, b.z - a.z, c.z - a.z)
+  const r = new Vector(ra, rb - ra, rc - ra)
+
+  return Point.fromVector(new Vector(y.cross(z).dot(r), z.cross(x).dot(r), x.cross(y).dot(r)).mul(0.5))
+}
+
+/**
+ * Returns the true centroid of the spherical geodesic edge AB
+ * multiplied by the length of the edge AB. As with triangles, the true centroid
+ * of a collection of line segments may be computed simply by summing the result
+ * of this method for each segment.
+ *
+ * Note that the planar centroid of a line segment is simply 0.5 * (a + b),
+ * while the surface centroid is (a + b).Normalize(). However neither of
+ * these values is appropriate for computing the centroid of a collection of
+ * edges (such as a polyline).
+ *
+ * Also note that the result of this function is defined to be Point(0, 0, 0)
+ * if the edge is degenerate.
+ */
+export const edgeTrueCentroid = (a: Point, b: Point): Point => {
+  // The centroid (multiplied by length) is a vector toward the midpoint
+  // of the edge, whose length is twice the sine of half the angle between
+  // the two vertices. Defining theta to be this angle, we have:
+  const vDiff = a.vector.sub(b.vector) // Length == 2*sin(theta)
+  const vSum = a.vector.add(b.vector) // Length == 2*cos(theta)
+  const sin2 = vDiff.norm2()
+  const cos2 = vSum.norm2()
+
+  if (cos2 === 0) return new Point(0, 0, 0) // Ignore antipodal edges.
+
+  return Point.fromVector(vSum.mul(Math.sqrt(sin2 / cos2))) // Length == 2*sin(theta)
+}
+
+/**
+ * Returns the centroid of the planar triangle ABC. This can be
+ * normalized to unit length to obtain the "surface centroid" of the corresponding
+ * spherical triangle, i.e. the intersection of the three medians. However, note
+ * that for large spherical triangles the surface centroid may be nowhere near
+ * the intuitive "center".
+ */
+export const planarCentroid = (a: Point, b: Point, c: Point): Point => {
+  return Point.fromVector(
+    a.vector
+      .add(b.vector)
+      .add(c.vector)
+      .mul(1 / 3)
+  )
+}

s2/centroids_test.ts

@@ -0,0 +1,91 @@
+import { test, describe } from 'node:test'
+import { equal, ok } from 'node:assert/strict'
+import { Point } from './Point'
+import { edgeTrueCentroid, planarCentroid, trueCentroid } from './centroids'
+import { randomFrame, randomFrameAtPoint, randomPoint } from './testing'
+import * as matrix from './matrix3x3'
+
+describe('s2.centroids', () => {
+  test('planarCentroid', () => {
+    const tests = [
+      {
+        name: 'xyz axis',
+        p0: new Point(0, 0, 1),
+        p1: new Point(0, 1, 0),
+        p2: new Point(1, 0, 0),
+        want: new Point(1 / 3, 1 / 3, 1 / 3)
+      },
+      {
+        name: 'Same point',
+        p0: new Point(1, 0, 0),
+        p1: new Point(1, 0, 0),
+        p2: new Point(1, 0, 0),
+        want: new Point(1, 0, 0)
+      }
+    ]
+
+    for (const { name, p0, p1, p2, want } of tests) {
+      const got = planarCentroid(p0, p1, p2)
+      ok(got.approxEqual(want), `${name}: PlanarCentroid(${p0}, ${p1}, ${p2}) = ${got}, want ${want}`)
+    }
+  })
+
+  test('trueCentroid', () => {
+    for (let i = 0; i < 100; i++) {
+      const f = randomFrame()
+      const p = matrix.col(f, 0)
+      const x = matrix.col(f, 1)
+      const y = matrix.col(f, 2)
+      const d = 1e-4 * Math.pow(1e-4, Math.random())
+
+      let p0 = Point.fromVector(p.vector.sub(x.vector.mul(d)).normalize())
+      let p1 = Point.fromVector(p.vector.add(x.vector.mul(d)).normalize())
+      let p2 = Point.fromVector(p.vector.add(y.vector.mul(d * 3)).normalize())
+      let want = Point.fromVector(p.vector.add(y.vector.mul(d)).normalize())
+
+      let got = trueCentroid(p0, p1, p2).vector.normalize()
+      ok(got.distance(want.vector) < 2e-8, `TrueCentroid(${p0}, ${p1}, ${p2}).Normalize() = ${got}, want ${want}`)
+
+      p0 = Point.fromVector(p.vector)
+      p1 = Point.fromVector(p.vector.add(x.vector.mul(d * 3)).normalize())
+      p2 = Point.fromVector(p.vector.add(y.vector.mul(d * 6)).normalize())
+      want = Point.fromVector(p.vector.add(x.vector.add(y.vector.mul(2)).mul(d)).normalize())
+
+      got = trueCentroid(p0, p1, p2).vector.normalize()
+      ok(got.distance(want.vector) < 2e-8, `TrueCentroid(${p0}, ${p1}, ${p2}).Normalize() = ${got}, want ${want}`)
+    }
+  })
+
+  test('edgeTrueCentroid semi-circles', () => {
+    const a = Point.fromCoords(0, -1, 0)
+    const b = Point.fromCoords(1, 0, 0)
+    const c = Point.fromCoords(0, 1, 0)
+    const centroid = Point.fromVector(edgeTrueCentroid(a, b).vector.add(edgeTrueCentroid(b, c).vector))
+
+    ok(
+      b.approxEqual(Point.fromVector(centroid.vector.normalize())),
+      `EdgeTrueCentroid(${a}, ${b}) + EdgeTrueCentroid(${b}, ${c}) = ${centroid}, want ${b}`
+    )
+    equal(centroid.vector.norm(), 2.0, `${centroid}.Norm() = ${centroid.vector.norm()}, want 2.0`)
+  })
+
+  test('edgeTrueCentroid great-circles', () => {
+    for (let iter = 0; iter < 100; iter++) {
+      const f = randomFrameAtPoint(randomPoint())
+      const x = matrix.col(f, 0)
+      const y = matrix.col(f, 1)
+
+      let centroid = new Point(0, 0, 0)
+
+      let v0 = x
+      for (let theta = 0.0; theta < 2 * Math.PI; theta += Math.pow(Math.random(), 10)) {
+        const v1 = Point.fromVector(x.vector.mul(Math.cos(theta)).add(y.vector.mul(Math.sin(theta))))
+        centroid = Point.fromVector(centroid.vector.add(edgeTrueCentroid(v0, v1).vector))
+        v0 = v1
+      }
+
+      centroid = Point.fromVector(centroid.vector.add(edgeTrueCentroid(v0, x).vector))
+      ok(centroid.vector.norm() <= 2e-14, `${centroid}.Norm() = ${centroid.vector.norm()}, want <= 2e-14`)
+    }
+  })
+})