com.google.common.geometry.S2EdgeUtil Maven / Gradle / Ivy
The newest version!
/*
* 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() {
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy