You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2020/01/19 02:53:55 UTC
[commons-geometry] 01/03: GEOMETRY-71: adding additional unit tests
and documentation regarding the computation of spherical barycenters
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 <ma...@hotmail.com>
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);