This is an automated email from the ASF dual-hosted git repository. erans pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-geometry.git
commit b1cf35815c01a861172eead1a498a7bb4b5ff69f Author: Matt Juntunen <[email protected]> AuthorDate: Sat Jan 18 08:42:12 2020 -0500 GEOMETRY-71: adding additional unit tests and documentation regarding the computation of spherical barycenters --- .../org/apache/commons/geometry/core/Region.java | 4 +- .../geometry/spherical/twod/ConvexArea2S.java | 17 +- .../geometry/spherical/twod/RegionBSPTree2S.java | 5 +- .../geometry/spherical/twod/ConvexArea2STest.java | 33 ++-- .../spherical/twod/RegionBSPTree2STest.java | 195 ++++++++++++++++++--- 5 files changed, 202 insertions(+), 52 deletions(-) diff --git a/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java b/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java index 233cd4b..2c15111 100644 --- a/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java +++ b/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/Region.java @@ -49,8 +49,8 @@ public interface Region<P extends Point<P>> { */ double getBoundarySize(); - /** Get the barycenter of the region or null if none exists. A barycenter - * will not exist for empty or infinite regions. + /** Get the barycenter of the region or null if no single, unique barycenter exists. + * A barycenter will not exist for empty or infinite regions. * @return the barycenter of the region or null if none exists */ P getBarycenter(); 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 c187cbb..4c56473 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 @@ -24,7 +24,6 @@ import java.util.List; import java.util.stream.Collectors; import java.util.stream.Stream; -import org.apache.commons.numbers.angle.PlaneAngleRadians; import org.apache.commons.geometry.core.Transform; import org.apache.commons.geometry.core.partitioning.AbstractConvexHyperplaneBoundedRegion; import org.apache.commons.geometry.core.partitioning.ConvexSubHyperplane; @@ -32,6 +31,7 @@ import org.apache.commons.geometry.core.partitioning.Hyperplane; import org.apache.commons.geometry.core.partitioning.Split; import org.apache.commons.geometry.core.precision.DoublePrecisionContext; import org.apache.commons.geometry.euclidean.threed.Vector3D; +import org.apache.commons.numbers.angle.PlaneAngleRadians; /** Class representing a convex area in 2D spherical space. The boundaries of this * area, if any, are composed of convex great circle arcs. @@ -128,16 +128,19 @@ public final class ConvexArea2S extends AbstractConvexHyperplaneBoundedRegion<Po /** {@inheritDoc} */ @Override public Point2S getBarycenter() { - Vector3D weighted = getWeightedBarycenter(); + Vector3D weighted = getWeightedBarycenterVector(); return weighted == null ? null : Point2S.from(weighted); } - /** - * Returns the weighted vector for barycenter. - * - * @return Weighted barycenter + /** Returns the weighted vector for the barycenter. This vector is computed by scaling the + * pole vector of the great circle of each boundary arc by the size of the arc and summing + * the results. By combining the weighted barycenter vectors of multiple areas, a single + * barycenter can be computed for the whole group. + * @return weighted barycenter vector. + * @see <a href="https://archive.org/details/centroidinertiat00broc";> + * <em>The Centroid and Inertia Tensor for a Spherical Triangle</em> - John E. Brock</a> */ - Vector3D getWeightedBarycenter() { + Vector3D getWeightedBarycenterVector() { List<GreatArc> arcs = getBoundaries(); switch (arcs.size()) { case 0: 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 17ddd45..c79d894 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 @@ -167,13 +167,12 @@ public class RegionBSPTree2S extends AbstractRegionBSPTree<Point2S, RegionBSPTre double sizeSum = 0; Vector3D barycenterVector = Vector3D.ZERO; - double size; for (ConvexArea2S area : areas) { sizeSum += area.getSize(); - barycenterVector = barycenterVector.add(area.getWeightedBarycenter()); + barycenterVector = barycenterVector.add(area.getWeightedBarycenterVector()); } - Point2S barycenter = barycenterVector.eq(Vector3D.ZERO, precision) ? + final Point2S barycenter = barycenterVector.eq(Vector3D.ZERO, precision) ? null : Point2S.from(barycenterVector); diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java index 0307804..759e9db 100644 --- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java +++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java @@ -42,10 +42,6 @@ public class ConvexArea2STest { private static final DoublePrecisionContext TEST_PRECISION = new EpsilonDoublePrecisionContext(TEST_EPS); - // epsilon value for use when comparing computed barycenter locations; - // this must currently be set much higher than the other epsilon - private static final double BARYCENTER_EPS = 1e-2; - @Test public void testFull() { // act @@ -257,7 +253,7 @@ public class ConvexArea2STest { Assert.assertEquals(size, area.getSize(), TEST_EPS); Point2S expectedBarycenter = triangleBarycenter(p1, p2, p3); - SphericalTestUtils.assertPointsEq(expectedBarycenter, area.getBarycenter(), BARYCENTER_EPS); + SphericalTestUtils.assertPointsEq(expectedBarycenter, area.getBarycenter(), TEST_EPS); checkBarycenterConsistency(area); @@ -759,11 +755,18 @@ public class ConvexArea2STest { } private static Point2S triangleBarycenter(Point2S p1, Point2S p2, Point2S p3) { - // compute the barycenter using intersection mid point arcs - GreatCircle c1 = GreatCircle.fromPoints(p1, p2.slerp(p3, 0.5), TEST_PRECISION); - GreatCircle c2 = GreatCircle.fromPoints(p2, p1.slerp(p3, 0.5), TEST_PRECISION); + // compute the barycenter as the sum of the cross product of each point pair weighted by + // the angle between the points + Vector3D v1 = p1.getVector(); + Vector3D v2 = p2.getVector(); + Vector3D v3 = p3.getVector(); - return c1.intersection(c2); + Vector3D sum = Vector3D.ZERO; + sum = sum.add(v1.cross(v2).withNorm(v1.angle(v2))); + sum = sum.add(v2.cross(v3).withNorm(v2.angle(v3))); + sum = sum.add(v3.cross(v1).withNorm(v3.angle(v1))); + + return Point2S.from(sum); } private static void checkArc(GreatArc arc, Point2S start, Point2S end) { @@ -802,21 +805,15 @@ public class ConvexArea2STest { ConvexArea2S minus = split.getMinus(); double minusSize = minus.getSize(); - Point2S minusBc = minus.getBarycenter(); - - Vector3D weightedMinus = minusBc.getVector() - .multiply(minus.getSize()); ConvexArea2S plus = split.getPlus(); double plusSize = plus.getSize(); - Point2S plusBc = plus.getBarycenter(); - Vector3D weightedPlus = plusBc.getVector() - .multiply(plus.getSize()); - Point2S computedBarycenter = Point2S.from(weightedMinus.add(weightedPlus)); + Point2S computedBarycenter = Point2S.from(minus.getWeightedBarycenterVector() + .add(plus.getWeightedBarycenterVector())); Assert.assertEquals(size, minusSize + plusSize, TEST_EPS); - SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, BARYCENTER_EPS); + SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, TEST_EPS); } } } 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 3e0e290..a7040e9 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 @@ -39,13 +39,13 @@ public class RegionBSPTree2STest { private static final double TEST_EPS = 1e-10; + // alternative epsilon value for checking the barycenters of complex + // or very small regions + private static final double BARYCENTER_EPS = 1e-5; + private static final DoublePrecisionContext TEST_PRECISION = new EpsilonDoublePrecisionContext(TEST_EPS); - // epsilon value for use when comparing computed barycenter locations; - // this must currently be set much higher than the other epsilon - private static final double BARYCENTER_EPS = 1e-2; - private static final GreatCircle EQUATOR = GreatCircle.fromPoleAndU( Vector3D.Unit.PLUS_Z, Vector3D.Unit.PLUS_X, TEST_PRECISION); @@ -339,7 +339,7 @@ public class RegionBSPTree2STest { } @Test - public void testToConvex_doubleLune_comlement() { + public void testToConvex_doubleLune_complement() { // arrange RegionBSPTree2S tree = GreatArcPath.builder(TEST_PRECISION) .append(EQUATOR.arc(0, PlaneAngleRadians.PI)) @@ -523,12 +523,12 @@ public class RegionBSPTree2STest { Assert.assertEquals(3.5 * PlaneAngleRadians.PI, tree.getSize(), TEST_EPS); Assert.assertEquals(1.5 * PlaneAngleRadians.PI, tree.getBoundarySize(), TEST_EPS); -// Point2S center = Point2S.from(Point2S.MINUS_K.getVector() -// .add(Point2S.PLUS_I.getVector()) -// .add(Point2S.MINUS_J.getVector())); -// SphericalTestUtils.assertPointsEq(center.antipodal(), tree.getBarycenter(), TEST_EPS); -// -// checkBarycenterConsistency(tree); + Point2S center = Point2S.from(Point2S.MINUS_K.getVector() + .add(Point2S.PLUS_I.getVector()) + .add(Point2S.MINUS_J.getVector())); + SphericalTestUtils.assertPointsEq(center.antipodal(), tree.getBarycenter(), TEST_EPS); + + checkBarycenterConsistency(tree); List<GreatArcPath> paths = tree.getBoundaryPaths(); Assert.assertEquals(1, paths.size()); @@ -543,6 +543,105 @@ public class RegionBSPTree2STest { } @Test + public void testGeometricProperties_polygonWithHole() { + // arrange + Point2S center = Point2S.of(0.5, 2); + + double outerRadius = 1; + double innerRadius = 0.5; + + RegionBSPTree2S outer = buildDiamond(center, outerRadius); + RegionBSPTree2S inner = buildDiamond(center, innerRadius); + + // rotate the inner diamond a quarter turn to become a square + inner.transform(Transform2S.createRotation(center, 0.25 * Math.PI)); + + // act + RegionBSPTree2S tree = RegionBSPTree2S.empty(); + tree.difference(outer, inner); + + // assert + double area = 4 * (rightTriangleArea(outerRadius, outerRadius) - rightTriangleArea(innerRadius, innerRadius)); + Assert.assertEquals(area, tree.getSize(), TEST_EPS); + + double outerSideLength = sphericalHypot(outerRadius, outerRadius); + double innerSideLength = sphericalHypot(innerRadius, innerRadius); + double boundarySize = 4 * (outerSideLength + innerSideLength); + Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS); + + SphericalTestUtils.assertPointsEq(center, tree.getBarycenter(), TEST_EPS); + checkBarycenterConsistency(tree); + + SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE, center); + } + + @Test + public void testGeometricProperties_polygonWithHole_complex() { + // arrange + Point2S center = Point2S.of(0.5, 2); + + double outerRadius = 2; + double midRadius = 1; + double innerRadius = 0.5; + + RegionBSPTree2S outer = buildDiamond(center, outerRadius); + RegionBSPTree2S mid = buildDiamond(center, midRadius); + RegionBSPTree2S inner = buildDiamond(center, innerRadius); + + // rotate the middle diamond a quarter turn to become a square + mid.transform(Transform2S.createRotation(center, 0.25 * Math.PI)); + + // act + RegionBSPTree2S tree = RegionBSPTree2S.empty(); + tree.difference(outer, mid); + tree.union(inner); + tree.complement(); + + // assert + // compute the area, adjusting the first computation for the fact that the triangles comprising the + // outer diamond have lengths greater than pi/2 + double nonComplementedArea = 4 * ((PlaneAngleRadians.PI - rightTriangleArea(outerRadius, outerRadius) - + rightTriangleArea(midRadius, midRadius) + rightTriangleArea(innerRadius, innerRadius))); + double area = (4 * PlaneAngleRadians.PI) - nonComplementedArea; + Assert.assertEquals(area, tree.getSize(), TEST_EPS); + + double outerSideLength = sphericalHypot(outerRadius, outerRadius); + double midSideLength = sphericalHypot(midRadius, midRadius); + double innerSideLength = sphericalHypot(innerRadius, innerRadius); + double boundarySize = 4 * (outerSideLength + midSideLength + innerSideLength); + Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS); + + SphericalTestUtils.assertPointsEq(center.antipodal(), tree.getBarycenter(), TEST_EPS); + checkBarycenterConsistency(tree); + + SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE, center); + } + + @Test + public void testGeometricProperties_equalAndOppositeRegions() { + // arrange + Point2S center = Point2S.PLUS_I; + double radius = 0.25 * Math.PI; + + RegionBSPTree2S a = buildDiamond(center, radius); + RegionBSPTree2S b = buildDiamond(center.antipodal(), radius); + + // act + RegionBSPTree2S tree = RegionBSPTree2S.empty(); + tree.union(a, b); + + // assert + double area = 8 * rightTriangleArea(radius, radius); + Assert.assertEquals(area, tree.getSize(), TEST_EPS); + + double boundarySize = 8 * sphericalHypot(radius, radius); + Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS); + + // should be null since no unique barycenter exists + Assert.assertNull(tree.getBarycenter()); + } + + @Test public void testSplit_both() { // arrange GreatCircle c1 = GreatCircle.fromPole(Vector3D.Unit.MINUS_X, TEST_PRECISION); @@ -700,6 +799,7 @@ public class RegionBSPTree2STest { // assert Assert.assertEquals(0.6316801448267251, france.getBoundarySize(), TEST_EPS); Assert.assertEquals(0.013964220234478741, france.getSize(), TEST_EPS); + SphericalTestUtils.assertPointsEq(Point2S.of(0.04368552749392928, 0.7590839905197961), france.getBarycenter(), BARYCENTER_EPS); @@ -809,8 +909,6 @@ public class RegionBSPTree2STest { Point2S barycenter = region.getBarycenter(); double size = region.getSize(); - SphericalTestUtils.checkClassify(region, RegionLocation.INSIDE, barycenter); - GreatCircle circle = GreatCircle.fromPole(barycenter.getVector(), TEST_PRECISION); for (double az = 0; az <= PlaneAngleRadians.TWO_PI; az += 0.2) { Point2S pt = circle.toSpace(Point1S.of(az)); @@ -822,24 +920,77 @@ public class RegionBSPTree2STest { RegionBSPTree2S minus = split.getMinus(); double minusSize = minus.getSize(); - Point2S minusBc = minus.getBarycenter(); - - Vector3D weightedMinus = minusBc.getVector() - .multiply(minus.getSize()); RegionBSPTree2S plus = split.getPlus(); double plusSize = plus.getSize(); - Point2S plusBc = plus.getBarycenter(); - Vector3D weightedPlus = plusBc.getVector() - .multiply(plus.getSize()); - Point2S computedBarycenter = Point2S.from(weightedMinus.add(weightedPlus)); + Point2S computedBarycenter = Point2S.from(weightedBarycenterVector(minus) + .add(weightedBarycenterVector(plus))); Assert.assertEquals(size, minusSize + plusSize, TEST_EPS); - SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, BARYCENTER_EPS); + SphericalTestUtils.assertPointsEq(barycenter, computedBarycenter, TEST_EPS); } } + private static Vector3D weightedBarycenterVector(RegionBSPTree2S tree) { + Vector3D sum = Vector3D.ZERO; + for (ConvexArea2S convex : tree.toConvex()) { + sum = sum.add(convex.getWeightedBarycenterVector()); + } + + return sum; + } + + private static RegionBSPTree2S buildDiamond(Point2S center, double radius) { + Vector3D u = center.getVector(); + Vector3D w = u.orthogonal(Vector3D.Unit.PLUS_Z); + Vector3D v = w.cross(u); + + Transform2S rotV = Transform2S.createRotation(v, radius); + Transform2S rotW = Transform2S.createRotation(w, radius); + + Point2S top = rotV.inverse().apply(center); + Point2S bottom = rotV.apply(center); + + Point2S right = rotW.apply(center); + Point2S left = rotW.inverse().apply(center); + + return GreatArcPath.fromVertexLoop(Arrays.asList(top, left, bottom, right), TEST_PRECISION) + .toTree(); + } + + /** Solve for the hypotenuse of a spherical right triangle, given the lengths of the + * other two side. The sides must have lengths less than pi/2. + * @param a first side; must be less than pi/2 + * @param b second side; must be less than pi/2 + * @return the hypotenuse of the spherical right triangle with sides of the given lengths + */ + private static double sphericalHypot(double a, double b) { + // use the spherical law of cosines and the fact that cos(pi/2) = 0 + // https://en.wikipedia.org/wiki/Spherical_trigonometry#Cosine_rules + return Math.acos(Math.cos(a) * Math.cos(b)); + } + + /** + * Compute the area of the spherical right triangle with the given sides. The sides must have lengths + * less than pi/2. + * @param a first side; must be less than pi/2 + * @param b second side; must be less than pi/2 + * @return the area of the spherical right triangle + */ + private static double rightTriangleArea(double a, double b) { + double c = sphericalHypot(a, b); + + // use the spherical law of sines to determine the interior angles + // https://en.wikipedia.org/wiki/Spherical_trigonometry#Sine_rules + double sinC = Math.sin(c); + double angleA = Math.asin(Math.sin(a) / sinC); + double angleB = Math.asin(Math.sin(b) / sinC); + + // use Girard's theorem + return angleA + angleB - (0.5 * PlaneAngleRadians.PI); + } + private static RegionBSPTree2S circleToPolygon(Point2S center, double radius, int numPts, boolean clockwise) { List<Point2S> pts = new ArrayList<>(numPts);
