diff --git a/s2-geometry.iml b/s2-geometry.iml
new file mode 100644
index 00000000..8abc0aea
--- /dev/null
+++ b/s2-geometry.iml
@@ -0,0 +1,23 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/src/com/google/common/geometry/S2PolylineSimplifier.java b/src/com/google/common/geometry/S2PolylineSimplifier.java
new file mode 100644
index 00000000..cb14a599
--- /dev/null
+++ b/src/com/google/common/geometry/S2PolylineSimplifier.java
@@ -0,0 +1,258 @@
+package com.google.common.geometry;
+
+import static com.google.common.geometry.S2.M_PI;
+
+/**
+ * This is a helper class for simplifying polylines. It allows you to compute a maximal edge that intersects a sequence
+ * of discs, and that optionally avoids a different sequence of discs. The results are conservative in that the edge
+ * is guaranteed to intersect or avoid the specified discs using exact arithmetic (see S2Predicates).
+ *
+ * Note that S2Builder can also simplify polylines and supports more features (e.g., snapping to S2CellId centers),
+ * so it is only recommended to use this class if S2Builder does not meet your needs.
+ *
+ * Here is a simple example showing how to simplify a polyline into a sequence of edges that stay within "max_error"
+ * of the original edges:
+ *
+ *
+ * List v = ...
+ * List simplifiedPolyline = new ArrayList();
+ * S2PolylineSimplifier simplifier = new S2PolylineSimplifier();
+ * simplifier.init(v.get(0));
+ * simplifiedPolyline.add(v.get(0));
+ * for (int i = 1 ; i < v.size() ; ++i) {
+ * if (simplifier.extend(v[i])) {
+ * simplifiedPolyline.add(v[i]);
+ * simplifier.init(v[i])
+ * simplifier.targetDisc(v[i], max_error)
+ * }
+ * }
+ *
+ *
+ * Note that the points targeted by TargetDisc do not need to be the same as the candidate endpoints passed to extend.
+ * So for example, you could target the original vertices of a polyline, but only consider endpoints that are snapped to
+ * E7 coordinates or S2CellId centers.
+ *
+ * Please be aware that this class works by maintaining a range of acceptable angles (bearings) from the start vertex
+ * to the hypothetical destination vertex. It does not keep track of distances to any of the discs to be targeted or
+ * avoided. Therefore to use this class correctly, constraints should be added in increasing order of distance.
+ * (The actual requirement is slightly weaker than this, which is why it is not enforced, but basically you should only
+ * call targetDisc() and avoidDisc() with arguments that you want to constrain the immediately following call to
+ * extend().)
+ *
+ */
+class S2PolylineSimplifier {
+
+ private S2Point src;
+ private S2Point xDir = new S2Point();
+ private S2Point yDir = new S2Point();
+ private S1Interval window;
+
+ /**
+ * Starts a new simplified edge at "src".
+ *
+ * @param src The source point.
+ */
+ public void init(S2Point src) {
+ this.src = src;
+ this.window = S1Interval.full();
+
+ // Precompute basis vectors for the tangent space at "src". This is similar
+ // to GetFrame() except that we don't normalize the vectors. As it turns
+ // out, the two basis vectors below have the same magnitude (up to the
+ // length error in S2Point::Normalize).
+
+ int i, j, k;
+ // Find the index of the component whose magnitude is smallest.
+ S2Point tmp = src.fabs();
+ if (tmp.x < tmp.y) {
+ if (tmp.x < tmp.z) i = 0;
+ else i = 2;
+ } else if (tmp.y < tmp.z) i = 1;
+ else i = 2;
+
+ // We define the "y" basis vector as the cross product of "src" and the
+ // basis vector for axis "i". Let "j" and "k" be the indices of the other
+ // two components in cyclic order.
+ if (i == 2) j = 0;
+ else j = (i + 1);
+ if (i == 0) k = 2;
+ else k = (i - 1);
+
+ double[] yCoords = new double[3];
+ yCoords[i] = 0.0;
+ yCoords[j] = src.get(k);
+ yCoords[k] = -src.get(j);
+ yDir = new S2Point(yCoords[0], yCoords[1], yCoords[2]);
+
+ // Compute the cross product of "y_dir" and "src". We write out the cross
+ // product here mainly for documentation purposes; it also happens to save a
+ // few multiplies because unfortunately the optimizer does *not* get rid of
+ // multiplies by zero (since these multiplies propagate NaN, for example).
+ double[] xCoords = new double[3];
+ xCoords[i] = src.get(j) * src.get(j) + src.get(k) * src.get(k);
+ xCoords[j] = -src.get(j) * src.get(i);
+ xCoords[k] = -src.get(k) * src.get(i);
+ xDir = new S2Point(xCoords[0], xCoords[1], xCoords[2]);
+ }
+
+ /** Returns the source vertex of the output edge. */
+ private S2Point getSrc() {
+ return src;
+ }
+
+ /**
+ * Returns true if the edge (src, dst) satisfies all of the targeting requirements so far.
+ * Returns false if the edge would be longer than 90 degrees (such edges are not supported).
+ */
+ public boolean extend(S2Point dst) {
+ // We limit the maximum edge length to 90 degrees in order to simplify the
+ // error bounds. (The error gets arbitrarily large as the edge length
+ // approaches 180 degrees.)
+ if (new S1ChordAngle(src, dst).getLength2() > S1ChordAngle.RIGHT.getLength2()) return false;
+
+ // Otherwise check whether this vertex is in the acceptable angle range.
+ return window.contains(getAngle(dst));
+ }
+
+ /** Requires that the output edge must pass through the given disc. */
+ public boolean targetDisc(S2Point point, S1ChordAngle radius) {
+ // Shrink the target interval by the maximum error from all sources. This
+ // guarantees that the output edge will intersect the given disc.
+ double semiwidth = getSemiwidth(point, radius, -1 /*round down*/);
+ if (semiwidth >= M_PI) {
+ // The target disc contains "src", so there is nothing to do.
+ return true;
+ }
+ if (semiwidth < 0) {
+ window = S1Interval.empty();
+ return false;
+ }
+ // Otherwise compute the angle interval corresponding to the target disc and
+ // intersect it with the current window.
+ double center = getAngle(point);
+ S1Interval target = S1Interval.fromPoint(center).expanded(semiwidth);
+ window = window.intersection(target);
+ return !window.isEmpty();
+ }
+
+ /**
+ * Requires that the output edge must avoid the given disc. "disc_on_left"
+ * specifies whether the disc must be to the left or right of the edge.
+ * (This feature allows the simplified edge to preserve the topology of the
+ * original polyline with respect to other nearby points.)
+ *
+ * If your input is a polyline, you can compute "disc_on_left" as follows.
+ * Let the polyline be ABCDE and assume that it already avoids a set of
+ * points X_i. Suppose that you have already added ABC to the simplifer, and
+ * now want to extend the edge chain to D. First find the X_i that are near
+ * the edge CD, then discard the ones such that AX_i <= AC or AX_i >= AD
+ * (since these points have either already been considered or aren't
+ * relevant yet). Now X_i is to the left of the polyline if and only if
+ * s2pred::OrderedCCW(A, D, X, C) (in other words, if X_i is to the left of
+ * the angle wedge ACD).
+ */
+ public boolean avoidDisc(S2Point point, S1ChordAngle radius, boolean discOnLeft) {
+ // Expand the interval by the maximum error from all sources. This
+ // guarantees that the final output edge will avoid the given disc.
+ double semiwidth = getSemiwidth(point, radius, 1 /*round up*/);
+ if (semiwidth >= M_PI) {
+ // The avoidance disc contains "src", so it is impossible to avoid.
+ window = S1Interval.empty();
+ return false;
+ }
+ double center = getAngle(point);
+ double opposite;
+ if (center > 0) opposite = center - M_PI;
+ else opposite = center + M_PI;
+ S1Interval target;
+ if (discOnLeft) target = new S1Interval(opposite, center);
+ else target = new S1Interval(center, opposite);
+ window = window.intersection(target.expanded(-semiwidth));
+ return !window.isEmpty();
+ }
+
+ private double getAngle(S2Point p) {
+ return Math.atan2(p.dotProd(yDir), p.dotProd(xDir));
+ }
+
+ private double getSemiwidth(S2Point p, S1ChordAngle r, int roundDirection) {
+ double err = 0.5 * S2.DBL_EPSILON;
+
+ // Using spherical trigonometry,
+ //
+ // sin(semiwidth) = sin(r) / sin(a)
+ //
+ // where "a" is the angle between "src" and "p". Rather than measuring
+ // these angles, instead we measure the squared chord lengths through the
+ // interior of the sphere (i.e., Cartersian distance). Letting "r2" be the
+ // squared chord distance corresponding to "r", and "a2" be the squared
+ // chord distance corresponding to "a", we use the relationships
+ //
+ // sin^2(r) = r2 (1 - r2 / 4)
+ // sin^2(a) = d2 (1 - d2 / 4)
+ //
+ // which follow from the fact that r2 = (2 * sin(r / 2)) ^ 2, etc.
+
+ // "a2" has a relative error up to 5 * DBL_ERR, plus an absolute error of up
+ // to 64 * DBL_ERR^2 (because "src" and "p" may differ from unit length by
+ // up to 4 * DBL_ERR). We can correct for the relative error later, but for
+ // the absolute error we use "round_direction" to account for it now.
+ double r2 = r.getLength2();
+ double a2 = new S1ChordAngle(src, p).getLength2();
+ a2 -= 64 * err * err * roundDirection;
+ if (a2 <= r2) return M_PI; // The given disc contains "src".
+
+ double sin2_r = r2 * (1 - 0.25 * r2);
+ double sin2_a = a2 * (1 - 0.25 * a2);
+ double semiwidth = Math.asin(Math.sqrt(sin2_r / sin2_a));
+
+ // We compute bounds on the errors from all sources:
+ //
+ // - The call to GetSemiwidth (this call).
+ // - The call to GetAngle that computes the center of the interval.
+ // - The call to GetAngle in Extend that tests whether a given point
+ // is an acceptable destination vertex.
+ //
+ // Summary of the errors in GetAngle:
+ //
+ // y_dir_ has no error.
+ //
+ // x_dir_ has a relative error of DBL_ERR in two components, a relative
+ // error of 2 * DBL_ERR in the other component, plus an overall relative
+ // length error of 4 * DBL_ERR (compared to y_dir_) because "src" is assumed
+ // to be normalized only to within the tolerances of S2Point::Normalize().
+ //
+ // p.DotProd(y_dir_) has a relative error of 1.5 * DBL_ERR and an
+ // absolute error of 1.5 * DBL_ERR * y_dir_.Norm().
+ //
+ // p.DotProd(x_dir_) has a relative error of 5.5 * DBL_ERR and an absolute
+ // error of 3.5 * DBL_ERR * y_dir_.Norm() (noting that x_dir_ and y_dir_
+ // have the same length to within a relative error of 4 * DBL_ERR).
+ //
+ // It's possible to show by taking derivatives that these errors can affect
+ // the angle atan2(y, x) by up 7.093 * DBL_ERR radians. Rounding up and
+ // including the call to atan2 gives a final error bound of 10 * DBL_ERR.
+ //
+ // Summary of the errors in GetSemiwidth:
+ //
+ // The distance a2 has a relative error of 5 * DBL_ERR plus an absolute
+ // error of 64 * DBL_ERR^2 because the points "src" and "p" may differ from
+ // unit length (by up to 4 * DBL_ERR). We have already accounted for the
+ // absolute error above, leaving only the relative error.
+ //
+ // sin2_r has a relative error of 2 * DBL_ERR.
+ //
+ // sin2_a has a relative error of 12 * DBL_ERR assuming that a2 <= 2,
+ // i.e. distance(src, p) <= 90 degrees. (The relative error gets
+ // arbitrarily larger as this distance approaches 180 degrees.)
+ //
+ // semiwidth has a relative error of 17 * DBL_ERR.
+ //
+ // Finally, (center +/- semiwidth) has a rounding error of up to 4 * DBL_ERR
+ // because in theory, the result magnitude may be as large as 1.5 * M_PI
+ // which is larger than 4.0. This gives a total error of:
+ double error = (2 * 10 + 4) * err + 17 * err * semiwidth;
+ return semiwidth + roundDirection * error;
+ }
+
+}
diff --git a/tests/com/google/common/geometry/S2PolylineSimplifierTest.java b/tests/com/google/common/geometry/S2PolylineSimplifierTest.java
new file mode 100644
index 00000000..1171712b
--- /dev/null
+++ b/tests/com/google/common/geometry/S2PolylineSimplifierTest.java
@@ -0,0 +1,162 @@
+package com.google.common.geometry;
+
+import com.google.common.collect.Lists;
+
+import java.util.ArrayList;
+import java.util.List;
+
+public class S2PolylineSimplifierTest extends GeometryTestCase {
+
+ private void checkSimplify(String src, String dst, String target, String avoid,
+ List disc_on_left, double radiusDegrees, boolean expectedResult) {
+ S1ChordAngle radius = S1ChordAngle.fromS1Angle(S1Angle.degrees(radiusDegrees));
+ S2PolylineSimplifier s = new S2PolylineSimplifier();
+ s.init(makePoint(src));
+ List targets = new ArrayList<>();
+ parseVertices(target, targets);
+ for (S2Point p : targets) {
+ s.targetDisc(p, radius);
+ }
+ int i = 0;
+ List avoids = new ArrayList<>();
+ parseVertices(avoid, avoids);
+ for (S2Point p : avoids) {
+ s.avoidDisc(p, radius, disc_on_left.get(i++));
+ }
+ assertEquals(expectedResult, s.extend(makePoint(dst)));
+ }
+
+ public void testReuse() {
+ // Check that Init() can be called more than once.
+ S2PolylineSimplifier s = new S2PolylineSimplifier();
+ S1ChordAngle radius = S1ChordAngle.fromS1Angle(S1Angle.degrees(10));
+ s.init(new S2Point(1, 0, 0));
+ assertTrue(s.targetDisc(new S2Point(1, 1, 0).normalize(), radius));
+ assertTrue(s.targetDisc(new S2Point(1.0, 1.0, 0.1).normalize(), radius));
+ assertFalse(s.extend(new S2Point(1.0, 1.0, 0.4).normalize()));
+
+ s.init(new S2Point(0, 1, 0));
+ assertTrue(s.targetDisc(new S2Point(1.0, 1.0, 0.3).normalize(), radius));
+ assertTrue(s.targetDisc(new S2Point(1.0, 1.0, 0.2).normalize(), radius));
+ assertFalse(s.extend(new S2Point(1.0, 1.0, 0.0).normalize()));
+ }
+
+ public void testNoConstraints() {
+ // No constraints, dst == src.
+ checkSimplify("0:1", "0:1", "", "", new ArrayList<>(), 0.0, true);
+
+ // No constraints, dst != src.
+ checkSimplify("0:1", "1:0", "", "", new ArrayList<>(), 0.0, true);
+
+ // No constraints, (src, dst) longer than 90 degrees (not supported).
+ checkSimplify("0:0", "0:91", "", "", new ArrayList<>(), 0.0, false);
+ }
+
+ public void testTargetOnePoint() {
+ // Three points on a straight line. In theory zero tolerance should work,
+ // but in practice there are floating point errors.
+ checkSimplify("0:0", "0:2", "0:1", "", new ArrayList<>(), 1e-10, true);
+
+ // Three points where the middle point is too far away.
+ checkSimplify("0:0", "0:2", "1:1", "", new ArrayList<>(), 0.9, false);
+
+ // A target disc that contains the source vertex.
+ checkSimplify("0:0", "0:2", "0:0.1", "", new ArrayList<>(), 1.0, true);
+
+ // A target disc that contains the destination vertex.
+ checkSimplify("0:0", "0:2", "0:2.1", "", new ArrayList<>(), 1.0, true);
+ }
+
+ public void testAvoidOnePoint() {
+ // Three points on a straight line, attempting to avoid the middle point.
+ checkSimplify("0:0", "0:2", "", "0:1", Lists.newArrayList(true), 1e-10, false);
+ // Three points where the middle point can be successfully avoided.
+ checkSimplify("0:0", "0:2", "", "1:1", Lists.newArrayList(true), 0.9, true);
+ // Three points where the middle point is on the left, but where the client
+ // requires the point to be on the right of the edge.
+ checkSimplify("0:0", "0:2", "", "1:1", Lists.newArrayList(false), 1e-10, false);
+ }
+ /*
+yDir=(0.0, -1.0, 0.0)
+xDir=(-0.0, -0.0, 1.0)
+getAngle: p=(0.9996954135095479, 0.017449748351250485, 0.01745240643728351), xDir=(-0.0, -0.0, 1.0), yDir=(0.0, -1.0, 0.0)
+center=-0.7853220051761581
+opposite=2.356270648413635
+target=[2.356270648413635, -0.7853220051761581]
+window=[3.046158530131703, -1.4752098868942265]
+ */
+
+ public void testTargetAndAvoid() {
+ // Target several points that are separated from the proposed edge by about
+ // 0.7 degrees, and avoid several points that are separated from the
+ // proposed edge by about 1.4 degrees.
+ checkSimplify(
+ "0:0", "10:10", "2:3, 4:3, 7:8",
+ "4:2, 7:5, 7:9", Lists.newArrayList(true, true, false), 1.0, true
+ );
+
+ // The same example, but one point to be targeted is 1.4 degrees away.
+ checkSimplify(
+ "0:0", "10:10", "2:3, 4:6, 7:8",
+ "4:2, 7:5, 7:9", Lists.newArrayList(true, true, false), 1.0, false
+ );
+
+ // The same example, but one point to be avoided is 0.7 degrees away.
+ checkSimplify(
+ "0:0", "10:10", "2:3, 4:3, 7:8",
+ "4:2, 6:5, 7:9", Lists.newArrayList(true, true, false), 1.0, false
+ );
+ }
+
+ public void testPrecision() {
+ // This is a rough upper bound on both the error in constructing the disc
+ // locations (i.e., S2::InterpolateAtDistance, etc), and also on the
+ // padding that S2PolylineSimplifier uses to ensure that its results are
+ // conservative (i.e., the error calculated by GetSemiwidth).
+ S1Angle kMaxError = S1Angle.radians(25 * S2.DBL_EPSILON);
+
+ // We repeatedly generate a random edge. We then target several discs that
+ // barely overlap the edge, and avoid several discs that barely miss the
+ // edge. About half the time, we choose one disc and make it slightly too
+ // large or too small so that targeting fails.
+ int kIters = 1000; // Passes with 1 million iterations.
+ S2PolylineSimplifier simplifier = new S2PolylineSimplifier();
+
+ for (int iter = 0; iter < kIters; iter++) {
+ rand.setSeed(iter + 1); // Easier to reproduce a specific case.
+ S2Point src = randomPoint();
+ simplifier.init(src);
+ S2Point dst = S2EdgeUtil.interpolateAtDistance(
+ S1Angle.radians(rand.nextDouble()),
+ src,
+ randomPoint()
+ );
+ S2Point n = S2.robustCrossProd(src, dst).normalize();
+
+ // If bad_disc >= 0, then we make targeting fail for that disc.
+ int kNumDiscs = 5;
+ int bad_disc = rand.nextInt(2 * kNumDiscs) - kNumDiscs;
+ for (int i = 0; i < kNumDiscs; i++) {
+ double f = rand.nextDouble();
+ S2Point a = ((src.mul(1 - f)).add(dst.mul(f))).normalize();
+ S1Angle r = S1Angle.radians(rand.nextDouble());
+ boolean on_left = oneIn(2);
+ S2Point x = S2EdgeUtil.interpolateAtDistance(r, a, on_left ? n : n.neg());
+ // We grow the radius slightly if we want to target the disc and shrink
+ // it otherwise, *unless* we want targeting to fail for this disc, in
+ // which case these actions are reversed.
+ boolean avoid = oneIn(2);
+ boolean grow_radius = (avoid == (i == bad_disc));
+ S1ChordAngle radius = S1ChordAngle.fromS1Angle(grow_radius ? (r.add(kMaxError)) : (r.sub(kMaxError)));
+ if (avoid) {
+ simplifier.avoidDisc(x, radius, on_left);
+ } else {
+ simplifier.targetDisc(x, radius);
+ }
+ }
+ // The result is true iff all the disc constraints were satisfiable.
+ assertEquals(bad_disc < 0, simplifier.extend(dst));
+ }
+ }
+
+}