| /* |
| * Copyright 2006 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.base.Preconditions; |
| |
| /** |
| * This class contains various utility functions related to edges. It collects |
| * together common code that is needed to implement polygonal geometry such as |
| * polylines, loops, and general polygons. |
| * |
| */ |
| public strictfp class S2EdgeUtil { |
| /** |
| * IEEE floating-point operations have a maximum error of 0.5 ULPS (units in |
| * the last place). For double-precision numbers, this works out to 2**-53 |
| * (about 1.11e-16) times the magnitude of the result. It is possible to |
| * analyze the calculation done by getIntersection() and work out the |
| * worst-case rounding error. I have done a rough version of this, and my |
| * estimate is that the worst case distance from the intersection point X to |
| * the great circle through (a0, a1) is about 12 ULPS, or about 1.3e-15. This |
| * needs to be increased by a factor of (1/0.866) to account for the |
| * edgeSpliceFraction() in S2PolygonBuilder. Note that the maximum error |
| * measured by the unittest in 1,000,000 trials is less than 3e-16. |
| */ |
| public static final S1Angle DEFAULT_INTERSECTION_TOLERANCE = S1Angle.radians(1.5e-15); |
| |
| /** |
| * This class allows a vertex chain v0, v1, v2, ... to be efficiently tested |
| * for intersection with a given fixed edge AB. |
| */ |
| public static class EdgeCrosser { |
| // The fields below are all constant. |
| |
| private final S2Point a; |
| private final S2Point b; |
| private final S2Point aCrossB; |
| |
| // The fields below are updated for each vertex in the chain. |
| |
| // Previous vertex in the vertex chain. |
| private S2Point c; |
| // The orientation of the triangle ACB. |
| private int acb; |
| |
| /** |
| * AB is the given fixed edge, and C is the first vertex of the vertex |
| * chain. All parameters must point to fixed storage that persists for the |
| * lifetime of the EdgeCrosser object. |
| */ |
| public EdgeCrosser(S2Point a, S2Point b, S2Point c) { |
| this.a = a; |
| this.b = b; |
| this.aCrossB = S2Point.crossProd(a, b); |
| restartAt(c); |
| } |
| |
| /** |
| * Call this function when your chain 'jumps' to a new place. |
| */ |
| public void restartAt(S2Point c) { |
| this.c = c; |
| acb = -S2.robustCCW(a, b, c, aCrossB); |
| } |
| |
| /** |
| * This method is equivalent to calling the S2EdgeUtil.robustCrossing() |
| * function (defined below) on the edges AB and CD. It returns +1 if there |
| * is a crossing, -1 if there is no crossing, and 0 if two points from |
| * different edges are the same. Returns 0 or -1 if either edge is |
| * degenerate. As a side effect, it saves vertex D to be used as the next |
| * vertex C. |
| */ |
| public int robustCrossing(S2Point d) { |
| // For there to be an edge crossing, the triangles ACB, CBD, BDA, DAC must |
| // all be oriented the same way (CW or CCW). We keep the orientation |
| // of ACB as part of our state. When each new point D arrives, we |
| // compute the orientation of BDA and check whether it matches ACB. |
| // This checks whether the points C and D are on opposite sides of the |
| // great circle through AB. |
| |
| // Recall that robustCCW is invariant with respect to rotating its |
| // arguments, i.e. ABC has the same orientation as BDA. |
| int bda = S2.robustCCW(a, b, d, aCrossB); |
| int result; |
| |
| if (bda == -acb && bda != 0) { |
| // Most common case -- triangles have opposite orientations. |
| result = -1; |
| } else if ((bda & acb) == 0) { |
| // At least one value is zero -- two vertices are identical. |
| result = 0; |
| } else { |
| // assert (bda == acb && bda != 0); |
| result = robustCrossingInternal(d); // Slow path. |
| } |
| // Now save the current vertex D as the next vertex C, and also save the |
| // orientation of the new triangle ACB (which is opposite to the current |
| // triangle BDA). |
| c = d; |
| acb = -bda; |
| return result; |
| } |
| |
| /** |
| * This method is equivalent to the S2EdgeUtil.edgeOrVertexCrossing() method |
| * defined below. It is similar to robustCrossing, but handles cases where |
| * two vertices are identical in a way that makes it easy to implement |
| * point-in-polygon containment tests. |
| */ |
| public boolean edgeOrVertexCrossing(S2Point d) { |
| // We need to copy c since it is clobbered by robustCrossing(). |
| S2Point c2 = new S2Point(c.get(0), c.get(1), c.get(2)); |
| |
| int crossing = robustCrossing(d); |
| if (crossing < 0) { |
| return false; |
| } |
| if (crossing > 0) { |
| return true; |
| } |
| |
| return vertexCrossing(a, b, c2, d); |
| } |
| |
| /** |
| * This function handles the "slow path" of robustCrossing(). |
| */ |
| private int robustCrossingInternal(S2Point d) { |
| // ACB and BDA have the appropriate orientations, so now we check the |
| // triangles CBD and DAC. |
| S2Point cCrossD = S2Point.crossProd(c, d); |
| int cbd = -S2.robustCCW(c, d, b, cCrossD); |
| if (cbd != acb) { |
| return -1; |
| } |
| |
| int dac = S2.robustCCW(c, d, a, cCrossD); |
| return (dac == acb) ? 1 : -1; |
| } |
| } |
| |
| /** |
| * This class computes a bounding rectangle that contains all edges defined by |
| * a vertex chain v0, v1, v2, ... All vertices must be unit length. Note that |
| * the bounding rectangle of an edge can be larger than the bounding rectangle |
| * of its endpoints, e.g. consider an edge that passes through the north pole. |
| */ |
| public static class RectBounder { |
| // The previous vertex in the chain. |
| private S2Point a; |
| |
| // The corresponding latitude-longitude. |
| private S2LatLng aLatLng; |
| |
| // The current bounding rectangle. |
| private S2LatLngRect bound; |
| |
| public RectBounder() { |
| this.bound = S2LatLngRect.empty(); |
| } |
| |
| /** |
| * This method is called to add each vertex to the chain. 'b' must point to |
| * fixed storage that persists for the lifetime of the RectBounder. |
| */ |
| public void addPoint(S2Point b) { |
| // assert (S2.isUnitLength(b)); |
| |
| S2LatLng bLatLng = new S2LatLng(b); |
| |
| if (bound.isEmpty()) { |
| bound = bound.addPoint(bLatLng); |
| } else { |
| // We can't just call bound.addPoint(bLatLng) here, since we need to |
| // ensure that all the longitudes between "a" and "b" are included. |
| bound = bound.union(S2LatLngRect.fromPointPair(aLatLng, bLatLng)); |
| |
| // Check whether the min/max latitude occurs in the edge interior. |
| // We find the normal to the plane containing AB, and then a vector |
| // "dir" in this plane that also passes through the equator. We use |
| // RobustCrossProd to ensure that the edge normal is accurate even |
| // when the two points are very close together. |
| S2Point aCrossB = S2.robustCrossProd(a, b); |
| S2Point dir = S2Point.crossProd(aCrossB, new S2Point(0, 0, 1)); |
| double da = dir.dotProd(a); |
| double db = dir.dotProd(b); |
| |
| if (da * db < 0) { |
| // Minimum/maximum latitude occurs in the edge interior. This affects |
| // the latitude bounds but not the longitude bounds. |
| double absLat = Math.acos(Math.abs(aCrossB.get(2) / aCrossB.norm())); |
| R1Interval lat = bound.lat(); |
| if (da < 0) { |
| // It's possible that absLat < lat.lo() due to numerical errors. |
| lat = new R1Interval(lat.lo(), Math.max(absLat, bound.lat().hi())); |
| } else { |
| lat = new R1Interval(Math.min(-absLat, bound.lat().lo()), lat.hi()); |
| } |
| bound = new S2LatLngRect(lat, bound.lng()); |
| } |
| } |
| a = b; |
| aLatLng = bLatLng; |
| } |
| |
| /** |
| * Return the bounding rectangle of the edge chain that connects the |
| * vertices defined so far. |
| */ |
| public S2LatLngRect getBound() { |
| return bound; |
| } |
| |
| } |
| |
| /** |
| * The purpose of this class is to find edges that intersect a given XYZ |
| * bounding box. It can be used as an efficient rejection test when attempting to |
| * find edges that intersect a given region. It accepts a vertex chain v0, v1, |
| * v2, ... and returns a boolean value indicating whether each edge intersects |
| * the specified bounding box. |
| * |
| * We use XYZ intervals instead of something like longitude intervals because |
| * it is cheap to collect from S2Point lists and any slicing strategy should |
| * give essentially equivalent results. See S2Loop for an example of use. |
| */ |
| public static class XYZPruner { |
| private S2Point lastVertex; |
| |
| // The region to be tested against. |
| private boolean boundSet; |
| private double xmin; |
| private double ymin; |
| private double zmin; |
| private double xmax; |
| private double ymax; |
| private double zmax; |
| private double maxDeformation; |
| |
| public XYZPruner() { |
| boundSet = false; |
| } |
| |
| /** |
| * Accumulate a bounding rectangle from provided edges. |
| * |
| * @param from start of edge |
| * @param to end of edge. |
| */ |
| public void addEdgeToBounds(S2Point from, S2Point to) { |
| if (!boundSet) { |
| boundSet = true; |
| xmin = xmax = from.x; |
| ymin = ymax = from.y; |
| zmin = zmax = from.z; |
| } |
| xmin = Math.min(xmin, Math.min(to.x, from.x)); |
| ymin = Math.min(ymin, Math.min(to.y, from.y)); |
| zmin = Math.min(zmin, Math.min(to.z, from.z)); |
| xmax = Math.max(xmax, Math.max(to.x, from.x)); |
| ymax = Math.max(ymax, Math.max(to.y, from.y)); |
| zmax = Math.max(zmax, Math.max(to.z, from.z)); |
| |
| // Because our arcs are really geodesics on the surface of the earth |
| // an edge can have intermediate points outside the xyz bounds implicit |
| // in the end points. Based on the length of the arc we compute a |
| // generous bound for the maximum amount of deformation. For small edges |
| // it will be very small but for some large arcs (ie. from (1N,90W) to |
| // (1N,90E) the path can be wildly deformed. I did a bunch of |
| // experiments with geodesics to get safe bounds for the deformation. |
| double approxArcLen = |
| Math.abs(from.x - to.x) + Math.abs(from.y - to.y) + Math.abs(from.z - to.z); |
| if (approxArcLen < 0.025) { // less than 2 degrees |
| maxDeformation = Math.max(maxDeformation, approxArcLen * 0.0025); |
| } else if (approxArcLen < 1.0) { // less than 90 degrees |
| maxDeformation = Math.max(maxDeformation, approxArcLen * 0.11); |
| } else { |
| maxDeformation = approxArcLen * 0.5; |
| } |
| } |
| |
| public void setFirstIntersectPoint(S2Point v0) { |
| xmin = xmin - maxDeformation; |
| ymin = ymin - maxDeformation; |
| zmin = zmin - maxDeformation; |
| xmax = xmax + maxDeformation; |
| ymax = ymax + maxDeformation; |
| zmax = zmax + maxDeformation; |
| this.lastVertex = v0; |
| } |
| |
| /** |
| * Returns true if the edge going from the last point to this point passes |
| * through the pruner bounding box, otherwise returns false. So the |
| * method returns false if we are certain there is no intersection, but it |
| * may return true when there turns out to be no intersection. |
| */ |
| public boolean intersects(S2Point v1) { |
| boolean result = true; |
| |
| if ((v1.x < xmin && lastVertex.x < xmin) || (v1.x > xmax && lastVertex.x > xmax)) { |
| result = false; |
| } else if ((v1.y < ymin && lastVertex.y < ymin) || (v1.y > ymax && lastVertex.y > ymax)) { |
| result = false; |
| } else if ((v1.z < zmin && lastVertex.z < zmin) || (v1.z > zmax && lastVertex.z > zmax)) { |
| result = false; |
| } |
| |
| lastVertex = v1; |
| return result; |
| } |
| } |
| |
| /** |
| * The purpose of this class is to find edges that intersect a given longitude |
| * interval. It can be used as an efficient rejection test when attempting to |
| * find edges that intersect a given region. It accepts a vertex chain v0, v1, |
| * v2, ... and returns a boolean value indicating whether each edge intersects |
| * the specified longitude interval. |
| * |
| * This class is not currently used as the XYZPruner is preferred for |
| * S2Loop, but this should be usable in similar circumstances. Be wary |
| * of the cost of atan2() in conversions from S2Point to longitude! |
| */ |
| public static class LongitudePruner { |
| // The interval to be tested against. |
| private S1Interval interval; |
| |
| // The longitude of the next v0. |
| private double lng0; |
| |
| /** |
| *'interval' is the longitude interval to be tested against, and 'v0' is |
| * the first vertex of edge chain. |
| */ |
| public LongitudePruner(S1Interval interval, S2Point v0) { |
| this.interval = interval; |
| this.lng0 = S2LatLng.longitude(v0).radians(); |
| } |
| |
| /** |
| * Returns true if the edge (v0, v1) intersects the given longitude |
| * interval, and then saves 'v1' to be used as the next 'v0'. |
| */ |
| public boolean intersects(S2Point v1) { |
| double lng1 = S2LatLng.longitude(v1).radians(); |
| boolean result = interval.intersects(S1Interval.fromPointPair(lng0, lng1)); |
| lng0 = lng1; |
| return result; |
| } |
| } |
| |
| /** |
| * A wedge relation's test method accepts two edge chains A=(a0,a1,a2) and |
| * B=(b0,b1,b2) where a1==b1, and returns either -1, 0, or 1 to indicate the |
| * relationship between the region to the left of A and the region to the left |
| * of B. Wedge relations are used to determine the local relationship between |
| * two polygons that share a common vertex. |
| * |
| * All wedge relations require that a0 != a2 and b0 != b2. Other degenerate |
| * cases (such as a0 == b2) are handled as expected. The parameter "ab1" |
| * denotes the common vertex a1 == b1. |
| */ |
| public interface WedgeRelation { |
| int test(S2Point a0, S2Point ab1, S2Point a2, S2Point b0, S2Point b2); |
| } |
| |
| public static class WedgeContains implements WedgeRelation { |
| /** |
| * Given two edge chains (see WedgeRelation above), this function returns +1 |
| * if the region to the left of A contains the region to the left of B, and |
| * 0 otherwise. |
| */ |
| @Override |
| public int test(S2Point a0, S2Point ab1, S2Point a2, S2Point b0, S2Point b2) { |
| // For A to contain B (where each loop interior is defined to be its left |
| // side), the CCW edge order around ab1 must be a2 b2 b0 a0. We split |
| // this test into two parts that test three vertices each. |
| return S2.orderedCCW(a2, b2, b0, ab1) && S2.orderedCCW(b0, a0, a2, ab1) ? 1 : 0; |
| } |
| } |
| |
| public static class WedgeIntersects implements WedgeRelation { |
| /** |
| * Given two edge chains (see WedgeRelation above), this function returns -1 |
| * if the region to the left of A intersects the region to the left of B, |
| * and 0 otherwise. Note that regions are defined such that points along a |
| * boundary are contained by one side or the other, not both. So for |
| * example, if A,B,C are distinct points ordered CCW around a vertex O, then |
| * the wedges BOA, AOC, and COB do not intersect. |
| */ |
| @Override |
| public int test(S2Point a0, S2Point ab1, S2Point a2, S2Point b0, S2Point b2) { |
| // For A not to intersect B (where each loop interior is defined to be |
| // its left side), the CCW edge order around ab1 must be a0 b2 b0 a2. |
| // Note that it's important to write these conditions as negatives |
| // (!OrderedCCW(a,b,c,o) rather than Ordered(c,b,a,o)) to get correct |
| // results when two vertices are the same. |
| return (S2.orderedCCW(a0, b2, b0, ab1) && S2.orderedCCW(b0, a2, a0, ab1) ? 0 : -1); |
| } |
| } |
| |
| public static class WedgeContainsOrIntersects implements WedgeRelation { |
| /** |
| * Given two edge chains (see WedgeRelation above), this function returns +1 |
| * if A contains B, 0 if A and B are disjoint, and -1 if A intersects but |
| * does not contain B. |
| */ |
| @Override |
| public int test(S2Point a0, S2Point ab1, S2Point a2, S2Point b0, S2Point b2) { |
| // This is similar to WedgeContainsOrCrosses, except that we want to |
| // distinguish cases (1) [A contains B], (3) [A and B are disjoint], |
| // and (2,4,5,6) [A intersects but does not contain B]. |
| |
| if (S2.orderedCCW(a0, a2, b2, ab1)) { |
| // We are in case 1, 5, or 6, or case 2 if a2 == b2. |
| return S2.orderedCCW(b2, b0, a0, ab1) ? 1 : -1; // Case 1 vs. 2,5,6. |
| } |
| // We are in cases 2, 3, or 4. |
| if (!S2.orderedCCW(a2, b0, b2, ab1)) { |
| return 0; // Case 3. |
| } |
| |
| // We are in case 2 or 4, or case 3 if a2 == b0. |
| return (a2.equals(b0)) ? 0 : -1; // Case 3 vs. 2,4. |
| } |
| } |
| |
| public static class WedgeContainsOrCrosses implements WedgeRelation { |
| /** |
| * Given two edge chains (see WedgeRelation above), this function returns +1 |
| * if A contains B, 0 if B contains A or the two wedges do not intersect, |
| * and -1 if the edge chains A and B cross each other (i.e. if A intersects |
| * both the interior and exterior of the region to the left of B). In |
| * degenerate cases where more than one of these conditions is satisfied, |
| * the maximum possible result is returned. For example, if A == B then the |
| * result is +1. |
| */ |
| @Override |
| public int test(S2Point a0, S2Point ab1, S2Point a2, S2Point b0, S2Point b2) { |
| // There are 6 possible edge orderings at a shared vertex (all |
| // of these orderings are circular, i.e. abcd == bcda): |
| // |
| // (1) a2 b2 b0 a0: A contains B |
| // (2) a2 a0 b0 b2: B contains A |
| // (3) a2 a0 b2 b0: A and B are disjoint |
| // (4) a2 b0 a0 b2: A and B intersect in one wedge |
| // (5) a2 b2 a0 b0: A and B intersect in one wedge |
| // (6) a2 b0 b2 a0: A and B intersect in two wedges |
| // |
| // In cases (4-6), the boundaries of A and B cross (i.e. the boundary |
| // of A intersects the interior and exterior of B and vice versa). |
| // Thus we want to distinguish cases (1), (2-3), and (4-6). |
| // |
| // Note that the vertices may satisfy more than one of the edge |
| // orderings above if two or more vertices are the same. The tests |
| // below are written so that we take the most favorable |
| // interpretation, i.e. preferring (1) over (2-3) over (4-6). In |
| // particular note that if orderedCCW(a,b,c,o) returns true, it may be |
| // possible that orderedCCW(c,b,a,o) is also true (if a == b or b == c). |
| |
| if (S2.orderedCCW(a0, a2, b2, ab1)) { |
| // The cases with this vertex ordering are 1, 5, and 6, |
| // although case 2 is also possible if a2 == b2. |
| if (S2.orderedCCW(b2, b0, a0, ab1)) { |
| return 1; // Case 1 (A contains B) |
| } |
| |
| // We are in case 5 or 6, or case 2 if a2 == b2. |
| return (a2.equals(b2)) ? 0 : -1; // Case 2 vs. 5,6. |
| } |
| // We are in case 2, 3, or 4. |
| return S2.orderedCCW(a0, b0, a2, ab1) ? 0 : -1; // Case 2,3 vs. 4. |
| } |
| } |
| |
| /** |
| * 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); |
| double acb = -(ab.dotProd(c)); |
| double bda = ab.dotProd(d); |
| if (acb * bda <= 0) { |
| return false; |
| } |
| |
| S2Point cd = S2Point.crossProd(c, d); |
| double cbd = -(cd.dotProd(b)); |
| double dac = cd.dotProd(a); |
| return (acb * cbd > 0) && (acb * dac > 0); |
| } |
| |
| /** |
| * Like SimpleCrossing, except that points that lie exactly on a line are |
| * arbitrarily classified as being on one side or the other (according to the |
| * rules of S2.robustCCW). It returns +1 if there is a crossing, -1 if there |
| * is no crossing, and 0 if any two vertices from different edges are the |
| * same. Returns 0 or -1 if either edge is degenerate. Properties of |
| * robustCrossing: |
| * |
| * (1) robustCrossing(b,a,c,d) == robustCrossing(a,b,c,d) (2) |
| * robustCrossing(c,d,a,b) == robustCrossing(a,b,c,d) (3) |
| * robustCrossing(a,b,c,d) == 0 if a==c, a==d, b==c, b==d (3) |
| * robustCrossing(a,b,c,d) <= 0 if a==b or c==d |
| * |
| * Note that if you want to check an edge against a *chain* of other edges, |
| * it is much more efficient to use an EdgeCrosser (above). |
| */ |
| public static int robustCrossing(S2Point a, S2Point b, S2Point c, S2Point d) { |
| // For there to be a crossing, the triangles ACB, CBD, BDA, DAC must |
| // all have the same orientation (clockwise or counterclockwise). |
| // |
| // First we compute the orientation of ACB and BDA. We permute the |
| // arguments to robustCCW so that we can reuse the cross-product of A and B. |
| // Recall that when the arguments to robustCCW are permuted, the sign of the |
| // result changes according to the sign of the permutation. Thus ACB and |
| // ABC are oppositely oriented, while BDA and ABD are the same. |
| S2Point aCrossB = S2Point.crossProd(a, b); |
| int acb = -S2.robustCCW(a, b, c, aCrossB); |
| int bda = S2.robustCCW(a, b, d, aCrossB); |
| |
| // If any two vertices are the same, the result is degenerate. |
| if ((bda & acb) == 0) { |
| return 0; |
| } |
| |
| // If ABC and BDA have opposite orientations (the most common case), |
| // there is no crossing. |
| if (bda != acb) { |
| return -1; |
| } |
| |
| // Otherwise we compute the orientations of CBD and DAC, and check whether |
| // their orientations are compatible with the other two triangles. |
| S2Point cCrossD = S2Point.crossProd(c, d); |
| int cbd = -S2.robustCCW(c, d, b, cCrossD); |
| if (cbd != acb) { |
| return -1; |
| } |
| |
| int dac = S2.robustCCW(c, d, a, cCrossD); |
| return (dac == acb) ? 1 : -1; |
| } |
| |
| /** |
| * Given two edges AB and CD where at least two vertices are identical (i.e. |
| * robustCrossing(a,b,c,d) == 0), this function defines whether the two edges |
| * "cross" in a such a way that point-in-polygon containment tests can be |
| * implemented by counting the number of edge crossings. The basic rule is |
| * that a "crossing" occurs if AB is encountered after CD during a CCW sweep |
| * around the shared vertex starting from a fixed reference point. |
| * |
| * Note that according to this rule, if AB crosses CD then in general CD does |
| * not cross AB. However, this leads to the correct result when counting |
| * polygon edge crossings. For example, suppose that A,B,C are three |
| * consecutive vertices of a CCW polygon. If we now consider the edge |
| * crossings of a segment BP as P sweeps around B, the crossing number changes |
| * parity exactly when BP crosses BA or BC. |
| * |
| * Useful properties of VertexCrossing (VC): |
| * |
| * (1) VC(a,a,c,d) == VC(a,b,c,c) == false (2) VC(a,b,a,b) == VC(a,b,b,a) == |
| * true (3) VC(a,b,c,d) == VC(a,b,d,c) == VC(b,a,c,d) == VC(b,a,d,c) (3) If |
| * exactly one of a,b equals one of c,d, then exactly one of VC(a,b,c,d) and |
| * VC(c,d,a,b) is true |
| * |
| * It is an error to call this method with 4 distinct vertices. |
| */ |
| public static boolean vertexCrossing(S2Point a, S2Point b, S2Point c, S2Point d) { |
| // If A == B or C == D there is no intersection. We need to check this |
| // case first in case 3 or more input points are identical. |
| if (a.equals(b) || c.equals(d)) { |
| return false; |
| } |
| |
| // If any other pair of vertices is equal, there is a crossing if and only |
| // if orderedCCW() indicates that the edge AB is further CCW around the |
| // shared vertex than the edge CD. |
| if (a.equals(d)) { |
| return S2.orderedCCW(S2.ortho(a), c, b, a); |
| } |
| if (b.equals(c)) { |
| return S2.orderedCCW(S2.ortho(b), d, a, b); |
| } |
| if (a.equals(c)) { |
| return S2.orderedCCW(S2.ortho(a), d, b, a); |
| } |
| if (b.equals(d)) { |
| return S2.orderedCCW(S2.ortho(b), c, a, b); |
| } |
| |
| // assert (false); |
| return false; |
| } |
| |
| /** |
| * A convenience function that calls robustCrossing() to handle cases where |
| * all four vertices are distinct, and VertexCrossing() to handle cases where |
| * two or more vertices are the same. This defines a crossing function such |
| * that point-in-polygon containment tests can be implemented by simply |
| * counting edge crossings. |
| */ |
| public static boolean edgeOrVertexCrossing(S2Point a, S2Point b, S2Point c, S2Point d) { |
| int crossing = robustCrossing(a, b, c, d); |
| if (crossing < 0) { |
| return false; |
| } |
| if (crossing > 0) { |
| return true; |
| } |
| return vertexCrossing(a, b, c, d); |
| } |
| |
| static class CloserResult { |
| private double dmin2; |
| private S2Point vmin; |
| |
| public double getDmin2() { |
| return dmin2; |
| } |
| |
| public S2Point getVmin() { |
| return vmin; |
| } |
| |
| public CloserResult(double dmin2, S2Point vmin) { |
| this.dmin2 = dmin2; |
| this.vmin = vmin; |
| } |
| |
| public void replaceIfCloser(S2Point x, S2Point y) { |
| // If the squared distance from x to y is less than dmin2, then replace |
| // vmin by y and update dmin2 accordingly. |
| double d2 = S2Point.minus(x, y).norm2(); |
| if (d2 < dmin2 || (d2 == dmin2 && y.lessThan(vmin))) { |
| dmin2 = d2; |
| vmin = y; |
| } |
| } |
| } |
| |
| /* |
| * Given two edges AB and CD such that robustCrossing() is true, return their |
| * intersection point. Useful properties of getIntersection (GI): |
| * |
| * (1) GI(b,a,c,d) == GI(a,b,d,c) == GI(a,b,c,d) (2) GI(c,d,a,b) == |
| * GI(a,b,c,d) |
| * |
| * The returned intersection point X is guaranteed to be close to the edges AB |
| * and CD, but if the edges intersect at a very small angle then X may not be |
| * close to the true mathematical intersection point P. See the description of |
| * "DEFAULT_INTERSECTION_TOLERANCE" below for details. |
| */ |
| public static S2Point getIntersection(S2Point a0, S2Point a1, S2Point b0, S2Point b1) { |
| Preconditions.checkArgument(robustCrossing(a0, a1, b0, b1) > 0, |
| "Input edges a0a1 and b0b1 muct have a true robustCrossing."); |
| |
| // We use robustCrossProd() to get accurate results even when two endpoints |
| // are close together, or when the two line segments are nearly parallel. |
| S2Point aNorm = S2Point.normalize(S2.robustCrossProd(a0, a1)); |
| S2Point bNorm = S2Point.normalize(S2.robustCrossProd(b0, b1)); |
| S2Point x = S2Point.normalize(S2.robustCrossProd(aNorm, bNorm)); |
| |
| // Make sure the intersection point is on the correct side of the sphere. |
| // Since all vertices are unit length, and edges are less than 180 degrees, |
| // (a0 + a1) and (b0 + b1) both have positive dot product with the |
| // intersection point. We use the sum of all vertices to make sure that the |
| // result is unchanged when the edges are reversed or exchanged. |
| if (x.dotProd(S2Point.add(S2Point.add(a0, a1), S2Point.add(b0, b1))) < 0) { |
| x = S2Point.neg(x); |
| } |
| |
| // The calculation above is sufficient to ensure that "x" is within |
| // DEFAULT_INTERSECTION_TOLERANCE of the great circles through (a0,a1) and |
| // (b0,b1). |
| // However, if these two great circles are very close to parallel, it is |
| // possible that "x" does not lie between the endpoints of the given line |
| // segments. In other words, "x" might be on the great circle through |
| // (a0,a1) but outside the range covered by (a0,a1). In this case we do |
| // additional clipping to ensure that it does. |
| |
| if (S2.orderedCCW(a0, x, a1, aNorm) && S2.orderedCCW(b0, x, b1, bNorm)) { |
| return x; |
| } |
| |
| // Find the acceptable endpoint closest to x and return it. An endpoint is |
| // acceptable if it lies between the endpoints of the other line segment. |
| CloserResult r = new CloserResult(10, x); |
| if (S2.orderedCCW(b0, a0, b1, bNorm)) { |
| r.replaceIfCloser(x, a0); |
| } |
| if (S2.orderedCCW(b0, a1, b1, bNorm)) { |
| r.replaceIfCloser(x, a1); |
| } |
| if (S2.orderedCCW(a0, b0, a1, aNorm)) { |
| r.replaceIfCloser(x, b0); |
| } |
| if (S2.orderedCCW(a0, b1, a1, aNorm)) { |
| r.replaceIfCloser(x, b1); |
| } |
| return r.getVmin(); |
| } |
| |
| /** |
| * Given a point X and an edge AB, return the distance ratio AX / (AX + BX). |
| * If X happens to be on the line segment AB, this is the fraction "t" such |
| * that X == Interpolate(A, B, t). Requires that A and B are distinct. |
| */ |
| public static double getDistanceFraction(S2Point x, S2Point a0, S2Point a1) { |
| Preconditions.checkArgument(!a0.equals(a1)); |
| double d0 = x.angle(a0); |
| double d1 = x.angle(a1); |
| return d0 / (d0 + d1); |
| } |
| |
| /** |
| * Return the minimum distance from X to any point on the edge AB. The result |
| * is very accurate for small distances but may have some numerical error if |
| * the distance is large (approximately Pi/2 or greater). The case A == B is |
| * handled correctly. Note: x, a and b must be of unit length. Throws |
| * IllegalArgumentException if this is not the case. |
| */ |
| public static S1Angle getDistance(S2Point x, S2Point a, S2Point b) { |
| return getDistance(x, a, b, S2.robustCrossProd(a, b)); |
| } |
| |
| /** |
| * A slightly more efficient version of getDistance() where the cross product |
| * of the two endpoints has been precomputed. The cross product does not need |
| * to be normalized, but should be computed using S2.robustCrossProd() for the |
| * most accurate results. |
| */ |
| public static S1Angle getDistance(S2Point x, S2Point a, S2Point b, S2Point aCrossB) { |
| Preconditions.checkArgument(S2.isUnitLength(x)); |
| Preconditions.checkArgument(S2.isUnitLength(a)); |
| Preconditions.checkArgument(S2.isUnitLength(b)); |
| |
| // There are three cases. If X is located in the spherical wedge defined by |
| // A, B, and the axis A x B, then the closest point is on the segment AB. |
| // Otherwise the closest point is either A or B; the dividing line between |
| // these two cases is the great circle passing through (A x B) and the |
| // midpoint of AB. |
| |
| if (S2.simpleCCW(aCrossB, a, x) && S2.simpleCCW(x, b, aCrossB)) { |
| // The closest point to X lies on the segment AB. We compute the distance |
| // to the corresponding great circle. The result is accurate for small |
| // distances but not necessarily for large distances (approaching Pi/2). |
| |
| double sinDist = Math.abs(x.dotProd(aCrossB)) / aCrossB.norm(); |
| return S1Angle.radians(Math.asin(Math.min(1.0, sinDist))); |
| } |
| |
| // Otherwise, the closest point is either A or B. The cheapest method is |
| // just to compute the minimum of the two linear (as opposed to spherical) |
| // distances and convert the result to an angle. Again, this method is |
| // accurate for small but not large distances (approaching Pi). |
| |
| double linearDist2 = Math.min(S2Point.minus(x, a).norm2(), S2Point.minus(x, b).norm2()); |
| return S1Angle.radians(2 * Math.asin(Math.min(1.0, 0.5 * Math.sqrt(linearDist2)))); |
| } |
| |
| /** |
| * Returns the point on edge AB closest to X. x, a and b must be of unit |
| * length. Throws IllegalArgumentException if this is not the case. |
| * |
| */ |
| public static S2Point getClosestPoint(S2Point x, S2Point a, S2Point b) { |
| Preconditions.checkArgument(S2.isUnitLength(x)); |
| Preconditions.checkArgument(S2.isUnitLength(a)); |
| Preconditions.checkArgument(S2.isUnitLength(b)); |
| |
| S2Point crossProd = S2.robustCrossProd(a, b); |
| // Find the closest point to X along the great circle through AB. |
| S2Point p = S2Point.minus(x, S2Point.mul(crossProd, x.dotProd(crossProd) / crossProd.norm2())); |
| |
| // If p is on the edge AB, then it's the closest point. |
| if (S2.simpleCCW(crossProd, a, p) && S2.simpleCCW(p, b, crossProd)) { |
| return S2Point.normalize(p); |
| } |
| // Otherwise, the closest point is either A or B. |
| return S2Point.minus(x, a).norm2() <= S2Point.minus(x, b).norm2() ? a : b; |
| } |
| |
| /** Constructor is private so that this class is never instantiated. */ |
| private S2EdgeUtil() { |
| } |
| } |