GEOMETRY-71 Fix spherical barycenter algorithm
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
index d87dceb..c187cbb 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
@@ -128,29 +128,35 @@
/** {@inheritDoc} */
@Override
public Point2S getBarycenter() {
- List<GreatArc> arcs = getBoundaries();
- int numSides = arcs.size();
+ Vector3D weighted = getWeightedBarycenter();
+ return weighted == null ? null : Point2S.from(weighted);
+ }
- if (numSides == 0) {
+ /**
+ * Returns the weighted vector for barycenter.
+ *
+ * @return Weighted barycenter
+ */
+ Vector3D getWeightedBarycenter() {
+ List<GreatArc> arcs = getBoundaries();
+ switch (arcs.size()) {
+ case 0:
// full space; no barycenter
return null;
- } else if (numSides == 1) {
+ case 1:
// hemisphere; barycenter is the pole of the hemisphere
- return arcs.get(0).getCircle().getPolePoint();
- } else {
+ GreatArc singleArc = arcs.get(0);
+ return singleArc.getCircle().getPole().withNorm(singleArc.getSize());
+ default:
// 2 or more sides; use an extension of the approach outlined here:
// https://archive.org/details/centroidinertiat00broc
// In short, the barycenter is the sum of the pole vectors of each side
// multiplied by their arc lengths.
Vector3D barycenter = Vector3D.ZERO;
-
for (GreatArc arc : getBoundaries()) {
- barycenter = Vector3D.linearCombination(
- 1, barycenter,
- arc.getSize(), arc.getCircle().getPole());
+ barycenter = barycenter.add(arc.getCircle().getPole().withNorm(arc.getSize()));
}
-
- return Point2S.from(barycenter);
+ return barycenter;
}
}
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
index 8c721aa..17ddd45 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
@@ -169,12 +169,8 @@
double size;
for (ConvexArea2S area : areas) {
- size = area.getSize();
-
- sizeSum += size;
- barycenterVector = Vector3D.linearCombination(
- 1, barycenterVector,
- size, area.getBarycenter().getVector());
+ sizeSum += area.getSize();
+ barycenterVector = barycenterVector.add(area.getWeightedBarycenter());
}
Point2S barycenter = barycenterVector.eq(Vector3D.ZERO, precision) ?
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
index 7a7943d..3e0e290 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
@@ -706,6 +706,54 @@
checkBarycenterConsistency(france);
}
+ @Test
+ public void testCircleToPolygonBarycenter() {
+ double radius = 0.0001;
+ Point2S center = Point2S.of(1.0, 1.0);
+ int numPts = 200;
+
+ // counterclockwise
+ RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
+ SphericalTestUtils.assertPointsEq(center, ccw.getBarycenter(), BARYCENTER_EPS);
+
+ // clockwise; barycenter should just be antipodal for the circle center
+ RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+ SphericalTestUtils.assertPointsEq(center.antipodal(), cw.getBarycenter(), BARYCENTER_EPS);
+ }
+
+ @Test
+ public void testCircleToPolygonSize() {
+ double radius = 0.0001;
+ Point2S center = Point2S.of(1.0, 1.0);
+ int numPts = 200;
+
+ // https://en.wikipedia.org/wiki/Spherical_cap
+ double ccwArea = 4.0 * PlaneAngleRadians.PI * Math.pow(Math.sin(radius / 2.0), 2.0);
+ double cwArea = 4.0 * PlaneAngleRadians.PI - ccwArea;
+
+ RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
+ Assert.assertEquals("Counterclockwise size", ccwArea, ccw.getSize(), TEST_EPS);
+
+ RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+ Assert.assertEquals("Clockwise size", cwArea, cw.getSize(), TEST_EPS);
+ }
+
+ @Test
+ public void testCircleToPolygonBoundarySize() {
+ double radius = 0.0001;
+ Point2S center = Point2S.of(1.0, 1.0);
+ int numPts = 200;
+
+ // boundary size is independent from winding
+ double boundary = PlaneAngleRadians.TWO_PI * Math.sin(radius);
+
+ RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
+ Assert.assertEquals("Counterclockwise boundary size", boundary, ccw.getBoundarySize(), 1.0e-7);
+
+ RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+ Assert.assertEquals("Clockwise boundary size", boundary, cw.getBoundarySize(), 1.0e-7);
+ }
+
/**
* Insert convex subhyperplanes defining the positive quadrant area.
* @param tree
@@ -791,4 +839,22 @@
SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, BARYCENTER_EPS);
}
}
+
+ private static RegionBSPTree2S circleToPolygon(Point2S center, double radius, int numPts, boolean clockwise) {
+ List<Point2S> pts = new ArrayList<>(numPts);
+
+ // get an arbitrary point on the circle boundary
+ pts.add(Transform2S.createRotation(center.getVector().orthogonal(), radius).apply(center));
+
+ // create the list of boundary points by rotating the previous point around the circle center
+ double span = PlaneAngleRadians.TWO_PI / numPts;
+
+ // negate the span for clockwise winding
+ Transform2S rotate = Transform2S.createRotation(center, clockwise ? -span : span);
+ for (int i = 1; i < numPts; ++i) {
+ pts.add(rotate.apply(pts.get(i - 1)));
+ }
+
+ return GreatArcPath.fromVertexLoop(pts, TEST_PRECISION).toTree();
+ }
}