singhbaljit commented on a change in pull request #97:
URL: https://github.com/apache/commons-geometry/pull/97#discussion_r453927651
##########
File path:
commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
##########
@@ -975,7 +1022,13 @@ private static double rightTriangleArea(final double a,
final double b) {
return angleA + angleB - (0.5 * PlaneAngleRadians.PI);
}
- private static RegionBSPTree2S circleToPolygon(final Point2S center, final
double radius, final int numPts, final boolean clockwise) {
+ private static RegionBSPTree2S circleToPolygon(final Point2S center, final
double radius, final int numPts,
Review comment:
Any reason for this method? I think its easy enough to inline
`circleToPolygon(center, radius, numPts, clockwise, TEST_PRECISION)`.
##########
File path:
commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
##########
@@ -838,6 +843,48 @@ public void testCircleToPolygonBoundarySize() {
Assert.assertEquals("Clockwise boundary size", boundary,
cw.getBoundarySize(), 1.0e-7);
}
+ @Test
+ public void testSmallCircleToPolygon() {
+ // arrange
+ final double radius = 5.0e-8;
+ final Point2S center = Point2S.of(0.5, 1.5);
+ final int numPts = 100;
+
+ // act
+ final RegionBSPTree2S circle = circleToPolygon(center, radius, numPts,
false);
+
+ // assert
+ // https://en.wikipedia.org/wiki/Spherical_cap
+ final double area = 4.0 * PlaneAngleRadians.PI *
Math.pow(Math.sin(radius / 2.0), 2.0);
+ final double boundary = PlaneAngleRadians.TWO_PI * Math.sin(radius);
+
+ SphericalTestUtils.assertPointsEq(center, circle.getCentroid(),
TEST_EPS);
+ Assert.assertEquals(area, circle.getSize(), TEST_EPS);
+ Assert.assertEquals(boundary, circle.getBoundarySize(), TEST_EPS);
+ }
+
+ @Test
+ public void testSmallGeographicalRectangle() {
+ // arrange
+ final double[][] vertices = {
+ {42.656216727628696, -70.61919768884546},
+ {42.65612858998112, -70.61938607250165},
+ {42.65579098923594, -70.61909615581666},
+ {42.655879126692355, -70.61890777301083}
+ };
+
+ // act
+ final RegionBSPTree2S rectangle = latLongToTree(TEST_PRECISION,
vertices);
+
+ // assert
+ // approximate the centroid as average of vertices
+ final double avgLat = Stream.of(vertices).mapToDouble(v ->
v[0]).average().getAsDouble();
+ final double avgLon = Stream.of(vertices).mapToDouble(v ->
v[1]).average().getAsDouble();
+ final Point2S centroid = latLongToPoint(avgLat, avgLon);
+
+ SphericalTestUtils.assertPointsEq(centroid, rectangle.getCentroid(),
TEST_EPS);
Review comment:
We should add tests for size and boundary size.
##########
File path:
commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
##########
@@ -132,35 +139,167 @@ public Point2S getCentroid() {
return weighted == null ? null : Point2S.from(weighted);
}
- /** Returns the weighted vector for the centroid. 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 centroid vectors of multiple
areas, a single
- * centroid can be computed for the whole group.
+ /** Return the weighted centroid vector of the area. The returned vector
points in the direction of the
+ * centroid point on the surface of the unit sphere with the length of the
vector proportional to the
+ * effective mass of the area at the centroid. By adding the weighted
centroid vectors of multiple
+ * convex areas, a single centroid can be computed for the combined area.
* @return weighted centroid 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 getWeightedCentroidVector() {
final List<GreatArc> arcs = getBoundaries();
- switch (arcs.size()) {
+ final int numBoundaries = arcs.size();
+
+ switch (numBoundaries) {
case 0:
// full space; no centroid
return null;
case 1:
- // hemisphere; centroid is the pole of the hemisphere
- final GreatArc singleArc = arcs.get(0);
- return
singleArc.getCircle().getPole().withNorm(singleArc.getSize());
+ // hemisphere
+ return computeHemisphereWeightedCentroidVector(arcs.get(0));
+ case 2:
+ // lune
+ return computeLuneWeightedCentroidVector(arcs.get(0), arcs.get(1));
default:
- // 2 or more sides; use an extension of the approach outlined here:
- // https://archive.org/details/centroidinertiat00broc
- // In short, the centroid is the sum of the pole vectors of each
side
- // multiplied by their arc lengths.
- Vector3D centroid = Vector3D.ZERO;
- for (final GreatArc arc : getBoundaries()) {
- centroid =
centroid.add(arc.getCircle().getPole().withNorm(arc.getSize()));
+ // triangle or other convex polygon
+ final double arcLengthSum = arcs.stream()
+ .mapToDouble(GreatArc::getSize)
+ .sum();
+
+ if (arcLengthSum < TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD) {
+ return computeTriangleFanWeightedCentroidVector(arcs);
+ }
+
+ return computeArcPoleWeightedCentroidVector(arcs);
+ }
+ }
+
+ /** Compute the weighted centroid vector for the hemisphere formed by the
given arc.
+ * @param arc arc defining the hemisphere
+ * @return the weighted centroid vector for the hemisphere
+ * @see #getWeightedCentroidVector()
+ */
+ private Vector3D computeHemisphereWeightedCentroidVector(final GreatArc
arc) {
+ return arc.getCircle().getPole().withNorm(HALF_SIZE);
+ }
+
+ /** Compute the weighted centroid vector for the lune formed by the given
arcs.
+ * @param a first arc for the lune
+ * @param b second arc for the lune
+ * @return the weighted centroid vector for the lune
+ * @see #getWeightedCentroidVector()
+ */
+ private Vector3D computeLuneWeightedCentroidVector(final GreatArc a, final
GreatArc b) {
+ final Point2S aMid = a.getCentroid();
+ final Point2S bMid = b.getCentroid();
+
+ // compute the centroid vector as the exact center of the lune to
avoid inaccurate
+ // results with very small regions
+ final Vector3D.Unit centroid = aMid.slerp(bMid, 0.5).getVector();
+
+ // compute the weight using the reverse of the algorithm from
computeArcPoleWeightedCentroidVector()
+ final double weight =
+ (a.getSize() * centroid.dot(a.getCircle().getPole())) +
+ (b.getSize() * centroid.dot(b.getCircle().getPole()));
+
+ return centroid.withNorm(weight);
+ }
+
+ /** Compute the weighted centroid vector for the triangle or polygon
formed by the given arcs
+ * by adding together the arc pole vectors multiplied by their respective
arc lengths. This
+ * algorithm is described in the paper <a
href="https://archive.org/details/centroidinertiat00broc";>
+ * <em>The Centroid and Inertia Tensor for a Spherical Triangle</em></a>
by John E Brock.
+ *
+ * <p>Note: This algorithm works well in general but is susceptible to
floating point errors
+ * on very small areas. In these cases, the computed centroid may not be
in the expected location
+ * and may even be outside of the area. The {@link
#computeTriangleFanWeightedCentroidVector(List)}
+ * method can produce more accurate results in these cases.</p>
+ * @param arcs boundary arcs for the area
+ * @return the weighted centroid vector for the area
+ * @see #computeTriangleFanWeightedCentroidVector(List)
+ */
+ private Vector3D computeArcPoleWeightedCentroidVector(final List<GreatArc>
arcs) {
+ Vector3D centroid = Vector3D.ZERO;
+
+ Vector3D arcContribution;
Review comment:
Any reason to declare this outside the loop?
##########
File path:
commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
##########
@@ -171,16 +171,32 @@ public RegionBSPTree2S toTree() {
final DoublePrecisionContext precision = ((GreatArc)
getRoot().getCut()).getPrecision();
double sizeSum = 0;
- Vector3D centroidVector = Vector3D.ZERO;
+ Vector3D centroidVectorSum = Vector3D.ZERO;
+
+ Vector3D centroidVector;
+ double maxCentroidVectorWeightSq = 0.0;
for (final ConvexArea2S area : areas) {
sizeSum += area.getSize();
- centroidVector =
centroidVector.add(area.getWeightedCentroidVector());
+
+ centroidVector = area.getWeightedCentroidVector();
+ maxCentroidVectorWeightSq = Math.max(maxCentroidVectorWeightSq,
centroidVector.normSq());
+
+ centroidVectorSum = centroidVectorSum.add(centroidVector);
}
- final Point2S centroid = centroidVector.eq(Vector3D.ZERO, precision) ?
- null :
- Point2S.from(centroidVector);
+ // Convert the weighted centroid vector to a point on the sphere
surface. If the centroid vector
Review comment:
I think you'll disagree with usage of `Optional`, but thought I
recommend it anyways. 😃
```java
final Point2S centroid = Optional.of(centroidVectorSum)
.filter(v -> !(v.normSq() < maxCentroidVectorWeightSq &&
v.eq(Vector3D.ZERO, precision)))
.map(Point2S::from)
.orElse(null);
```
##########
File path:
commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
##########
@@ -171,16 +171,32 @@ public RegionBSPTree2S toTree() {
final DoublePrecisionContext precision = ((GreatArc)
getRoot().getCut()).getPrecision();
double sizeSum = 0;
- Vector3D centroidVector = Vector3D.ZERO;
+ Vector3D centroidVectorSum = Vector3D.ZERO;
+
+ Vector3D centroidVector;
Review comment:
Any reason the declaration is not inside the loop? It doesn't look to be
used elsewhere.
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
[email protected]
