Skip to content
Permalink
Browse files
GEOMETRY-143: defining strict rules in CutAngle for classification of…
… points near zero; initial fix provided by Ron Goldman
  • Loading branch information
darkma773r committed Jan 29, 2022
1 parent adec3f1 commit e13f15213d9dc46c96759d87660f8dbb05155811
Showing 6 changed files with 206 additions and 35 deletions.
@@ -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 @@ private AngularInterval(final CutAngle minBoundary, final CutAngle maxBoundary)

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 RegionLocation classify(final Point1S pt) {
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 Point1S project(final Point1S pt) {
* @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
@@ -33,24 +33,46 @@
* 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 double offset(final Point1S pt) {
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} */
@@ -210,6 +210,22 @@ public boolean eq(final Point1S other, final Precision.DoubleEquivalence precisi
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
@@ -48,6 +48,35 @@ void testOf_doubles() {
checkFull(AngularInterval.of(0, Angle.TWO_PI, TEST_PRECISION));
}

@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
@@ -802,6 +831,8 @@ private static void checkInterval(final AngularInterval interval, final double m
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());
@@ -167,6 +167,29 @@ void testClassify() {
checkClassify(negPiPos, HyperplaneLocation.PLUS, -0.5, -Angle.PI_OVER_TWO);
}

@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
@@ -266,6 +266,51 @@ void testEq_wrapAround() {
Assertions.assertTrue(c.eq(a, precision));
}

@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

0 comments on commit e13f152

Please sign in to comment.