This is an automated email from the ASF dual-hosted git repository.
mattjuntunen pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-geometry.git
The following commit(s) were added to refs/heads/master by this push:
new e13f152 GEOMETRY-143: defining strict rules in CutAngle for
classification of points near zero; initial fix provided by Ron Goldman
e13f152 is described below
commit e13f15213d9dc46c96759d87660f8dbb05155811
Author: Matt Juntunen <mattjuntu...@apache.org>
AuthorDate: Sun Jan 23 00:15:25 2022 -0500
GEOMETRY-143: defining strict rules in CutAngle for classification of
points near zero; initial fix provided by Ron Goldman
---
.../geometry/spherical/oned/AngularInterval.java | 38 ++++++----
.../commons/geometry/spherical/oned/CutAngle.java | 88 +++++++++++++++++-----
.../commons/geometry/spherical/oned/Point1S.java | 16 ++++
.../spherical/oned/AngularIntervalTest.java | 31 ++++++++
.../geometry/spherical/oned/CutAngleTest.java | 23 ++++++
.../geometry/spherical/oned/Point1STest.java | 45 +++++++++++
6 files changed, 206 insertions(+), 35 deletions(-)
diff --git
a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/AngularInterval.java
b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/AngularInterval.java
index 906ac58..ef7e457 100644
---
a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/AngularInterval.java
+++
b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/AngularInterval.java
@@ -47,6 +47,9 @@ public class AngularInterval implements
HyperplaneBoundedRegion<Point1S> {
/** Point halfway between the min and max boundaries. */
private final Point1S midpoint;
+ /** Flag set to true if the interval wraps around the {@code 0/2pi} point.
*/
+ private final boolean wraps;
+
/** Construct a new instance representing the angular region between the
given
* min and max azimuth boundaries. The arguments must be either all finite
or all
* null (to indicate the full space). If the boundaries are finite, then
the min
@@ -59,9 +62,22 @@ public class AngularInterval implements
HyperplaneBoundedRegion<Point1S> {
this.minBoundary = minBoundary;
this.maxBoundary = maxBoundary;
- this.midpoint = (minBoundary != null && maxBoundary != null) ?
- Point1S.of(0.5 * (minBoundary.getAzimuth() +
maxBoundary.getAzimuth())) :
- null;
+
+ Point1S midpointVal = null;
+ boolean wrapsVal = false;
+
+ if (minBoundary != null && maxBoundary != null) {
+ midpointVal = Point1S.of(0.5 * (minBoundary.getAzimuth() +
maxBoundary.getAzimuth()));
+
+ // The interval wraps zero if the max boundary lies on the other
side of zero than
+ // the min. This is a more reliable way to compute the wrapping
flag than direct
+ // comparison of the normalized azimuths since this approach takes
into account
+ // azimuths that are equivalent to zero.
+ wrapsVal = minBoundary.classify(maxBoundary.getPoint()) ==
HyperplaneLocation.PLUS;
+ }
+
+ this.midpoint = midpointVal;
+ this.wraps = wrapsVal;
}
/** Get the minimum azimuth angle for the interval, or {@code 0}
@@ -163,13 +179,11 @@ public class AngularInterval implements
HyperplaneBoundedRegion<Point1S> {
final HyperplaneLocation minLoc = minBoundary.classify(pt);
final HyperplaneLocation maxLoc = maxBoundary.classify(pt);
- final boolean wraps = wrapsZero();
-
- if ((!wraps && (minLoc == HyperplaneLocation.PLUS || maxLoc ==
HyperplaneLocation.PLUS)) ||
+ if (minLoc == HyperplaneLocation.ON || maxLoc ==
HyperplaneLocation.ON) {
+ return RegionLocation.BOUNDARY;
+ } else if ((!wraps && (minLoc == HyperplaneLocation.PLUS || maxLoc
== HyperplaneLocation.PLUS)) ||
(wraps && minLoc == HyperplaneLocation.PLUS && maxLoc ==
HyperplaneLocation.PLUS)) {
return RegionLocation.OUTSIDE;
- } else if (minLoc == HyperplaneLocation.ON || maxLoc ==
HyperplaneLocation.ON) {
- return RegionLocation.BOUNDARY;
}
}
return RegionLocation.INSIDE;
@@ -195,13 +209,7 @@ public class AngularInterval implements
HyperplaneBoundedRegion<Point1S> {
* @return true if the interval wraps around the zero/{@code 2pi} point
*/
public boolean wrapsZero() {
- if (!isFull()) {
- final double minNormAz =
minBoundary.getPoint().getNormalizedAzimuth();
- final double maxNormAz =
maxBoundary.getPoint().getNormalizedAzimuth();
-
- return maxNormAz < minNormAz;
- }
- return false;
+ return wraps;
}
/** Return a new instance transformed by the argument. If the transformed
size
diff --git
a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/CutAngle.java
b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/CutAngle.java
index ad760a4..dcb993f 100644
---
a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/CutAngle.java
+++
b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/CutAngle.java
@@ -33,24 +33,46 @@ import org.apache.commons.numbers.core.Precision;
* meaning an azimuth angle and a direction (increasing or decreasing angles)
* along the circle.
*
+ * <p><strong>Point Classification</strong></p>
* <p>Hyperplanes split the spaces they are embedded in into three distinct
parts:
* the hyperplane itself, a plus side and a minus side. However, since
spherical
* space wraps around, a single oriented point is not sufficient to partition
the space;
* any point could be classified as being on the plus or minus side of a
hyperplane
* depending on the direction that the circle is traversed. The approach taken
in this
- * class to address this issue is to (1) define a second, implicit cut point
at {@code 0pi} and
+ * class to address this issue is to (1) define a second, implicit cut point
at \(0\pi\) and
* (2) define the domain of hyperplane points (for partitioning purposes) to
be the
- * range {@code [0, 2pi)}. Each hyperplane then splits the space into the
intervals
- * {@code [0, x]} and {@code [x, 2pi)}, where {@code x} is the location of the
hyperplane.
+ * range \([0, 2\pi)\). Each hyperplane then splits the space into the
intervals
+ * \([0, x]\) and \([x, 2\pi)\), where \(x\) is the location of the hyperplane.
* One way to visualize this is to picture the circle as a cake that has
already been
- * cut at {@code 0pi}. Each hyperplane then specifies the location of the
second
+ * cut at \(0\pi\). Each hyperplane then specifies the location of the second
* cut of the cake, with the plus and minus sides being the pieces thus cut.
* </p>
*
- * <p>Note that with the hyperplane partitioning rules described above, the
hyperplane
- * at {@code 0pi} is unique in that it has the entire space on one side (minus
the hyperplane
- * itself) and no points whatsoever on the other. This is very different from
hyperplanes in
- * Euclidean space, which always have infinitely many points on both sides.</p>
+ * <p>Note that with the hyperplane partitioning approach described above, the
hyperplane
+ * at \(0\pi\) is unique in that it has the entire space on one side (except
for the
+ * hyperplane itself) and no points whatsoever on the other. This is very
different from
+ * hyperplanes in Euclidean space, which always have infinitely many points on
both sides.</p>
+ *
+ * <p>Due to the unique status of the \(0\pi\) point, special care must be
given to azimuths very
+ * close to this value. The rules below define exactly how points are
classified when the hyperplane
+ * point, the test point, or both lie very close to \(0\pi\). In what follows,
\(H\) represents the
+ * hyperplane point, \(P\) the point being classified, and the symbol
\(\approx\) means equivalent as
+ * evaluated by the instance's {@link #getPrecision() precision context}. Note
that points are considered
+ * equivalent to \(0\pi\) if their normalized azimuths are close to either
\(0\pi\) or \(2\pi\).
+ * <ul>
+ * <li>\(H \approx P\) \(\implies\) \(P\) is classified as {@link
HyperplaneLocation#ON ON}.</li>
+ * <li>\(H \approx 0\) and \(P \approx 0\) \(\implies\) \(P\) is classified as
+ * {@link HyperplaneLocation#ON ON}.</li>
+ * <li>\(H \approx 0\) and \(P \neq 0\) \(\implies\) \(P\) is classified as
{@link HyperplaneLocation#PLUS PLUS}
+ * if the cut is positive facing and {@link HyperplaneLocation#MINUS
MINUS} if negative facing.</li>
+ * <li>\(H \neq 0\) and \(P \approx 0\) \(\implies\) \(P\) is classified as
{@link HyperplaneLocation#MINUS MINUS}
+ * if the cut is positive facing and {@link HyperplaneLocation#PLUS PLUS}
if negative facing.</li>
+ * <li>\(H \neq 0\) and \(P \neq 0\) \(\implies\) The normalized azimuths of
\(H\) and \(P\) are compared and the
+ * standard rules applied. If \(P \gt H\), then \(P\) is classified as
{@link HyperplaneLocation#PLUS PLUS}
+ * if the cut is positive facing and {@link HyperplaneLocation#MINUS
MINUS} if negative facing. If
+ * \(P \lt H\), then \(P\) is classified as {@link
HyperplaneLocation#MINUS MINUS} if the cut is
+ * positive facing and {@link HyperplaneLocation#PLUS PLUS} if negative
facing.</li>
+ * </ul>
*
* <p>Instances of this class are guaranteed to be immutable.</p>
* @see CutAngles
@@ -140,25 +162,51 @@ public final class CutAngle extends
AbstractHyperplane<Point1S> {
return positiveFacing ? +dist : -dist;
}
- /** {@inheritDoc} */
+ /** {@inheritDoc}
+ *
+ * <p>Special classification rules are applied in this method due to the
unique nature
+ * of spherical space. See the class level documentation for details.</p>
+ */
@Override
public HyperplaneLocation classify(final Point1S pt) {
+ int cmp = classifyPositiveFacing(pt);
+ if (!positiveFacing) {
+ cmp = -cmp;
+ }
+
+ if (cmp < 0) {
+ return HyperplaneLocation.MINUS;
+ } else if (cmp > 0) {
+ return HyperplaneLocation.PLUS;
+ }
+ return HyperplaneLocation.ON;
+ }
+
+ /** Classify the given point assuming a positive facing cut.
+ * @param pt point to classify
+ * @return int value indicating the classification region of {@code pt}
+ */
+ private int classifyPositiveFacing(final Point1S pt) {
final Precision.DoubleEquivalence precision = getPrecision();
- final Point1S compPt = Point1S.ZERO.eq(pt, precision) ?
- Point1S.ZERO :
- pt;
+ final double az = pt.getNormalizedAzimuth();
+ final double base = this.point.getNormalizedAzimuth();
- final double offsetValue = offset(compPt);
- final double cmp = precision.signum(offsetValue);
+ final int cmp = precision.compare(az, base);
- if (cmp > 0) {
- return HyperplaneLocation.PLUS;
- } else if (cmp < 0) {
- return HyperplaneLocation.MINUS;
- }
+ if (cmp != 0) {
+ final boolean azIsZero = pt.eqZero(precision);
+ final boolean baseIsZero = this.point.eqZero(precision);
- return HyperplaneLocation.ON;
+ if (baseIsZero) {
+ return azIsZero ?
+ 0 :
+ +1;
+ } else if (azIsZero) {
+ return -1;
+ }
+ }
+ return cmp;
}
/** {@inheritDoc} */
diff --git
a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/Point1S.java
b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/Point1S.java
index b5a1935..916d21b 100644
---
a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/Point1S.java
+++
b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/oned/Point1S.java
@@ -210,6 +210,22 @@ public final class Point1S implements Point<Point1S> {
return precision.eqZero(dist);
}
+ /** Return true if this instance is equivalent to the {@link Point1S#ZERO
zero point},
+ * meaning that the shortest angular distance of this point from the zero
point is equal
+ * to zero as evaluated by the given precision context. This method is
functionally equivalent
+ * to {@code pt.eq(Point1S.ZERO, precision)} but with simplified logic due
to the comparison
+ * point being zero.
+ * @param precision precision context used for floating point comparisons
+ * @return true if this instance is equivalent to the zero point
+ * @see #eq(Point1S,
org.apache.commons.numbers.core.Precision.DoubleEquivalence)
+ */
+ public boolean eqZero(final Precision.DoubleEquivalence precision) {
+ final double closest = normalizedAzimuth > Math.PI ?
+ Angle.TWO_PI :
+ 0d;
+ return precision.eq(normalizedAzimuth, closest);
+ }
+
/**
* Get a hashCode for the point. Points normally must have exactly the
* same azimuth angles in order to have the same hash code. Points
diff --git
a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/AngularIntervalTest.java
b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/AngularIntervalTest.java
index 28361af..f245abd 100644
---
a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/AngularIntervalTest.java
+++
b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/AngularIntervalTest.java
@@ -49,6 +49,35 @@ class AngularIntervalTest {
}
@Test
+ void testOf_endPointsCloseToZero() {
+ // arrange
+ final double pi = Math.PI;
+
+ final double belowZero = -5e-11;
+ final double aboveZero = 5e-11;
+
+ final double belowTwoPi = Angle.TWO_PI - 5e-11;
+ final double aboveTwoPi = Angle.TWO_PI + 5e-11;
+
+ // act/assert
+ checkInterval(AngularInterval.of(belowZero, pi, TEST_PRECISION),
belowZero, pi);
+ checkInterval(AngularInterval.of(aboveZero, pi, TEST_PRECISION),
aboveZero, pi);
+
+ checkInterval(AngularInterval.of(belowTwoPi, pi, TEST_PRECISION),
belowTwoPi, pi + Angle.TWO_PI);
+ checkInterval(AngularInterval.of(aboveTwoPi, pi, TEST_PRECISION),
aboveTwoPi, pi + Angle.TWO_PI);
+
+ checkInterval(AngularInterval.of(pi, belowZero, TEST_PRECISION), pi,
belowZero + Angle.TWO_PI);
+ checkInterval(AngularInterval.of(pi, aboveZero, TEST_PRECISION), pi,
aboveZero + Angle.TWO_PI);
+
+ checkInterval(AngularInterval.of(pi, belowTwoPi, TEST_PRECISION), pi,
belowTwoPi);
+ checkInterval(AngularInterval.of(pi, aboveTwoPi, TEST_PRECISION), pi,
aboveTwoPi);
+
+ // from GEOMETRY-143
+ checkInterval(AngularInterval.of(6,
Double.parseDouble("0x1.921fb54442c8ep2"), TEST_PRECISION),
+ 6, Double.parseDouble("0x1.921fb54442c8ep2"));
+ }
+
+ @Test
void testOf_doubles_invalidArgs() {
// act/assert
Assertions.assertThrows(IllegalArgumentException.class, () ->
AngularInterval.of(Double.NEGATIVE_INFINITY, 0, TEST_PRECISION));
@@ -802,6 +831,8 @@ class AngularIntervalTest {
Assertions.assertEquals(0, interval.getBoundarySize(), TEST_EPS);
Assertions.assertEquals(max - min, interval.getSize(), TEST_EPS);
+ checkClassify(interval, RegionLocation.BOUNDARY,
interval.getMinBoundary().getPoint());
+
checkClassify(interval, RegionLocation.INSIDE, interval.getMidPoint());
checkClassify(interval, RegionLocation.BOUNDARY,
interval.getMinBoundary().getPoint(),
interval.getMaxBoundary().getPoint());
diff --git
a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/CutAngleTest.java
b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/CutAngleTest.java
index 95a0ad2..48fc2fb 100644
---
a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/CutAngleTest.java
+++
b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/CutAngleTest.java
@@ -168,6 +168,29 @@ class CutAngleTest {
}
@Test
+ void testClassify_azimuthsCloseToZero() {
+ // arrange
+ final CutAngle belowZeroPos = CutAngles.createPositiveFacing(-5e-11,
TEST_PRECISION);
+ final CutAngle belowZeroNeg = CutAngles.createNegativeFacing(-5e-11,
TEST_PRECISION);
+
+ final CutAngle aboveZeroPos = CutAngles.createPositiveFacing(5e-11,
TEST_PRECISION);
+ final CutAngle aboveZeroNeg = CutAngles.createNegativeFacing(5e-11,
TEST_PRECISION);
+
+ // act/assert
+ checkClassify(belowZeroPos, HyperplaneLocation.PLUS, -1.6e-10,
1.2e-10, 1.6e-10);
+ checkClassify(belowZeroPos, HyperplaneLocation.ON, -1.2e-10, -8e-11,
-4e-11, 0, 4e-11, 8e-11);
+
+ checkClassify(belowZeroNeg, HyperplaneLocation.MINUS, -1.6e-10,
1.2e-10, 1.6e-10);
+ checkClassify(belowZeroNeg, HyperplaneLocation.ON, -1.2e-10, -8e-11,
-4e-11, 0, 4e-11, 8e-11);
+
+ checkClassify(aboveZeroPos, HyperplaneLocation.PLUS, -1.6e-10,
-1.2e-10, 1.6e-10);
+ checkClassify(aboveZeroPos, HyperplaneLocation.ON, -8e-11, -4e-11, 0,
4e-11, 8e-11, 1.2e-10);
+
+ checkClassify(aboveZeroNeg, HyperplaneLocation.MINUS, -1.6e-10,
-1.2e-10, 1.6e-10);
+ checkClassify(aboveZeroNeg, HyperplaneLocation.ON, -8e-11, -4e-11, 0,
4e-11, 8e-11, 1.2e-10);
+ }
+
+ @Test
void testContains() {
// arrange
final CutAngle pt = CutAngles.createNegativeFacing(Angle.PI_OVER_TWO,
TEST_PRECISION);
diff --git
a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/Point1STest.java
b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/Point1STest.java
index 54609b7..7105019 100644
---
a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/Point1STest.java
+++
b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/oned/Point1STest.java
@@ -267,6 +267,51 @@ class Point1STest {
}
@Test
+ void testEqZero() {
+ // arrange
+ final double eps = 1e-8;
+ final double delta = 1e-3 * eps;
+ final Precision.DoubleEquivalence precision =
Precision.doubleEquivalenceOfEpsilon(eps);
+
+ // act/assert
+ checkEqZero(0, precision);
+ checkEqZero(Angle.TWO_PI, precision);
+ checkEqZero(-Angle.TWO_PI, precision);
+
+ checkNotEqZero(Math.PI, precision);
+ checkNotEqZero(-Math.PI, precision);
+
+ checkNotEqZero(Angle.PI_OVER_TWO, precision);
+ checkNotEqZero(-Angle.PI_OVER_TWO, precision);
+
+ checkNotEqZero(-eps - delta, precision);
+ checkNotEqZero(-eps - delta - Angle.TWO_PI, precision);
+
+ for (double a = -eps + delta; a < eps; a += delta) {
+ checkEqZero(a, precision);
+ checkEqZero(a - Angle.TWO_PI, precision);
+ checkEqZero(a + Angle.TWO_PI, precision);
+ }
+
+ checkNotEqZero(eps + delta, precision);
+ checkNotEqZero(eps + delta + Angle.TWO_PI, precision);
+ }
+
+ private void checkEqZero(final double az, final
Precision.DoubleEquivalence precision) {
+ final Point1S pt = Point1S.of(az);
+
+ Assertions.assertTrue(pt.eqZero(precision));
+ Assertions.assertTrue(Point1S.ZERO.eq(pt, precision));
+ }
+
+ private void checkNotEqZero(final double az, final
Precision.DoubleEquivalence precision) {
+ final Point1S pt = Point1S.of(az);
+
+ Assertions.assertFalse(pt.eqZero(precision));
+ Assertions.assertFalse(Point1S.ZERO.eq(pt, precision));
+ }
+
+ @Test
void testDistance() {
// arrange
final Point1S a = Point1S.of(0.0);
