com.google.common.geometry.S2 Maven / Gradle / Ivy
/*
* Copyright 2005 Google Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package com.google.common.geometry;
import com.google.common.annotations.VisibleForTesting;
import com.google.common.base.Preconditions;
public final strictfp class S2 {
// Declare some frequently used constants
public static final double M_PI = Math.PI;
public static final double M_1_PI = 1.0 / Math.PI;
public static final double M_PI_2 = Math.PI / 2.0;
public static final double M_PI_4 = Math.PI / 4.0;
public static final double M_SQRT2 = Math.sqrt(2);
public static final double M_E = Math.E;
// Together these flags define a cell orientation. If SWAP_MASK
// is true, then canonical traversal order is flipped around the
// diagonal (i.e. i and j are swapped with each other). If
// INVERT_MASK is true, then the traversal order is rotated by 180
// degrees (i.e. the bits of i and j are inverted, or equivalently,
// the axis directions are reversed).
public static final int SWAP_MASK = 0x01;
public static final int INVERT_MASK = 0x02;
// Number of bits in the mantissa of a double.
private static final int EXPONENT_SHIFT = 52;
// Mask to extract the exponent from a double.
private static final long EXPONENT_MASK = 0x7ff0000000000000L;
/**
* If v is non-zero, return an integer {@code exp} such that
* {@code (0.5 <= |v|*2^(-exp) < 1)}. If v is zero, return 0.
*
* Note that this arguably a bad definition of exponent because it makes
* {@code exp(9) == 4}. In decimal this would be like saying that the
* exponent of 1234 is 4, when in scientific 'exponent' notation 1234 is
* {@code 1.234 x 10^3}.
*
* TODO(dbeaumont): Replace this with "DoubleUtils.getExponent(v) - 1" ?
*/
@VisibleForTesting
static int exp(double v) {
if (v == 0) {
return 0;
}
long bits = Double.doubleToLongBits(v);
return (int) ((EXPONENT_MASK & bits) >> EXPONENT_SHIFT) - 1022;
}
/** Mapping Hilbert traversal order to orientation adjustment mask. */
private static final int[] POS_TO_ORIENTATION =
{SWAP_MASK, 0, 0, INVERT_MASK + SWAP_MASK};
/**
* Returns an XOR bit mask indicating how the orientation of a child subcell
* is related to the orientation of its parent cell. The returned value can
* be XOR'd with the parent cell's orientation to give the orientation of
* the child cell.
*
* @param position the position of the subcell in the Hilbert traversal, in
* the range [0,3].
* @return a bit mask containing some combination of {@link #SWAP_MASK} and
* {@link #INVERT_MASK}.
* @throws IllegalArgumentException if position is out of bounds.
*/
public static int posToOrientation(int position) {
Preconditions.checkArgument(0 <= position && position < 4);
return POS_TO_ORIENTATION[position];
}
/** Mapping from cell orientation + Hilbert traversal to IJ-index. */
private static final int[][] POS_TO_IJ = {
// 0 1 2 3
{0, 1, 3, 2}, // canonical order: (0,0), (0,1), (1,1), (1,0)
{0, 2, 3, 1}, // axes swapped: (0,0), (1,0), (1,1), (0,1)
{3, 2, 0, 1}, // bits inverted: (1,1), (1,0), (0,0), (0,1)
{3, 1, 0, 2}, // swapped & inverted: (1,1), (0,1), (0,0), (1,0)
};
/**
* Return the IJ-index of the subcell at the given position in the Hilbert
* curve traversal with the given orientation. This is the inverse of
* {@link #ijToPos}.
*
* @param orientation the subcell orientation, in the range [0,3].
* @param position the position of the subcell in the Hilbert traversal, in
* the range [0,3].
* @return the IJ-index where {@code 0->(0,0), 1->(0,1), 2->(1,0), 3->(1,1)}.
* @throws IllegalArgumentException if either parameter is out of bounds.
*/
public static int posToIJ(int orientation, int position) {
Preconditions.checkArgument(0 <= orientation && orientation < 4);
Preconditions.checkArgument(0 <= position && position < 4);
return POS_TO_IJ[orientation][position];
}
/** Mapping from Hilbert traversal order + cell orientation to IJ-index. */
private static final int IJ_TO_POS[][] = {
// (0,0) (0,1) (1,0) (1,1)
{0, 1, 3, 2}, // canonical order
{0, 3, 1, 2}, // axes swapped
{2, 3, 1, 0}, // bits inverted
{2, 1, 3, 0}, // swapped & inverted
};
/**
* Returns the order in which a specified subcell is visited by the Hilbert
* curve. This is the inverse of {@link #posToIJ}.
*
* @param orientation the subcell orientation, in the range [0,3].
* @param ijIndex the subcell index where
* {@code 0->(0,0), 1->(0,1), 2->(1,0), 3->(1,1)}.
* @return the position of the subcell in the Hilbert traversal, in the range
* [0,3].
* @throws IllegalArgumentException if either parameter is out of bounds.
*/
public static final int ijToPos(int orientation, int ijIndex) {
Preconditions.checkArgument(0 <= orientation && orientation < 4);
Preconditions.checkArgument(0 <= ijIndex && ijIndex < 4);
return IJ_TO_POS[orientation][ijIndex];
}
/**
* Defines an area or a length cell metric.
*/
public static class Metric {
private final double deriv;
private final int dim;
/**
* Defines a cell metric of the given dimension (1 == length, 2 == area).
*/
public Metric(int dim, double deriv) {
this.deriv = deriv;
this.dim = dim;
}
/**
* The "deriv" value of a metric is a derivative, and must be multiplied by
* a length or area in (s,t)-space to get a useful value.
*/
public double deriv() {
return deriv;
}
/** Return the value of a metric for cells at the given level. */
public double getValue(int level) {
return StrictMath.scalb(deriv, dim * (1 - level));
}
/**
* Return the level at which the metric has approximately the given value.
* For example, S2::kAvgEdge.GetClosestLevel(0.1) returns the level at which
* the average cell edge length is approximately 0.1. The return value is
* always a valid level.
*/
public int getClosestLevel(double value) {
return getMinLevel(M_SQRT2 * value);
}
/**
* Return the minimum level such that the metric is at most the given value,
* or S2CellId::kMaxLevel if there is no such level. For example,
* S2::kMaxDiag.GetMinLevel(0.1) returns the minimum level such that all
* cell diagonal lengths are 0.1 or smaller. The return value is always a
* valid level.
*/
public int getMinLevel(double value) {
if (value <= 0) {
return S2CellId.MAX_LEVEL;
}
// This code is equivalent to computing a floating-point "level"
// value and rounding up.
int exponent = exp(value / ((1 << dim) * deriv));
int level = Math.max(0,
Math.min(S2CellId.MAX_LEVEL, -((exponent - 1) >> (dim - 1))));
// assert (level == S2CellId.MAX_LEVEL || getValue(level) <= value);
// assert (level == 0 || getValue(level - 1) > value);
return level;
}
/**
* Return the maximum level such that the metric is at least the given
* value, or zero if there is no such level. For example,
* S2.kMinWidth.GetMaxLevel(0.1) returns the maximum level such that all
* cells have a minimum width of 0.1 or larger. The return value is always a
* valid level.
*/
public int getMaxLevel(double value) {
if (value <= 0) {
return S2CellId.MAX_LEVEL;
}
// This code is equivalent to computing a floating-point "level"
// value and rounding down.
int exponent = exp((1 << dim) * deriv / value);
int level = Math.max(0,
Math.min(S2CellId.MAX_LEVEL, ((exponent - 1) >> (dim - 1))));
// assert (level == 0 || getValue(level) >= value);
// assert (level == S2CellId.MAX_LEVEL || getValue(level + 1) < value);
return level;
}
}
/**
* Return a unique "origin" on the sphere for operations that need a fixed
* reference point. It should *not* be a point that is commonly used in edge
* tests in order to avoid triggering code to handle degenerate cases. (This
* rules out the north and south poles.)
*/
public static S2Point origin() {
return new S2Point(0, 1, 0);
}
/**
* Return true if the given point is approximately unit length (this is mainly
* useful for assertions).
*/
public static boolean isUnitLength(S2Point p) {
return Math.abs(p.norm2() - 1) <= 1e-15;
}
/**
* Return true if edge AB crosses CD at a point that is interior to both
* edges. Properties:
*
* (1) SimpleCrossing(b,a,c,d) == SimpleCrossing(a,b,c,d) (2)
* SimpleCrossing(c,d,a,b) == SimpleCrossing(a,b,c,d)
*/
public static boolean simpleCrossing(S2Point a, S2Point b, S2Point c, S2Point d) {
// We compute SimpleCCW() for triangles ACB, CBD, BDA, and DAC. All
// of these triangles need to have the same orientation (CW or CCW)
// for an intersection to exist. Note that this is slightly more
// restrictive than the corresponding definition for planar edges,
// since we need to exclude pairs of line segments that would
// otherwise "intersect" by crossing two antipodal points.
S2Point ab = S2Point.crossProd(a, b);
S2Point cd = S2Point.crossProd(c, d);
double acb = -ab.dotProd(c);
double cbd = -cd.dotProd(b);
double bda = ab.dotProd(d);
double dac = cd.dotProd(a);
return (acb * cbd > 0) && (cbd * bda > 0) && (bda * dac > 0);
}
/**
* Return a vector "c" that is orthogonal to the given unit-length vectors "a"
* and "b". This function is similar to a.CrossProd(b) except that it does a
* better job of ensuring orthogonality when "a" is nearly parallel to "b",
* and it returns a non-zero result even when a == b or a == -b.
*
* It satisfies the following properties (RCP == RobustCrossProd):
*
* (1) RCP(a,b) != 0 for all a, b (2) RCP(b,a) == -RCP(a,b) unless a == b or
* a == -b (3) RCP(-a,b) == -RCP(a,b) unless a == b or a == -b (4) RCP(a,-b)
* == -RCP(a,b) unless a == b or a == -b
*/
public static S2Point robustCrossProd(S2Point a, S2Point b) {
// The direction of a.CrossProd(b) becomes unstable as (a + b) or (a - b)
// approaches zero. This leads to situations where a.CrossProd(b) is not
// very orthogonal to "a" and/or "b". We could fix this using Gram-Schmidt,
// but we also want b.RobustCrossProd(a) == -b.RobustCrossProd(a).
//
// The easiest fix is to just compute the cross product of (b+a) and (b-a).
// Given that "a" and "b" are unit-length, this has good orthogonality to
// "a" and "b" even if they differ only in the lowest bit of one component.
// assert (isUnitLength(a) && isUnitLength(b));
S2Point x = S2Point.crossProd(S2Point.add(b, a), S2Point.sub(b, a));
if (!x.equals(new S2Point(0, 0, 0))) {
return x;
}
// The only result that makes sense mathematically is to return zero, but
// we find it more convenient to return an arbitrary orthogonal vector.
return ortho(a);
}
/**
* Return a unit-length vector that is orthogonal to "a". Satisfies Ortho(-a)
* = -Ortho(a) for all a.
*/
public static S2Point ortho(S2Point a) {
// The current implementation in S2Point has the property we need,
// i.e. Ortho(-a) = -Ortho(a) for all a.
return a.ortho();
}
/**
* Return the area of triangle ABC. The method used is about twice as
* expensive as Girard's formula, but it is numerically stable for both large
* and very small triangles. The points do not need to be normalized. The area
* is always positive.
*
* The triangle area is undefined if it contains two antipodal points, and
* becomes numerically unstable as the length of any edge approaches 180
* degrees.
*/
static double area(S2Point a, S2Point b, S2Point c) {
// This method is based on l'Huilier's theorem,
//
// tan(E/4) = sqrt(tan(s/2) tan((s-a)/2) tan((s-b)/2) tan((s-c)/2))
//
// where E is the spherical excess of the triangle (i.e. its area),
// a, b, c, are the side lengths, and
// s is the semiperimeter (a + b + c) / 2 .
//
// The only significant source of error using l'Huilier's method is the
// cancellation error of the terms (s-a), (s-b), (s-c). This leads to a
// *relative* error of about 1e-16 * s / min(s-a, s-b, s-c). This compares
// to a relative error of about 1e-15 / E using Girard's formula, where E is
// the true area of the triangle. Girard's formula can be even worse than
// this for very small triangles, e.g. a triangle with a true area of 1e-30
// might evaluate to 1e-5.
//
// So, we prefer l'Huilier's formula unless dmin < s * (0.1 * E), where
// dmin = min(s-a, s-b, s-c). This basically includes all triangles
// except for extremely long and skinny ones.
//
// Since we don't know E, we would like a conservative upper bound on
// the triangle area in terms of s and dmin. It's possible to show that
// E <= k1 * s * sqrt(s * dmin), where k1 = 2*sqrt(3)/Pi (about 1).
// Using this, it's easy to show that we should always use l'Huilier's
// method if dmin >= k2 * s^5, where k2 is about 1e-2. Furthermore,
// if dmin < k2 * s^5, the triangle area is at most k3 * s^4, where
// k3 is about 0.1. Since the best case error using Girard's formula
// is about 1e-15, this means that we shouldn't even consider it unless
// s >= 3e-4 or so.
// We use volatile doubles to force the compiler to truncate all of these
// quantities to 64 bits. Otherwise it may compute a value of dmin > 0
// simply because it chose to spill one of the intermediate values to
// memory but not one of the others.
final double sa = b.angle(c);
final double sb = c.angle(a);
final double sc = a.angle(b);
final double s = 0.5 * (sa + sb + sc);
if (s >= 3e-4) {
// Consider whether Girard's formula might be more accurate.
double s2 = s * s;
double dmin = s - Math.max(sa, Math.max(sb, sc));
if (dmin < 1e-2 * s * s2 * s2) {
// This triangle is skinny enough to consider Girard's formula.
double area = girardArea(a, b, c);
if (dmin < s * (0.1 * area)) {
return area;
}
}
}
// Use l'Huilier's formula.
return 4
* Math.atan(
Math.sqrt(
Math.max(0.0,
Math.tan(0.5 * s) * Math.tan(0.5 * (s - sa)) * Math.tan(0.5 * (s - sb))
* Math.tan(0.5 * (s - sc)))));
}
/**
* Return the area of the triangle computed using Girard's formula. This is
* slightly faster than the Area() method above is not accurate for very small
* triangles.
*/
public static double girardArea(S2Point a, S2Point b, S2Point c) {
// This is equivalent to the usual Girard's formula but is slightly
// more accurate, faster to compute, and handles a == b == c without
// a special case.
S2Point ab = S2Point.crossProd(a, b);
S2Point bc = S2Point.crossProd(b, c);
S2Point ac = S2Point.crossProd(a, c);
return Math.max(0.0, ab.angle(ac) - ab.angle(bc) + bc.angle(ac));
}
/**
* Like Area(), but returns a positive value for counterclockwise triangles
* and a negative value otherwise.
*/
public static double signedArea(S2Point a, S2Point b, S2Point c) {
return area(a, b, c) * robustCCW(a, b, c);
}
// About centroids:
// ----------------
//
// There are several notions of the "centroid" of a triangle. First, there
// // is the planar centroid, which is simply the centroid of the ordinary
// (non-spherical) triangle defined by the three vertices. Second, there is
// the surface centroid, which is defined as the intersection of the three
// medians of the spherical triangle. It is possible to show that this
// point is simply the planar centroid projected to the surface of the
// sphere. Finally, there is the true centroid (mass centroid), which is
// defined as the area integral over the spherical triangle of (x,y,z)
// divided by the triangle area. This is the point that the triangle would
// rotate around if it was spinning in empty space.
//
// The best centroid for most purposes is the true centroid. Unlike the
// planar and surface centroids, the true centroid behaves linearly as
// regions are added or subtracted. That is, if you split a triangle into
// pieces and compute the average of their centroids (weighted by triangle
// area), the result equals the centroid of the original triangle. This is
// not true of the other centroids.
//
// Also note that the surface centroid may be nowhere near the intuitive
// "center" of a spherical triangle. For example, consider the triangle
// with vertices A=(1,eps,0), B=(0,0,1), C=(-1,eps,0) (a quarter-sphere).
// The surface centroid of this triangle is at S=(0, 2*eps, 1), which is
// within a distance of 2*eps of the vertex B. Note that the median from A
// (the segment connecting A to the midpoint of BC) passes through S, since
// this is the shortest path connecting the two endpoints. On the other
// hand, the true centroid is at M=(0, 0.5, 0.5), which when projected onto
// the surface is a much more reasonable interpretation of the "center" of
// this triangle.
/**
* Return the centroid of the planar triangle ABC. This can be normalized to
* unit length to obtain the "surface centroid" of the corresponding spherical
* triangle, i.e. the intersection of the three medians. However, note that
* for large spherical triangles the surface centroid may be nowhere near the
* intuitive "center" (see example above).
*/
public static S2Point planarCentroid(S2Point a, S2Point b, S2Point c) {
return new S2Point((a.x + b.x + c.x) / 3.0, (a.y + b.y + c.y) / 3.0, (a.z + b.z + c.z) / 3.0);
}
/**
* Returns the true centroid of the spherical triangle ABC multiplied by the
* signed area of spherical triangle ABC. The reasons for multiplying by the
* signed area are (1) this is the quantity that needs to be summed to compute
* the centroid of a union or difference of triangles, and (2) it's actually
* easier to calculate this way.
*/
public static S2Point trueCentroid(S2Point a, S2Point b, S2Point c) {
// I couldn't find any references for computing the true centroid of a
// spherical triangle... I have a truly marvellous demonstration of this
// formula which this margin is too narrow to contain :)
// assert (isUnitLength(a) && isUnitLength(b) && isUnitLength(c));
double sina = S2Point.crossProd(b, c).norm();
double sinb = S2Point.crossProd(c, a).norm();
double sinc = S2Point.crossProd(a, b).norm();
double ra = (sina == 0) ? 1 : (Math.asin(sina) / sina);
double rb = (sinb == 0) ? 1 : (Math.asin(sinb) / sinb);
double rc = (sinc == 0) ? 1 : (Math.asin(sinc) / sinc);
// Now compute a point M such that M.X = rX * det(ABC) / 2 for X in A,B,C.
S2Point x = new S2Point(a.x, b.x, c.x);
S2Point y = new S2Point(a.y, b.y, c.y);
S2Point z = new S2Point(a.z, b.z, c.z);
S2Point r = new S2Point(ra, rb, rc);
return new S2Point(0.5 * S2Point.crossProd(y, z).dotProd(r),
0.5 * S2Point.crossProd(z, x).dotProd(r), 0.5 * S2Point.crossProd(x, y).dotProd(r));
}
/**
* Return true if the points A, B, C are strictly counterclockwise. Return
* false if the points are clockwise or colinear (i.e. if they are all
* contained on some great circle).
*
* Due to numerical errors, situations may arise that are mathematically
* impossible, e.g. ABC may be considered strictly CCW while BCA is not.
* However, the implementation guarantees the following:
*
* If SimpleCCW(a,b,c), then !SimpleCCW(c,b,a) for all a,b,c.
*
* In other words, ABC and CBA are guaranteed not to be both CCW
*/
public static boolean simpleCCW(S2Point a, S2Point b, S2Point c) {
// We compute the signed volume of the parallelepiped ABC. The usual
// formula for this is (AxB).C, but we compute it here using (CxA).B
// in order to ensure that ABC and CBA are not both CCW. This follows
// from the following identities (which are true numerically, not just
// mathematically):
//
// (1) x.CrossProd(y) == -(y.CrossProd(x))
// (2) (-x).DotProd(y) == -(x.DotProd(y))
return S2Point.crossProd(c, a).dotProd(b) > 0;
}
/**
* WARNING! This requires arbitrary precision arithmetic to be truly robust.
* This means that for nearly colinear AB and AC, this function may return the
* wrong answer.
*
*
* Like SimpleCCW(), but returns +1 if the points are counterclockwise and -1
* if the points are clockwise. It satisfies the following conditions:
*
* (1) RobustCCW(a,b,c) == 0 if and only if a == b, b == c, or c == a (2)
* RobustCCW(b,c,a) == RobustCCW(a,b,c) for all a,b,c (3) RobustCCW(c,b,a)
* ==-RobustCCW(a,b,c) for all a,b,c
*
* In other words:
*
* (1) The result is zero if and only if two points are the same. (2)
* Rotating the order of the arguments does not affect the result. (3)
* Exchanging any two arguments inverts the result.
*
* This function is essentially like taking the sign of the determinant of
* a,b,c, except that it has additional logic to make sure that the above
* properties hold even when the three points are coplanar, and to deal with
* the limitations of floating-point arithmetic.
*
* Note: a, b and c are expected to be of unit length. Otherwise, the results
* are undefined.
*/
public static int robustCCW(S2Point a, S2Point b, S2Point c) {
return robustCCW(a, b, c, S2Point.crossProd(a, b));
}
/**
* A more efficient version of RobustCCW that allows the precomputed
* cross-product of A and B to be specified.
*
* Note: a, b and c are expected to be of unit length. Otherwise, the results
* are undefined
*/
public static int robustCCW(S2Point a, S2Point b, S2Point c, S2Point aCrossB) {
// assert (isUnitLength(a) && isUnitLength(b) && isUnitLength(c));
// There are 14 multiplications and additions to compute the determinant
// below. Since all three points are normalized, it is possible to show
// that the average rounding error per operation does not exceed 2**-54,
// the maximum rounding error for an operation whose result magnitude is in
// the range [0.5,1). Therefore, if the absolute value of the determinant
// is greater than 2*14*(2**-54), the determinant will have the same sign
// even if the arguments are rotated (which produces a mathematically
// equivalent result but with potentially different rounding errors).
final double kMinAbsValue = 1.6e-15; // 2 * 14 * 2**-54
double det = aCrossB.dotProd(c);
// Double-check borderline cases in debug mode.
// assert ((Math.abs(det) < kMinAbsValue) || (Math.abs(det) > 1000 * kMinAbsValue)
// || (det * expensiveCCW(a, b, c) > 0));
if (det > kMinAbsValue) {
return 1;
}
if (det < -kMinAbsValue) {
return -1;
}
return expensiveCCW(a, b, c);
}
/**
* A relatively expensive calculation invoked by RobustCCW() if the sign of
* the determinant is uncertain.
*/
private static int expensiveCCW(S2Point a, S2Point b, S2Point c) {
// Return zero if and only if two points are the same. This ensures (1).
if (a.equals(b) || b.equals(c) || c.equals(a)) {
return 0;
}
// Now compute the determinant in a stable way. Since all three points are
// unit length and we know that the determinant is very close to zero, this
// means that points are very nearly colinear. Furthermore, the most common
// situation is where two points are nearly identical or nearly antipodal.
// To get the best accuracy in this situation, it is important to
// immediately reduce the magnitude of the arguments by computing either
// A+B or A-B for each pair of points. Note that even if A and B differ
// only in their low bits, A-B can be computed very accurately. On the
// other hand we can't accurately represent an arbitrary linear combination
// of two vectors as would be required for Gaussian elimination. The code
// below chooses the vertex opposite the longest edge as the "origin" for
// the calculation, and computes the different vectors to the other two
// vertices. This minimizes the sum of the lengths of these vectors.
//
// This implementation is very stable numerically, but it still does not
// return consistent results in all cases. For example, if three points are
// spaced far apart from each other along a great circle, the sign of the
// result will basically be random (although it will still satisfy the
// conditions documented in the header file). The only way to return
// consistent results in all cases is to compute the result using
// arbitrary-precision arithmetic. I considered using the Gnu MP library,
// but this would be very expensive (up to 2000 bits of precision may be
// needed to store the intermediate results) and seems like overkill for
// this problem. The MP library is apparently also quite particular about
// compilers and compilation options and would be a pain to maintain.
// We want to handle the case of nearby points and nearly antipodal points
// accurately, so determine whether A+B or A-B is smaller in each case.
double sab = (a.dotProd(b) > 0) ? -1 : 1;
double sbc = (b.dotProd(c) > 0) ? -1 : 1;
double sca = (c.dotProd(a) > 0) ? -1 : 1;
S2Point vab = S2Point.add(a, S2Point.mul(b, sab));
S2Point vbc = S2Point.add(b, S2Point.mul(c, sbc));
S2Point vca = S2Point.add(c, S2Point.mul(a, sca));
double dab = vab.norm2();
double dbc = vbc.norm2();
double dca = vca.norm2();
// Sort the difference vectors to find the longest edge, and use the
// opposite vertex as the origin. If two difference vectors are the same
// length, we break ties deterministically to ensure that the symmetry
// properties guaranteed in the header file will be true.
double sign;
if (dca < dbc || (dca == dbc && a.lessThan(b))) {
if (dab < dbc || (dab == dbc && a.lessThan(c))) {
// The "sab" factor converts A +/- B into B +/- A.
sign = S2Point.crossProd(vab, vca).dotProd(a) * sab; // BC is longest
// edge
} else {
sign = S2Point.crossProd(vca, vbc).dotProd(c) * sca; // AB is longest
// edge
}
} else {
if (dab < dca || (dab == dca && b.lessThan(c))) {
sign = S2Point.crossProd(vbc, vab).dotProd(b) * sbc; // CA is longest
// edge
} else {
sign = S2Point.crossProd(vca, vbc).dotProd(c) * sca; // AB is longest
// edge
}
}
if (sign > 0) {
return 1;
}
if (sign < 0) {
return -1;
}
// The points A, B, and C are numerically indistinguishable from coplanar.
// This may be due to roundoff error, or the points may in fact be exactly
// coplanar. We handle this situation by perturbing all of the points by a
// vector (eps, eps**2, eps**3) where "eps" is an infinitesmally small
// positive number (e.g. 1 divided by a googolplex). The perturbation is
// done symbolically, i.e. we compute what would happen if the points were
// perturbed by this amount. It turns out that this is equivalent to
// checking whether the points are ordered CCW around the origin first in
// the Y-Z plane, then in the Z-X plane, and then in the X-Y plane.
int ccw =
planarOrderedCCW(new R2Vector(a.y, a.z), new R2Vector(b.y, b.z), new R2Vector(c.y, c.z));
if (ccw == 0) {
ccw =
planarOrderedCCW(new R2Vector(a.z, a.x), new R2Vector(b.z, b.x), new R2Vector(c.z, c.x));
if (ccw == 0) {
ccw = planarOrderedCCW(
new R2Vector(a.x, a.y), new R2Vector(b.x, b.y), new R2Vector(c.x, c.y));
// assert (ccw != 0);
}
}
return ccw;
}
public static int planarCCW(R2Vector a, R2Vector b) {
// Return +1 if the edge AB is CCW around the origin, etc.
double sab = (a.dotProd(b) > 0) ? -1 : 1;
R2Vector vab = R2Vector.add(a, R2Vector.mul(b, sab));
double da = a.norm2();
double db = b.norm2();
double sign;
if (da < db || (da == db && a.lessThan(b))) {
sign = a.crossProd(vab) * sab;
} else {
sign = vab.crossProd(b);
}
if (sign > 0) {
return 1;
}
if (sign < 0) {
return -1;
}
return 0;
}
public static int planarOrderedCCW(R2Vector a, R2Vector b, R2Vector c) {
int sum = 0;
sum += planarCCW(a, b);
sum += planarCCW(b, c);
sum += planarCCW(c, a);
if (sum > 0) {
return 1;
}
if (sum < 0) {
return -1;
}
return 0;
}
/**
* Return true if the edges OA, OB, and OC are encountered in that order while
* sweeping CCW around the point O. You can think of this as testing whether
* A <= B <= C with respect to a continuous CCW ordering around O.
*
* Properties:
*
* - If orderedCCW(a,b,c,o) && orderedCCW(b,a,c,o), then a == b
* - If orderedCCW(a,b,c,o) && orderedCCW(a,c,b,o), then b == c
* - If orderedCCW(a,b,c,o) && orderedCCW(c,b,a,o), then a == b == c
* - If a == b or b == c, then orderedCCW(a,b,c,o) is true
* - Otherwise if a == c, then orderedCCW(a,b,c,o) is false
*
*/
public static boolean orderedCCW(S2Point a, S2Point b, S2Point c, S2Point o) {
// The last inequality below is ">" rather than ">=" so that we return true
// if A == B or B == C, and otherwise false if A == C. Recall that
// RobustCCW(x,y,z) == -RobustCCW(z,y,x) for all x,y,z.
int sum = 0;
if (robustCCW(b, o, a) >= 0) {
++sum;
}
if (robustCCW(c, o, b) >= 0) {
++sum;
}
if (robustCCW(a, o, c) > 0) {
++sum;
}
return sum >= 2;
}
/**
* Return the angle at the vertex B in the triangle ABC. The return value is
* always in the range [0, Pi]. The points do not need to be normalized.
* Ensures that Angle(a,b,c) == Angle(c,b,a) for all a,b,c.
*
* The angle is undefined if A or C is diametrically opposite from B, and
* becomes numerically unstable as the length of edge AB or BC approaches 180
* degrees.
*/
public static double angle(S2Point a, S2Point b, S2Point c) {
return S2Point.crossProd(a, b).angle(S2Point.crossProd(c, b));
}
/**
* Return the exterior angle at the vertex B in the triangle ABC. The return
* value is positive if ABC is counterclockwise and negative otherwise. If you
* imagine an ant walking from A to B to C, this is the angle that the ant
* turns at vertex B (positive = left, negative = right). Ensures that
* TurnAngle(a,b,c) == -TurnAngle(c,b,a) for all a,b,c.
*
* @param a
* @param b
* @param c
* @return the exterior angle at the vertex B in the triangle ABC
*/
public static double turnAngle(S2Point a, S2Point b, S2Point c) {
// This is a bit less efficient because we compute all 3 cross products, but
// it ensures that turnAngle(a,b,c) == -turnAngle(c,b,a) for all a,b,c.
double outAngle = S2Point.crossProd(b, a).angle(S2Point.crossProd(c, b));
return (robustCCW(a, b, c) > 0) ? outAngle : -outAngle;
}
/**
* Return true if two points are within the given distance of each other
* (mainly useful for testing).
*/
public static boolean approxEquals(S2Point a, S2Point b, double maxError) {
return a.angle(b) <= maxError;
}
public static boolean approxEquals(S2Point a, S2Point b) {
return approxEquals(a, b, 1e-15);
}
public static boolean approxEquals(double a, double b, double maxError) {
return Math.abs(a - b) <= maxError;
}
public static boolean approxEquals(double a, double b) {
return approxEquals(a, b, 1e-15);
}
// Don't instantiate
private S2() {
}
}