cells and edges arrays. */
+func (e *EdgeIndex) sortIndex() {
+ // create an array of indices and sort based on the values in the parallel
+ // arrays at each index
+ indices := make([]int, len(e.cells))
+ for i := 0; i < len(indices); i++ {
+ indices[i] = i
+ }
+ sort.Sort(ByCellThenEdge{e, indices})
+ // copy the cells and edges in the order given by the sorted list of indices
+ newCells := make([]uint64, len(e.cells))
+ newEdges := make([]int, len(e.edges))
+ for i := 0; i < len(indices); i++ {
+ newCells[i] = e.cells[indices[i]]
+ newEdges[i] = e.edges[indices[i]]
+ }
+ // replace the cells and edges with the sorted arrays
+ e.cells = newCells
+ e.edges = newEdges
+}
+
+func (e *EdgeIndex) IsIndexComputed() bool {
+ return e.indexComputed
+}
+
+/**
+ * Tell the index that we just received a new request for candidates. Useful
+ * to compute when to switch to quad tree.
+ */
+func (e *EdgeIndex) incrementQueryCount() {
+ e.queryCount++
+}
+
+/**
+ * If the index hasn't been computed yet, looks at how much work has gone into
+ * iterating using the brute force method, and how much more work is planned
+ * as defined by 'cost'. If it were to have been cheaper to use a quad tree
+ * from the beginning, then compute it now. This guarantees that we will never
+ * use more than twice the time we would have used had we known in advance
+ * exactly how many edges we would have wanted to test. It is the theoretical
+ * best.
+ *
+ * The value 'n' is the number of iterators we expect to request from this
+ * edge index.
+ *
+ * If we have m data edges and n query edges, then the brute force cost is m
+ * * n * testCost where testCost is taken to be the cost of
+ * EdgeCrosser.robustCrossing, measured to be about 30ns at the time of this
+ * writing.
+ *
+ * If we compute the index, the cost becomes: m * costInsert + n *
+ * costFind(m)
+ *
+ * - costInsert can be expected to be reasonably stable, and was measured at
+ * 1200ns with the BM_QuadEdgeInsertionCost benchmark.
+ *
+ * - costFind depends on the length of the edge . For m=1000 edges, we got
+ * timings ranging from 1ms (edge the length of the polygon) to 40ms. The
+ * latter is for very long query edges, and needs to be optimized. We will
+ * assume for the rest of the discussion that costFind is roughly 3ms.
+ *
+ * When doing one additional query, the differential cost is m * testCost -
+ * costFind(m) With the numbers above, it is better to use the quad tree (if
+ * we have it) if m >= 100.
+ *
+ * If m = 100, 30 queries will give m*n*testCost = m * costInsert = 100ms,
+ * while the marginal cost to find is 3ms. Thus, this is a reasonable thing to
+ * do.
+ */
+func (e *EdgeIndex) PredictAdditionalCalls(n int) {
+ if e.indexComputed {
+ return
+ }
+ if e.getNumEdges() > 100 && (e.queryCount+n) > 30 {
+ e.ComputeIndex()
+ }
+}
+
+/**
+ * Appends to "candidateCrossings" all edge references which may cross the
+ * given edge. This is done by covering the edge and then finding all
+ * references of edges whose coverings overlap this covering. Parent cells are
+ * checked level by level. Child cells are checked all at once by taking
+ * advantage of the natural ordering of S2CellIds.
+ */
+func (e *EdgeIndex) findCandidateCrossings(a, b Point, candidateCrossings *[]int) {
+ // Preconditions.checkState(indexComputed);
+ cover := []CellID{}
+ e.getCovering(a, b, false, &cover)
+
+ // Edge references are inserted into the map once for each covering cell, so
+ // absorb duplicates here
+
+ uniqueSet := make(map[int]bool)
+ e.getEdgesInParentCells(cover, &uniqueSet)
+
+ // TODO(user): An important optimization for long query
+ // edges (Contains queries): keep a bounding cap and clip the query
+ // edge to the cap before starting the descent.
+ e.getEdgesInChildrenCells(a, b, &cover, &uniqueSet)
+
+ *candidateCrossings = make([]int, len(uniqueSet))
+ for k := range uniqueSet {
+ *candidateCrossings = append(*candidateCrossings, k)
+ }
+}
+
+/**
+ * Returns the smallest cell containing all four points, or
+ * {@link S2CellId#sentinel()} if they are not all on the same face. The
+ * points don't need to be normalized.
+ */
+func (e *EdgeIndex) containingCell4(pa, pb, pc, pd Point) CellID {
+ a := CellIDFromPoint(pa)
+ b := CellIDFromPoint(pb)
+ c := CellIDFromPoint(pc)
+ d := CellIDFromPoint(pd)
+
+ if a.Face() != b.Face() || a.Face() != c.Face() || a.Face() != d.Face() {
+ return CellIDSentinel()
+ }
+
+ for a != b || a != c || a != d {
+ a = a.immediateParent()
+ b = b.immediateParent()
+ c = c.immediateParent()
+ d = d.immediateParent()
+ }
+ return a
+}
+
+/**
+ * Returns the smallest cell containing both points, or Sentinel if they are
+ * not all on the same face. The points don't need to be normalized.
+ */
+func (e *EdgeIndex) containingCell2(pa, pb Point) CellID {
+ a := CellIDFromPoint(pa)
+ b := CellIDFromPoint(pb)
+
+ if a.Face() != b.Face() {
+ return CellIDSentinel()
+ }
+
+ for a != b {
+ a = a.immediateParent()
+ b = b.immediateParent()
+ }
+ return a
+}
+
+/**
+ * Computes a cell covering of an edge. Clears edgeCovering and returns the
+ * level of the s2 cells used in the covering (only one level is ever used for
+ * each call).
+ *
+ * If thickenEdge is true, the edge is thickened and extended by 1% of its
+ * length.
+ *
+ * It is guaranteed that no child of a covering cell will fully contain the
+ * covered edge.
+ */
+func (e *EdgeIndex) getCovering(a, b Point, thickenEdge bool, edgeCovering *[]CellID) int {
+ *edgeCovering = []CellID{}
+
+ // Selects the ideal s2 level at which to cover the edge, this will be the
+ // level whose S2 cells have a width roughly commensurate to the length of
+ // the edge. We multiply the edge length by 2*THICKENING to guarantee the
+ // thickening is honored (it's not a big deal if we honor it when we don't
+ // request it) when doing the covering-by-cap trick.
+ edgeLength := a.Angle(b.Vector).Radians()
+ idealLevel := S2_PROJECTION.MIN_WIDTH().getMaxLevel(edgeLength * (1 + 2*THICKENING))
+
+ var containingCellId CellID
+ if !thickenEdge {
+ containingCellId = e.containingCell2(a, b)
+ } else {
+ if idealLevel == MAX_LEVEL {
+ // If the edge is tiny, instabilities are more likely, so we
+ // want to limit the number of operations.
+ // We pretend we are in a cell much larger so as to trigger the
+ // 'needs covering' case, so we won't try to thicken the edge.
+ containingCellId = CellID(0xFFF0).Parent(3)
+ } else {
+ pq := b.Sub(a.Vector).Mul(THICKENING)
+ ortho := pq.Cross(a.Vector).Normalize().Mul(edgeLength * THICKENING)
+ p := a.Sub(pq)
+ q := b.Add(pq)
+ // If p and q were antipodal, the edge wouldn't be lengthened,
+ // and it could even flip! This is not a problem because
+ // idealLevel != 0 here. The farther p and q can be is roughly
+ // a quarter Earth away from each other, so we remain
+ // Theta(THICKENING).
+ containingCellId = e.containingCell4(
+ Point{p.Sub(ortho)},
+ Point{p.Add(ortho)},
+ Point{q.Sub(ortho)},
+ Point{q.Add(ortho)},
+ )
+ }
+ }
+
+ // Best case: edge is fully contained in a cell that's not too big.
+ if containingCellId != CellIDSentinel() && containingCellId.Level() >= idealLevel-2 {
+ *edgeCovering = append(*edgeCovering, containingCellId)
+ return containingCellId.Level()
+ }
+
+ if idealLevel == 0 {
+ // Edge is very long, maybe even longer than a face width, so the
+ // trick below doesn't work. For now, we will add the whole S2 sphere.
+ // TODO(user): Do something a tad smarter (and beware of the
+ // antipodal case).
+ for cellid := CellIDBegin(0); cellid != CellIDEnd(0); cellid = cellid.Next() {
+ *edgeCovering = append(*edgeCovering, cellid)
+ }
+ return 0
+ }
+ // TODO(user): Check trick below works even when vertex is at
+ // interface
+ // between three faces.
+
+ // Use trick as in S2PolygonBuilder.PointIndex.findNearbyPoint:
+ // Cover the edge by a cap centered at the edge midpoint, then cover
+ // the cap by four big-enough cells around the cell vertex closest to the
+ // cap center.
+ middle := Point{a.Add(b.Vector).Div(2).Normalize()}
+ actualLevel := min(idealLevel, MAX_LEVEL-1)
+ *edgeCovering = CellIDFromPoint(middle).VertexNeighbors(actualLevel)
+ return actualLevel
+}
+
+/**
+ * Filters a list of entries down to the inclusive range defined by the given
+ * cells, in O(log N) time.
+ *
+ * @param cell1 One side of the inclusive query range.
+ * @param cell2 The other side of the inclusive query range.
+ * @return An array of length 2, containing the start/end indices.
+ */
+func (e *EdgeIndex) getEdges(cell1, cell2 uint64) []int {
+ // ensure cell1 <= cell2
+ if cell1 > cell2 {
+ cell1, cell2 = cell2, cell1
+ }
+ // The binary search returns -N-1 to indicate an insertion point at index N,
+ // if an exact match cannot be found. Since the edge indices queried for are
+ // not valid edge indices, we will always get -N-1, so we immediately
+ // convert to N.
+ return []int{
+ -1 - e.binarySearch(cell1, math.MinInt32),
+ -1 - e.binarySearch(cell2, math.MaxInt32),
+ }
+}
+
+func (e *EdgeIndex) binarySearch(cell uint64, edge int) int {
+ low := 0
+ high := len(e.cells) - 1
+ for low <= high {
+ mid := (low + high) >> 1
+ cmp := compare(e.cells[mid], e.edges[mid], cell, edge)
+ if cmp < 0 {
+ low = mid + 1
+ } else if cmp > 0 {
+ high = mid - 1
+ } else {
+ return mid
+ }
+ }
+ return -(low + 1)
+}
+
+/**
+ * Adds to candidateCrossings all the edges present in any ancestor of any
+ * cell of cover, down to minimumS2LevelUsed. The cell->edge map is in the
+ * variable mapping.
+ */
+func (e *EdgeIndex) getEdgesInParentCells(cover []CellID, candidateCrossings *map[int]bool) {
+ // Find all parent cells of covering cells.
+ parentCells := make(map[CellID]bool)
+ for _, coverCell := range cover {
+ for parentLevel := coverCell.Level() - 1; parentLevel >= e.minimumS2LevelUsed; parentLevel-- {
+ if _, ok := parentCells[coverCell.Parent(parentLevel)]; ok {
+ break // cell is already in => parents are too.
+ }
+ parentCells[coverCell.Parent(parentLevel)] = true
+ }
+ }
+
+ // Put parent cell edge references into result.
+ for parentCell := range parentCells {
+ bounds := e.getEdges(uint64(parentCell), uint64(parentCell))
+ for i := bounds[0]; i < bounds[1]; i++ {
+ (*candidateCrossings)[e.edges[i]] = true
+ }
+ }
+}
+
+/**
+ * Returns true if ab possibly crosses cd, by clipping tiny angles to zero.
+ */
+func lenientCrossing(a, b, c, d Point) bool {
+ // assert (S2.isUnitLength(a));
+ // assert (S2.isUnitLength(b));
+ // assert (S2.isUnitLength(c));
+
+ acb := a.Cross(c.Vector).Dot(b.Vector)
+ bda := b.Cross(d.Vector).Dot(a.Vector)
+ if math.Abs(acb) < MAX_DET_ERROR || math.Abs(bda) < MAX_DET_ERROR {
+ return true
+ }
+ if acb*bda < 0 {
+ return false
+ }
+ cbd := c.Cross(b.Vector).Dot(d.Vector)
+ dac := c.Cross(a.Vector).Dot(c.Vector)
+ if math.Abs(cbd) < MAX_DET_ERROR || math.Abs(dac) < MAX_DET_ERROR {
+ return true
+ }
+ return (acb*cbd >= 0) && (acb*dac >= 0)
+}
+
+/**
+ * Returns true if the edge and the cell (including boundary) intersect.
+ */
+func edgeIntersectsCellBoundary(a, b Point, cell Cell) bool {
+ vertices := make([]Point, 4)
+ for i := 0; i < 4; i++ {
+ vertices[i] = cell.Vertex(i)
+ }
+ for i := 0; i < 4; i++ {
+ fromPoint := vertices[i]
+ toPoint := vertices[(i+1)%4]
+ if lenientCrossing(a, b, fromPoint, toPoint) {
+ return true
+ }
+ }
+ return false
+}
+
+/**
+ * Appends to candidateCrossings the edges that are fully contained in an S2
+ * covering of edge. The covering of edge used is initially cover, but is
+ * refined to eliminate quickly subcells that contain many edges but do not
+ * intersect with edge.
+ */
+func (e *EdgeIndex) getEdgesInChildrenCells(a, b Point, cover *[]CellID, candidateCrossings *map[int]bool) {
+ // Put all edge references of (covering cells + descendant cells) into
+ // result.
+ // This relies on the natural ordering of S2CellIds.
+ for len(*cover) > 0 {
+ cell := (*cover)[len(*cover)-1]
+ *cover = (*cover)[0 : len(*cover)-1]
+
+ bounds := e.getEdges(uint64(cell.RangeMin()), uint64(cell.RangeMax()))
+ if bounds[1]-bounds[0] <= 16 {
+ for i := bounds[0]; i < bounds[1]; i++ {
+ (*candidateCrossings)[e.edges[i]] = true
+ }
+ } else {
+ // Add cells at this level
+ bounds = e.getEdges(uint64(cell), uint64(cell))
+ for i := bounds[0]; i < bounds[1]; i++ {
+ (*candidateCrossings)[e.edges[i]] = true
+ }
+ // Recurse on the children -- hopefully some will be empty.
+ children := cell.Children()
+ for _, child := range children {
+ // TODO(user): Do the check for the four cells at once,
+ // as it is enough to check the four edges between the cells. At
+ // this time, we are checking 16 edges, 4 times too many.
+ //
+ // Note that given the guarantee of AppendCovering, it is enough
+ // to check that the edge intersect with the cell boundary as it
+ // cannot be fully contained in a cell.
+ if edgeIntersectsCellBoundary(a, b, CellFromCellID(child)) {
+ *cover = append(*cover, child)
+ }
+ }
+ }
+ }
+}
+
+/*
+ * An iterator on data edges that may cross a query edge (a,b). Create the
+ * iterator, call getCandidates(), then hasNext()/next() repeatedly.
+ *
+ * The current edge in the iteration has index index(), goes between from()
+ * and to().
+ */
+type DataEdgeIterator struct {
+ /**
+ * The structure containing the data edges.
+ */
+ edgeIndex *EdgeIndex
+
+ /**
+ * Tells whether getCandidates() obtained the candidates through brute force
+ * iteration or using the quad tree structure.
+ */
+ isBruteForce bool
+
+ /**
+ * Index of the current edge and of the edge before the last next() call.
+ */
+ currentIndex int
+
+ /**
+ * Cache of edgeIndex.getNumEdges() so that hasNext() doesn't make an extra
+ * call
+ */
+ numEdges int
+
+ /**
+ * All the candidates obtained by getCandidates() when we are using a
+ * quad-tree (i.e. isBruteForce = false).
+ */
+ candidates []int
+
+ /**
+ * Index within array above. We have: currentIndex =
+ * candidates.get(currentIndexInCandidates).
+ */
+ currentIndexInCandidates int
+}
+
+func NewDataEdgeIterator(edgeIndex *EdgeIndex) *DataEdgeIterator {
+ return &DataEdgeIterator{
+ edgeIndex: edgeIndex,
+ candidates: []int{},
+ }
+}
+
+/**
+ * Initializes the iterator to iterate over a set of candidates that may
+ * cross the edge (a,b).
+ */
+func (d *DataEdgeIterator) GetCandidates(a, b Point) {
+ d.edgeIndex.PredictAdditionalCalls(1)
+ d.isBruteForce = !d.edgeIndex.IsIndexComputed()
+ if d.isBruteForce {
+ d.edgeIndex.incrementQueryCount()
+ d.currentIndex = 0
+ d.numEdges = d.edgeIndex.getNumEdges()
+ } else {
+ d.candidates = []int{}
+ d.edgeIndex.findCandidateCrossings(a, b, &d.candidates)
+ d.currentIndexInCandidates = 0
+ if len(d.candidates) > 0 {
+ d.currentIndex = d.candidates[0]
+ }
+ }
+}
+
+/**
+ * Index of the current edge in the iteration.
+ */
+func (d *DataEdgeIterator) Index() int {
+ if !d.HasNext() {
+ panic("No next candidate, use HasNext")
+ }
+ return d.currentIndex
+}
+
+/**
+ * False if there are no more candidates; true otherwise.
+ */
+func (d *DataEdgeIterator) HasNext() bool {
+ if d.isBruteForce {
+ return d.currentIndex < d.numEdges
+ } else {
+ return d.currentIndexInCandidates < len(d.candidates)
+ }
+}
+
+/**
+ * Iterate to the next available candidate.
+ */
+func (d *DataEdgeIterator) Next() {
+ if !d.HasNext() {
+ panic("No next candidate, use HasNext")
+ }
+ if d.isBruteForce {
+ d.currentIndex++
+ } else {
+ d.currentIndexInCandidates++
+ if d.currentIndexInCandidates < len(d.candidates) {
+ d.currentIndex = d.candidates[d.currentIndexInCandidates]
+ }
+ }
+}
diff --git a/s2/edgeutil.go b/s2/edgeutil.go
new file mode 100644
index 00000000..5cf19609
--- /dev/null
+++ b/s2/edgeutil.go
@@ -0,0 +1,417 @@
+package s2
+
+import (
+ "math"
+
+ "github.com/golang/geo/r1"
+ "github.com/golang/geo/s1"
+)
+
+type EdgeCrosser struct {
+ // The fields below are all constant.
+ a Point
+ b Point
+ aCrossB Point
+
+ // The fields below are updated for each vertex in the chain.
+
+ // Previous vertex in the vertex chain.
+ c Point
+ // The orientation of the triangle ACB.
+ acb int
+}
+
+func NewEdgeCrosser(a, b, c Point) *EdgeCrosser {
+ ec := &EdgeCrosser{
+ a: a,
+ b: b,
+ aCrossB: Point{a.Cross(b.Vector)},
+ }
+ ec.RestartAt(c)
+ return ec
+}
+
+func (ec *EdgeCrosser) RestartAt(c Point) {
+ ec.c = c
+ ec.acb = -int(RobustCCWWithCross(ec.a, ec.b, ec.c, ec.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.
+ */
+func (ec *EdgeCrosser) RobustCrossing(d Point) int {
+ // 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.
+ bda := int(RobustCCWWithCross(ec.a, ec.b, d, ec.aCrossB))
+ var result int
+
+ if bda == -ec.acb && bda != 0 {
+ // Most common case -- triangles have opposite orientations.
+ result = -1
+ } else if (bda & ec.acb) == 0 {
+ // At least one value is zero -- two vertices are identical.
+ result = 0
+ } else {
+ // assert (bda == acb && bda != 0);
+ result = ec.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).
+ ec.c = d
+ ec.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.
+ */
+func (ec *EdgeCrosser) EdgeOrVertexCrossing(d Point) bool {
+ // We need to copy c since it is clobbered by robustCrossing().
+ c2 := PointFromCoordsRaw(ec.c.X, ec.c.Y, ec.c.Z)
+
+ crossing := ec.RobustCrossing(d)
+ if crossing < 0 {
+ return false
+ }
+ if crossing > 0 {
+ return true
+ }
+
+ return VertexCrossing(ec.a, ec.b, c2, d)
+}
+
+/**
+ * This function handles the "slow path" of robustCrossing().
+ */
+func (ec *EdgeCrosser) robustCrossingInternal(d Point) int {
+ // ACB and BDA have the appropriate orientations, so now we check the
+ // triangles CBD and DAC.
+ cCrossD := Point{ec.c.Cross(d.Vector)}
+ cbd := -int(RobustCCWWithCross(ec.c, d, ec.b, cCrossD))
+ if cbd != ec.acb {
+ return -1
+ }
+
+ dac := int(RobustCCWWithCross(ec.c, d, ec.a, cCrossD))
+ if dac == ec.acb {
+ return 1
+ } else {
+ return -1
+ }
+}
+
+type RectBounder struct {
+ a Point
+ aLatLng LatLng
+ bound Rect
+}
+
+func NewRectBounder() *RectBounder {
+ return &RectBounder{bound: EmptyRect()}
+}
+
+func (rb *RectBounder) AddPoint(b Point) {
+ bLatLng := LatLngFromPoint(b)
+
+ if rb.bound.IsEmpty() {
+ rb.bound = rb.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.
+ rb.bound = rb.bound.Union(RectFromLatLngPointPair(rb.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.
+ aCrossB := rb.a.PointCross(b)
+ dir := aCrossB.Cross(PointFromCoordsRaw(0, 0, 1).Vector)
+ da := dir.Dot(rb.a.Vector)
+ db := dir.Dot(b.Vector)
+
+ if da*db < 0 {
+ // Minimum/maximum latitude occurs in the edge interior. This affects
+ // the latitude bounds but not the longitude bounds.
+ absLat := math.Acos(math.Abs(aCrossB.Z / aCrossB.Norm()))
+ lat := rb.bound.Lat
+ if da < 0 {
+ // It's possible that absLat < lat.lo() due to numerical errors.
+ lat = r1.IntervalFromPointPair(lat.Lo, math.Max(absLat, rb.bound.Lat.Hi))
+ } else {
+ lat = r1.IntervalFromPointPair(math.Min(-absLat, rb.bound.Lat.Lo), lat.Hi)
+ }
+ rb.bound = Rect{lat, rb.bound.Lng}
+ }
+ }
+
+ rb.a = b
+ rb.aLatLng = bLatLng
+}
+
+func (rb *RectBounder) GetBound() Rect { return rb.bound }
+
+/**
+ * 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.
+ */
+func VertexCrossing(a, b, c, d Point) bool {
+ // 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 OrderedCCW(Point{a.Ortho()}, c, b, a)
+ }
+ if b.Equals(c) {
+ return OrderedCCW(Point{b.Ortho()}, d, a, b)
+ }
+ if a.Equals(c) {
+ return OrderedCCW(Point{a.Ortho()}, d, b, a)
+ }
+ if b.Equals(d) {
+ return OrderedCCW(Point{b.Ortho()}, c, a, b)
+ }
+
+ // assert (false);
+ return false
+}
+
+/**
+ * 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.
+ */
+type WedgeRelation interface {
+ Test(a0, ab1, a2, b0, b2 Point) int
+}
+
+/**
+ * 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.
+ */
+type WedgeContains struct{}
+
+func (w WedgeContains) Test(a0, ab1, a2, b0, b2 Point) int {
+ // 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.
+ if OrderedCCW(a2, b2, b0, ab1) && OrderedCCW(b0, a0, a2, ab1) {
+ return 1
+ }
+ return 0
+}
+
+/**
+ * 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.
+ */
+type WedgeIntersects struct{}
+
+func (w WedgeIntersects) Test(a0, ab1, a2, b0, b2 Point) int {
+ // 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.
+ if OrderedCCW(a0, b2, b0, ab1) && OrderedCCW(b0, a2, a0, ab1) {
+ return 0
+ }
+ return -1
+}
+
+/**
+ * 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.
+ */
+type WedgeContainsOrIntersects struct{}
+
+func (w WedgeContainsOrIntersects) Test(a0, ab1, a2, b0, b2 Point) int {
+ // 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 OrderedCCW(a0, a2, b2, ab1) {
+ // We are in case 1, 5, or 6, or case 2 if a2 == b2.
+ if OrderedCCW(b2, b0, a0, ab1) {
+ return 1 // Case 1
+ }
+ return -1 // Case 2,5,6.
+ }
+ // We are in cases 2, 3, or 4.
+ if !OrderedCCW(a2, b0, b2, ab1) {
+ return 0 // Case 3.
+ }
+
+ // We are in case 2 or 4, or case 3 if a2 == b0.
+ if a2.Equals(b0) {
+ return 0 // Case 3
+ }
+ return -1 // Case 2,4.
+}
+
+/**
+ * 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.
+ */
+type WedgeContainsOrCrosses struct{}
+
+func (w WedgeContainsOrCrosses) Test(a0, ab1, a2, b0, b2 Point) int {
+ // 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 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 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.
+ if a2.Equals(b2) {
+ return 0 // Case 2
+ }
+ return -1 // Case 5,6.
+ }
+ // We are in case 2, 3, or 4.
+ if OrderedCCW(a0, b0, a2, ab1) {
+ return 0 // Case 2,3
+ }
+ return -1 // Case 4.
+}
+
+/**
+ * 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.
+ */
+func getDistanceFraction(x, a0, a1 Point) float64 {
+ if a0.Equals(a1) {
+ panic("a0 and a1 are equal")
+ }
+ d0 := x.Distance(a0).Radians()
+ d1 := x.Distance(a1).Radians()
+ 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.
+ */
+func getDistance(x, a, b Point) s1.Angle {
+ return getDistanceWithCross(x, a, b, a.PointCross(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.
+ */
+func getDistanceWithCross(x, a, b, aCrossB Point) s1.Angle {
+ if !x.IsUnit() || !a.IsUnit() || !b.IsUnit() {
+ panic("x, a and b need to be unit length")
+ }
+
+ // 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 simpleCCW(aCrossB, a, x) && 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).
+
+ sinDist := math.Abs(x.Dot(aCrossB.Vector)) / aCrossB.Norm()
+ return s1.Angle(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).
+
+ linearDist2 := math.Min(x.Sub(a.Vector).Norm2(), x.Sub(b.Vector).Norm2())
+ return s1.Angle(2 * math.Asin(math.Min(1.0, 0.5*math.Sqrt(linearDist2))))
+}
diff --git a/s2/latlng.go b/s2/latlng.go
index 6111b9b4..b34efae8 100644
--- a/s2/latlng.go
+++ b/s2/latlng.go
@@ -40,6 +40,10 @@ func (ll LatLng) IsValid() bool {
func (ll LatLng) String() string { return fmt.Sprintf("[%v, %v]", ll.Lat, ll.Lng) }
+func (ll LatLng) StringDegrees() string {
+ return fmt.Sprintf("[%f, %f]", ll.Lat.Degrees(), ll.Lng.Degrees())
+}
+
// Distance returns the angle between two LatLngs.
func (ll LatLng) Distance(ll2 LatLng) s1.Angle {
// Haversine formula, as used in C++ S2LatLng::GetDistance.
@@ -48,26 +52,18 @@ func (ll LatLng) Distance(ll2 LatLng) s1.Angle {
dlat := math.Sin(0.5 * (lat2 - lat1))
dlng := math.Sin(0.5 * (lng2 - lng1))
x := dlat*dlat + dlng*dlng*math.Cos(lat1)*math.Cos(lat2)
- return s1.Angle(2*math.Atan2(math.Sqrt(x), math.Sqrt(math.Max(0, 1-x)))) * s1.Radian
+ return s1.Angle(2 * math.Atan2(math.Sqrt(x), math.Sqrt(math.Max(0, 1-x))))
}
// NOTE(mikeperrow): The C++ implementation publicly exposes latitude/longitude
// functions. Let's see if that's really necessary before exposing the same functionality.
func latitude(p Point) s1.Angle {
- return s1.Angle(math.Atan2(p.Z, math.Sqrt(p.X*p.X+p.Y*p.Y))) * s1.Radian
+ return s1.Angle(math.Atan2(p.Z, math.Sqrt(p.X*p.X+p.Y*p.Y)))
}
func longitude(p Point) s1.Angle {
- return s1.Angle(math.Atan2(p.Y, p.X)) * s1.Radian
-}
-
-// PointFromLatLng returns an Point for the given LatLng.
-func PointFromLatLng(ll LatLng) Point {
- phi := ll.Lat.Radians()
- theta := ll.Lng.Radians()
- cosphi := math.Cos(phi)
- return PointFromCoords(math.Cos(theta)*cosphi, math.Sin(theta)*cosphi, math.Sin(phi))
+ return s1.Angle(math.Atan2(p.Y, p.X))
}
// LatLngFromPoint returns an LatLng for a given Point.
diff --git a/s2/loop.go b/s2/loop.go
new file mode 100644
index 00000000..12ee4098
--- /dev/null
+++ b/s2/loop.go
@@ -0,0 +1,824 @@
+package s2
+
+import (
+ "fmt"
+ "math"
+
+ "github.com/golang/geo/r1"
+ "github.com/golang/geo/s1"
+)
+
+/**
+ *
+ * An S2Loop represents a simple spherical polygon. It consists of a single
+ * chain of vertices where the first vertex is implicitly connected to the last.
+ * All loops are defined to have a CCW orientation, i.e. the interior of the
+ * polygon is on the left side of the edges. This implies that a clockwise loop
+ * enclosing a small area is interpreted to be a CCW loop enclosing a very large
+ * area.
+ *
+ * Loops are not allowed to have any duplicate vertices (whether adjacent or
+ * not), and non-adjacent edges are not allowed to intersect. Loops must have at
+ * least 3 vertices. Although these restrictions are not enforced in optimized
+ * code, you may get unexpected results if they are violated.
+ *
+ * Point containment is defined such that if the sphere is subdivided into
+ * faces (loops), every point is contained by exactly one face. This implies
+ * that loops do not necessarily contain all (or any) of their vertices An
+ * S2LatLngRect represents a latitude-longitude rectangle. It is capable of
+ * representing the empty and full rectangles as well as single points.
+ *
+ */
+type Loop struct {
+ // Edge index used for performance-critical operations. For example,
+ // contains() can determine whether a point is inside a loop in nearly
+ // constant time, whereas without an edge index it is forced to compare the
+ // query point against every edge in the loop.
+ index *EdgeIndex
+
+ // Maps each S2Point to its order in the loop, from 1 to numVertices.
+ vertexToIndex map[Point]int
+
+ vertices []Point
+
+ // The index (into "vertices") of the vertex that comes first in the total
+ // ordering of all vertices in this loop.
+ firstLogicalVertex int
+
+ bound Rect
+ originInside bool
+ depth int
+}
+
+func LoopFromPoints(points []Point) *Loop {
+ l := &Loop{
+ vertices: points,
+ bound: FullRect(),
+ depth: 0,
+ }
+ l.initOrigin()
+ l.initBound()
+ l.initFirstLogicalVertex()
+ return l
+}
+
+func LoopFromCell(cell Cell) *Loop {
+ return LoopFromCellAndRect(cell, cell.RectBound())
+}
+
+func LoopFromCellAndRect(cell Cell, bound Rect) *Loop {
+ l := &Loop{
+ bound: bound,
+ vertices: make([]Point, 4),
+ }
+ for i := 0; i < 4; i++ {
+ l.vertices[i] = cell.Vertex(i)
+ }
+ l.initOrigin()
+ l.initFirstLogicalVertex()
+ return l
+}
+
+func (l *Loop) Depth() int { return l.depth }
+
+/**
+ * The depth of a loop is defined as its nesting level within its containing
+ * polygon. "Outer shell" loops have depth 0, holes within those loops have
+ * depth 1, shells within those holes have depth 2, etc. This field is only
+ * used by the S2Polygon implementation.
+ *
+ * @param depth
+ */
+func (l *Loop) SetDepth(depth int) { l.depth = depth }
+
+/**
+ * Return true if this loop represents a hole in its containing polygon.
+ */
+func (l *Loop) IsHole() bool {
+ return (l.depth & 1) != 0
+}
+
+/**
+ * The sign of a loop is -1 if the loop represents a hole in its containing
+ * polygon, and +1 otherwise.
+ */
+func (l *Loop) Sign() int {
+ if l.IsHole() {
+ return -1
+ } else {
+ return 1
+ }
+}
+
+func (l *Loop) NumVertices() int {
+ return len(l.vertices)
+}
+
+/**
+ * For convenience, we make two entire copies of the vertex list available:
+ * vertex(n..2*n-1) is mapped to vertex(0..n-1), where n == numVertices().
+ */
+func (l *Loop) Vertex(i int) Point {
+ if i >= l.NumVertices() {
+ i = i - l.NumVertices()
+ }
+ return l.vertices[i]
+}
+
+func (l *Loop) CompareTo(other *Loop) int {
+ if l.NumVertices() != other.NumVertices() {
+ return l.NumVertices() - other.NumVertices()
+ }
+ // Compare the two loops' vertices, starting with each loop's
+ // firstLogicalVertex. This allows us to always catch cases where logically
+ // identical loops have different vertex orderings (e.g. ABCD and BCDA).
+ maxVertices := l.NumVertices()
+ iThis := l.firstLogicalVertex
+ iOther := other.firstLogicalVertex
+ for i := 0; i < maxVertices; i++ {
+ compare := l.Vertex(iThis).CompareTo(other.Vertex(iOther))
+ if compare != 0 {
+ return compare
+ }
+ iThis++
+ iOther++
+ }
+ return 0
+}
+
+/**
+ * Calculates firstLogicalVertex, the vertex in this loop that comes first in
+ * a total ordering of all vertices (by way of S2Point's compareTo function).
+ */
+func (l *Loop) initFirstLogicalVertex() {
+ first := 0
+ for i := 1; i < l.NumVertices(); i++ {
+ if l.Vertex(i).CompareTo(l.Vertex(first)) < 0 {
+ first = i
+ }
+ }
+ l.firstLogicalVertex = first
+}
+
+/**
+ * Return true if the loop area is at most 2*Pi.
+ */
+func (l *Loop) IsNormalized() bool {
+ // We allow a bit of error so that exact hemispheres are
+ // considered normalized.
+ return l.GetArea() <= 2*math.Pi+1e-14
+}
+
+/**
+ * Invert the loop if necessary so that the area enclosed by the loop is at
+ * most 2*Pi.
+ */
+func (l *Loop) Normalize() {
+ if !l.IsNormalized() {
+ l.Invert()
+ }
+}
+
+/**
+ * Reverse the order of the loop vertices, effectively complementing the
+ * region represented by the loop.
+ */
+func (l *Loop) Invert() {
+ last := l.NumVertices() - 1
+ for i := (last - 1) / 2; i >= 0; i-- {
+ t := l.vertices[i]
+ l.vertices[i] = l.vertices[last-i]
+ l.vertices[last-i] = t
+ }
+ l.vertexToIndex = nil
+ l.index = nil
+ l.originInside = !l.originInside
+ if l.bound.Lat.Lo > -math.Pi/2 && l.bound.Lat.Hi < math.Pi/2 {
+ // The complement of this loop contains both poles.
+ l.bound = FullRect()
+ } else {
+ l.initBound()
+ }
+ l.initFirstLogicalVertex()
+}
+
+/**
+ * Helper method to get area and optionally centroid.
+ */
+func (l *Loop) getAreaCentroid(doCentroid bool) AreaCentroid {
+ var centroid *Point = nil
+ // Don't crash even if loop is not well-defined.
+ if l.NumVertices() < 3 {
+ return NewAreaCentroid(0, centroid)
+ }
+
+ // The triangle area calculation becomes numerically unstable as the length
+ // of any edge approaches 180 degrees. However, a loop may contain vertices
+ // that are 180 degrees apart and still be valid, e.g. a loop that defines
+ // the northern hemisphere using four points. We handle this case by using
+ // triangles centered around an origin that is slightly displaced from the
+ // first vertex. The amount of displacement is enough to get plenty of
+ // accuracy for antipodal points, but small enough so that we still get
+ // accurate areas for very tiny triangles.
+ //
+ // Of course, if the loop contains a point that is exactly antipodal from
+ // our slightly displaced vertex, the area will still be unstable, but we
+ // expect this case to be very unlikely (i.e. a polygon with two vertices on
+ // opposite sides of the Earth with one of them displaced by about 2mm in
+ // exactly the right direction). Note that the approximate point resolution
+ // using the E7 or S2CellId representation is only about 1cm.
+
+ origin := l.Vertex(0)
+ axis := (origin.LargestAbsComponent() + 1) % 3
+ slightlyDisplaced := origin.GetAxis(axis) + math.E*1e-10
+ switch axis {
+ case 0:
+ origin = PointFromCoordsRaw(slightlyDisplaced, origin.Y, origin.Z)
+ case 1:
+ origin = PointFromCoordsRaw(origin.X, slightlyDisplaced, origin.Z)
+ case 2:
+ origin = PointFromCoordsRaw(origin.X, origin.Y, slightlyDisplaced)
+ }
+ origin = Point{origin.Normalize()}
+
+ var areaSum float64 = 0
+ centroidSum := PointFromCoordsRaw(0, 0, 0)
+ for i := 1; i <= l.NumVertices(); i++ {
+ areaSum += SignedArea(origin, l.Vertex(i-1), l.Vertex(i))
+ if doCentroid {
+ // The true centroid is already premultiplied by the triangle area.
+ trueCentroid := TrueCentroid(origin, l.Vertex(i-1), l.Vertex(i))
+ centroidSum = Point{centroidSum.Add(trueCentroid.Vector)}
+ }
+ }
+ // The calculated area at this point should be between -4*Pi and 4*Pi,
+ // although it may be slightly larger or smaller than this due to
+ // numerical errors.
+ // assert (Math.abs(areaSum) <= 4 * S2.M_PI + 1e-12);
+
+ if areaSum < 0 {
+ // If the area is negative, we have computed the area to the right of the
+ // loop. The area to the left is 4*Pi - (-area). Amazingly, the centroid
+ // does not need to be changed, since it is the negative of the integral
+ // of position over the region to the right of the loop. This is the same
+ // as the integral of position over the region to the left of the loop,
+ // since the integral of position over the entire sphere is (0, 0, 0).
+ areaSum += 4 * math.Pi
+ }
+ // The loop's sign() does not affect the return result and should be taken
+ // into account by the caller.
+ if doCentroid {
+ centroid = &Point{centroidSum.Vector}
+ }
+ return NewAreaCentroid(areaSum, centroid)
+}
+
+/**
+ * Return the area of the loop interior, i.e. the region on the left side of
+ * the loop. The return value is between 0 and 4*Pi and the true centroid of
+ * the loop multiplied by the area of the loop (see S2.java for details on
+ * centroids). Note that the centroid may not be contained by the loop.
+ */
+func (l *Loop) GetAreaAndCentroid() AreaCentroid {
+ return l.getAreaCentroid(true)
+}
+
+/**
+ * Return the area of the polygon interior, i.e. the region on the left side
+ * of an odd number of loops. The return value is between 0 and 4*Pi.
+ */
+func (l *Loop) GetArea() float64 {
+ return l.getAreaCentroid(false).GetArea()
+}
+
+/**
+ * Return the true centroid of the polygon multiplied by the area of the
+ * polygon (see {@link S2} for details on centroids). Note that the centroid
+ * may not be contained by the polygon.
+ */
+func (l *Loop) GetCentroid() Point {
+ return l.getAreaCentroid(true).GetCentroid()
+}
+
+func (l *Loop) ContainsLoop(b *Loop) bool {
+ // For this loop A to contains the given loop B, all of the following must
+ // be true:
+ //
+ // (1) There are no edge crossings between A and B except at vertices.
+ //
+ // (2) At every vertex that is shared between A and B, the local edge
+ // ordering implies that A contains B.
+ //
+ // (3) If there are no shared vertices, then A must contain a vertex of B
+ // and B must not contain a vertex of A. (An arbitrary vertex may be
+ // chosen in each case.)
+ //
+ // The second part of (3) is necessary to detect the case of two loops whose
+ // union is the entire sphere, i.e. two loops that contains each other's
+ // boundaries but not each other's interiors.
+
+ if !l.bound.ContainsRect(b.RectBound()) {
+ return false
+ }
+
+ // Unless there are shared vertices, we need to check whether A contains a
+ // vertex of B. Since shared vertices are rare, it is more efficient to do
+ // this test up front as a quick rejection test.
+ if !l.ContainsPoint(b.Vertex(0)) && l.findVertex(b.Vertex(0)) < 0 {
+ return false
+ }
+
+ // Now check whether there are any edge crossings, and also check the loop
+ // relationship at any shared vertices.
+ if l.checkEdgeCrossings(b, WedgeContains{}) <= 0 {
+ return false
+ }
+
+ // At this point we know that the boundaries of A and B do not intersect,
+ // and that A contains a vertex of B. However we still need to check for
+ // the case mentioned above, where (A union B) is the entire sphere.
+ // Normally this check is very cheap due to the bounding box precondition.
+ if l.bound.Union(b.RectBound()).IsFull() {
+ if b.ContainsPoint(l.Vertex(0)) && b.findVertex(l.Vertex(0)) < 0 {
+ return false
+ }
+ }
+ return true
+}
+
+/**
+ * Return true if the region contained by this loop intersects the region
+ * contained by the given other loop.
+ */
+func (l *Loop) IntersectsLoop(b *Loop) bool {
+ // a->Intersects(b) if and only if !a->Complement()->Contains(b).
+ // This code is similar to Contains(), but is optimized for the case
+ // where both loops enclose less than half of the sphere.
+
+ if !l.bound.IntersectsRect(b.RectBound()) {
+ return false
+ }
+
+ // Normalize the arguments so that B has a smaller longitude span than A.
+ // This makes intersection tests much more efficient in the case where
+ // longitude pruning is used (see CheckEdgeCrossings).
+ if b.RectBound().Lng.Length() > l.bound.Lng.Length() {
+ return b.IntersectsLoop(l)
+ }
+
+ // Unless there are shared vertices, we need to check whether A contains a
+ // vertex of B. Since shared vertices are rare, it is more efficient to do
+ // this test up front as a quick acceptance test.
+ if l.ContainsPoint(b.Vertex(0)) && l.findVertex(b.Vertex(0)) < 0 {
+ return true
+ }
+
+ // Now check whether there are any edge crossings, and also check the loop
+ // relationship at any shared vertices.
+ if l.checkEdgeCrossings(b, WedgeIntersects{}) < 0 {
+ return true
+ }
+
+ // We know that A does not contain a vertex of B, and that there are no edge
+ // crossings. Therefore the only way that A can intersect B is if B
+ // entirely contains A. We can check this by testing whether B contains an
+ // arbitrary non-shared vertex of A. Note that this check is cheap because
+ // of the bounding box precondition and the fact that we normalized the
+ // arguments so that A's longitude span is at least as long as B's.
+ if b.RectBound().ContainsRect(l.bound) {
+ if b.ContainsPoint(l.Vertex(0)) && b.findVertex(l.Vertex(0)) < 0 {
+ return true
+ }
+ }
+
+ return false
+}
+
+/**
+ * Given two loops of a polygon, return true if A contains B. This version of
+ * contains() is much cheaper since it does not need to check whether the
+ * boundaries of the two loops cross.
+ */
+func (l *Loop) ContainsNested(b *Loop) bool {
+ if !l.bound.ContainsRect(b.RectBound()) {
+ return false
+ }
+
+ // We are given that A and B do not share any edges, and that either one
+ // loop contains the other or they do not intersect.
+ m := l.findVertex(b.Vertex(1))
+ if m < 0 {
+ // Since b->vertex(1) is not shared, we can check whether A contains it.
+ return l.ContainsPoint(b.Vertex(1))
+ }
+ // Check whether the edge order around b->vertex(1) is compatible with
+ // A containin B.
+ return (WedgeContains{}).Test(l.Vertex(m-1), l.Vertex(m), l.Vertex(m+1), b.Vertex(0), b.Vertex(2)) > 0
+}
+
+/**
+ * Return +1 if A contains B (i.e. the interior of B is a subset of the
+ * interior of A), -1 if the boundaries of A and B cross, and 0 otherwise.
+ * Requires that A does not properly contain the complement of B, i.e. A and B
+ * do not contain each other's boundaries. This method is used for testing
+ * whether multi-loop polygons contain each other.
+ */
+func (l *Loop) ContainsOrCrosses(b *Loop) int {
+ // There can be containment or crossing only if the bounds intersect.
+ if !l.bound.IntersectsRect(b.RectBound()) {
+ return 0
+ }
+
+ // Now check whether there are any edge crossings, and also check the loop
+ // relationship at any shared vertices. Note that unlike Contains() or
+ // Intersects(), we can't do a point containment test as a shortcut because
+ // we need to detect whether there are any edge crossings.
+ result := l.checkEdgeCrossings(b, WedgeContainsOrCrosses{})
+
+ // If there was an edge crossing or a shared vertex, we know the result
+ // already. (This is true even if the result is 1, but since we don't
+ // bother keeping track of whether a shared vertex was seen, we handle this
+ // case below.)
+ if result <= 0 {
+ return result
+ }
+
+ // At this point we know that the boundaries do not intersect, and we are
+ // given that (A union B) is a proper subset of the sphere. Furthermore
+ // either A contains B, or there are no shared vertices (due to the check
+ // above). So now we just need to distinguish the case where A contains B
+ // from the case where B contains A or the two loops are disjoint.
+ if !l.bound.ContainsRect(b.RectBound()) {
+ return 0
+ }
+ if !l.ContainsPoint(b.Vertex(0)) && l.findVertex(b.Vertex(0)) < 0 {
+ return 0
+ }
+
+ return 1
+}
+
+/**
+ * Returns true if two loops have the same boundary except for vertex
+ * perturbations. More precisely, the vertices in the two loops must be in the
+ * same cyclic order, and corresponding vertex pairs must be separated by no
+ * more than maxError. Note: This method mostly useful only for testing
+ * purposes.
+ */
+func (l *Loop) BoundaryApproxEquals(b *Loop, maxError float64) bool {
+ if l.NumVertices() != b.NumVertices() {
+ return false
+ }
+ maxVertices := l.NumVertices()
+ iThis := l.firstLogicalVertex
+ iOther := b.firstLogicalVertex
+ for i := 0; i < maxVertices; i++ {
+ if !l.Vertex(iThis).ApproxEquals(b.Vertex(iOther), maxError) {
+ return false
+ }
+ iThis++
+ iOther++
+ }
+ return true
+}
+
+// CapBound returns a bounding spherical cap. This is not guaranteed to be exact.
+func (l *Loop) CapBound() Cap {
+ return l.bound.CapBound()
+}
+
+// RectBound returns a bounding latitude-longitude rectangle that contains
+// the region. The bounds are not guaranteed to be tight.
+func (l *Loop) RectBound() Rect {
+ return l.bound
+}
+
+/**
+ * If this method returns true, the region completely contains the given cell.
+ * Otherwise, either the region does not contain the cell or the containment
+ * relationship could not be determined.
+ */
+func (l *Loop) ContainsCell(cell Cell) bool {
+ // It is faster to construct a bounding rectangle for an S2Cell than for
+ // a general polygon. A future optimization could also take advantage of
+ // the fact than an S2Cell is convex.
+
+ cellBound := cell.RectBound()
+ if !l.bound.ContainsRect(cellBound) {
+ return false
+ }
+ cellLoop := LoopFromCellAndRect(cell, cellBound)
+ return l.ContainsLoop(cellLoop)
+}
+
+/**
+ * If this method returns false, the region does not intersect the given cell.
+ * Otherwise, either region intersects the cell, or the intersection
+ * relationship could not be determined.
+ */
+func (l *Loop) IntersectsCell(cell Cell) bool {
+ // It is faster to construct a bounding rectangle for an S2Cell than for
+ // a general polygon. A future optimization could also take advantage of
+ // the fact than an S2Cell is convex.
+
+ cellBound := cell.RectBound()
+ if !l.bound.IntersectsRect(cellBound) {
+ return false
+ }
+ return LoopFromCellAndRect(cell, cellBound).IntersectsLoop(l)
+}
+
+/**
+ * The point 'p' does not need to be normalized.
+ */
+func (l *Loop) ContainsPoint(p Point) bool {
+ if !l.bound.ContainsLatLng(LatLngFromPoint(p)) {
+ return false
+ }
+
+ inside := l.originInside
+ origin := PointFromCoordsRaw(0, 1, 0)
+ crosser := NewEdgeCrosser(origin, p, l.Vertex(l.NumVertices()-1))
+
+ // The s2edgeindex library is not optimized yet for long edges,
+ // so the tradeoff to using it comes with larger loops.
+ if l.NumVertices() < 2000 {
+ for i := 0; i < l.NumVertices(); i++ {
+ inside = inside != crosser.EdgeOrVertexCrossing(l.Vertex(i))
+ }
+ } else {
+ it := l.getEdgeIterator(l.NumVertices())
+ previousIndex := -2
+ for it.GetCandidates(origin, p); it.HasNext(); it.Next() {
+ ai := it.Index()
+ if previousIndex != ai-1 {
+ crosser.RestartAt(l.Vertex(ai))
+ }
+ previousIndex = ai
+ inside = inside != crosser.EdgeOrVertexCrossing(l.Vertex(ai+1))
+ }
+ }
+
+ return inside
+}
+
+/**
+ * Returns the shortest distance from a point P to this loop, given as the
+ * angle formed between P, the origin and the nearest point on the loop to P.
+ * This angle in radians is equivalent to the arclength along the unit sphere.
+ */
+func (l *Loop) GetDistance(p Point) s1.Angle {
+ normalized := Point{p.Normalize()}
+
+ // The furthest point from p on the sphere is its antipode, which is an
+ // angle of PI radians. This is an upper bound on the angle.
+ minDistance := math.Pi
+ for i := 0; i < l.NumVertices(); i++ {
+ minDistance = math.Min(minDistance, getDistance(normalized, l.Vertex(i), l.Vertex(i+1)).Radians())
+ }
+ return s1.Angle(minDistance)
+}
+
+/**
+ * Creates an edge index over the vertices, which by itself takes no time.
+ * Then the expected number of queries is used to determine whether brute
+ * force lookups are likely to be slower than really creating an index, and if
+ * so, we do so. Finally an iterator is returned that can be used to perform
+ * edge lookups.
+ */
+func (l *Loop) getEdgeIterator(expectedQueries int) *DataEdgeIterator {
+ if l.index == nil {
+ l.index = NewEdgeIndex(
+ func() int { return l.NumVertices() },
+ func(i int) Point { return l.Vertex(i) },
+ func(i int) Point { return l.Vertex(i + 1) },
+ )
+ }
+ l.index.PredictAdditionalCalls(expectedQueries)
+ return NewDataEdgeIterator(l.index)
+}
+
+/** Return true if this loop is valid. */
+func (l *Loop) IsValid() bool {
+ if l.NumVertices() < 3 {
+ fmt.Println("Degenerate loop")
+ return false
+ }
+
+ // All vertices must be unit length.
+ for i := 0; i < l.NumVertices(); i++ {
+ if !l.Vertex(i).IsUnit() {
+ fmt.Printf("Vertex %d is not unit length\n", i)
+ return false
+ }
+ }
+
+ // Loops are not allowed to have any duplicate vertices.
+ vmap := make(map[Point]int)
+ for i := 0; i < l.NumVertices(); i++ {
+ if previousVertexIndex, ok := vmap[l.Vertex(i)]; ok {
+ fmt.Printf("Duplicate vertices: %d and %d\n", previousVertexIndex, i)
+ return false
+ }
+ vmap[l.Vertex(i)] = i
+ }
+
+ // Non-adjacent edges are not allowed to intersect.
+ MAX_INTERSECTION_ERROR := 1e-15
+ crosses := false
+ it := l.getEdgeIterator(l.NumVertices())
+ for a1 := 0; a1 < l.NumVertices(); a1++ {
+ a2 := (a1 + 1) % l.NumVertices()
+ crosser := NewEdgeCrosser(l.Vertex(a1), l.Vertex(a2), l.Vertex(0))
+ previousIndex := -2
+ for it.GetCandidates(l.Vertex(a1), l.Vertex(a2)); it.HasNext(); it.Next() {
+ b1 := it.Index()
+ b2 := (b1 + 1) % l.NumVertices()
+ // If either 'a' index equals either 'b' index, then these two edges
+ // share a vertex. If a1==b1 then it must be the case that a2==b2, e.g.
+ // the two edges are the same. In that case, we skip the test, since we
+ // don't want to test an edge against itself. If a1==b2 or b1==a2 then
+ // we have one edge ending at the start of the other, or in other words,
+ // the edges share a vertex -- and in S2 space, where edges are always
+ // great circle segments on a sphere, edges can only intersect at most
+ // once, so we don't need to do further checks in that case either.
+ if a1 != b2 && a2 != b1 && a1 != b1 {
+ // WORKAROUND(shakusa, ericv): S2.robustCCW() currently
+ // requires arbitrary-precision arithmetic to be truly robust. That
+ // means it can give the wrong answers in cases where we are trying
+ // to determine edge intersections. The workaround is to ignore
+ // intersections between edge pairs where all four points are
+ // nearly colinear.
+ abc := angle(l.Vertex(a1), l.Vertex(a2), l.Vertex(b1))
+ abcNearlyLinear := approxEqualsNumber(abc, 0, MAX_INTERSECTION_ERROR) || approxEqualsNumber(abc, math.Pi, MAX_INTERSECTION_ERROR)
+ abd := angle(l.Vertex(a1), l.Vertex(a2), l.Vertex(b2))
+ abdNearlyLinear := approxEqualsNumber(abd, 0, MAX_INTERSECTION_ERROR) || approxEqualsNumber(abd, math.Pi, MAX_INTERSECTION_ERROR)
+ if abcNearlyLinear && abdNearlyLinear {
+ continue
+ }
+
+ if previousIndex != b1 {
+ crosser.RestartAt(l.Vertex(b1))
+ }
+
+ // Beware, this may return the loop is valid if there is a
+ // "vertex crossing".
+ // TODO(user): Fix that.
+ crosses = crosser.RobustCrossing(l.Vertex(b2)) > 0
+ previousIndex = b2
+ if crosses {
+ fmt.Printf("Edges %d and %d cross\n", a1, b1)
+ fmt.Printf("Edge locations in degrees: %s-%s and %s-%s\n",
+ LatLngFromPoint(l.Vertex(a1)).StringDegrees(),
+ LatLngFromPoint(l.Vertex(a2)).StringDegrees(),
+ LatLngFromPoint(l.Vertex(b1)).StringDegrees(),
+ LatLngFromPoint(l.Vertex(b2)).StringDegrees())
+ return false
+ }
+ }
+ }
+ }
+
+ return true
+}
+
+/**
+ * Static version of isValid(), to be used only when an S2Loop instance is not
+ * available, but validity of the points must be checked.
+ *
+ * @return true if the given loop is valid. Creates an instance of S2Loop and
+ * defers this call to {@link #isValid()}.
+ */
+func LoopIsValid(vertices []Point) bool {
+ return LoopFromPoints(vertices).IsValid()
+}
+
+func (l *Loop) ToString() string {
+ result := fmt.Sprintf("Loop, %d points. [", l.NumVertices())
+ for _, v := range l.vertices {
+ result = result + v.String() + " "
+ }
+ return result + "]"
+}
+
+func (l *Loop) initOrigin() {
+ // The bounding box does not need to be correct before calling this
+ // function, but it must at least contain vertex(1) since we need to
+ // do a Contains() test on this point below.
+ if !l.bound.ContainsLatLng(LatLngFromPoint(l.Vertex(1))) {
+ panic("Bounds needs to at least contain Vertex(1)")
+ }
+
+ // To ensure that every point is contained in exactly one face of a
+ // subdivision of the sphere, all containment tests are done by counting the
+ // edge crossings starting at a fixed point on the sphere (S2::Origin()).
+ // We need to know whether this point is inside or outside of the loop.
+ // We do this by first guessing that it is outside, and then seeing whether
+ // we get the correct containment result for vertex 1. If the result is
+ // incorrect, the origin must be inside the loop.
+ //
+ // A loop with consecutive vertices A,B,C contains vertex B if and only if
+ // the fixed vector R = S2::Ortho(B) is on the left side of the wedge ABC.
+ // The test below is written so that B is inside if C=R but not if A=R.
+
+ l.originInside = false // Initialize before calling Contains().
+ v1Inside := OrderedCCW(Point{l.Vertex(1).Ortho()}, l.Vertex(0), l.Vertex(2), l.Vertex(1))
+ if v1Inside != l.ContainsPoint(l.Vertex(1)) {
+ l.originInside = true
+ }
+}
+
+func (l *Loop) initBound() {
+ // The bounding rectangle of a loop is not necessarily the same as the
+ // bounding rectangle of its vertices. First, the loop may wrap entirely
+ // around the sphere (e.g. a loop that defines two revolutions of a
+ // candy-cane stripe). Second, the loop may include one or both poles.
+ // Note that a small clockwise loop near the equator contains both poles.
+
+ bounder := NewRectBounder()
+ for i := 0; i <= l.NumVertices(); i++ {
+ bounder.AddPoint(l.Vertex(i))
+ }
+ b := bounder.GetBound()
+ // Note that we need to initialize bound with a temporary value since
+ // contains() does a bounding rectangle check before doing anything else.
+ l.bound = FullRect()
+ if l.ContainsPoint(PointFromCoordsRaw(0, 0, 1)) {
+ b = Rect{r1.IntervalFromPointPair(b.Lat.Lo, math.Pi/2), s1.FullInterval()}
+ }
+ // If a loop contains the south pole, then either it wraps entirely
+ // around the sphere (full longitude range), or it also contains the
+ // north pole in which case b.lng().isFull() due to the test above.
+
+ if b.Lng.IsFull() && l.ContainsPoint(PointFromCoordsRaw(0, 0, -1)) {
+ b = Rect{r1.IntervalFromPointPair(-math.Pi/2, b.Lat.Hi), b.Lng}
+ }
+ l.bound = b
+}
+
+/**
+ * Return the index of a vertex at point "p", or -1 if not found. The return
+ * value is in the range 1..num_vertices_ if found.
+ */
+func (l *Loop) findVertex(p Point) int {
+ if l.vertexToIndex == nil {
+ l.vertexToIndex = make(map[Point]int)
+ for i := 1; i <= l.NumVertices(); i++ {
+ l.vertexToIndex[l.Vertex(i)] = i
+ }
+ }
+ if index, ok := l.vertexToIndex[p]; ok {
+ return index
+ }
+ return -1
+}
+
+/**
+ * This method encapsulates the common code for loop containment and
+ * intersection tests. It is used in three slightly different variations to
+ * implement contains(), intersects(), and containsOrCrosses().
+ *
+ * In a nutshell, this method checks all the edges of this loop (A) for
+ * intersection with all the edges of B. It returns -1 immediately if any edge
+ * intersections are found. Otherwise, if there are any shared vertices, it
+ * returns the minimum value of the given WedgeRelation for all such vertices
+ * (returning immediately if any wedge returns -1). Returns +1 if there are no
+ * intersections and no shared vertices.
+ */
+func (l *Loop) checkEdgeCrossings(b *Loop, relation WedgeRelation) int {
+ it := l.getEdgeIterator(b.NumVertices())
+ result := 1
+ // since 'this' usually has many more vertices than 'b', use the index on
+ // 'this' and loop over 'b'
+ for j := 0; j < b.NumVertices(); j++ {
+ crosser := NewEdgeCrosser(b.Vertex(j), b.Vertex(j+1), l.Vertex(0))
+ previousIndex := -2
+ for it.GetCandidates(b.Vertex(j), b.Vertex(j+1)); it.HasNext(); it.Next() {
+ i := it.Index()
+ if previousIndex != i-1 {
+ crosser.RestartAt(l.Vertex(i))
+ }
+ previousIndex = i
+ crossing := crosser.RobustCrossing(l.Vertex(i + 1))
+ if crossing < 0 {
+ continue
+ }
+ if crossing > 0 {
+ return -1 // There is a proper edge crossing.
+ }
+ if l.Vertex(i + 1).Equals(b.Vertex(j + 1)) {
+ result = min(result, relation.Test(
+ l.Vertex(i),
+ l.Vertex(i+1),
+ l.Vertex(i+2),
+ b.Vertex(j),
+ b.Vertex(j+2),
+ ))
+ if result < 0 {
+ return result
+ }
+ }
+ }
+ }
+ return result
+}
diff --git a/s2/loop_test.go b/s2/loop_test.go
new file mode 100644
index 00000000..c2a7f372
--- /dev/null
+++ b/s2/loop_test.go
@@ -0,0 +1,481 @@
+package s2
+
+import (
+ "math"
+ "math/rand"
+ "testing"
+
+ "github.com/golang/geo/r1"
+ "github.com/golang/geo/s1"
+)
+
+var (
+ // A stripe that slightly over-wraps the equator.
+ candyCane *Loop = makeLoop("-20:150, -20:-70, 0:70, 10:-150, 10:70, -10:-70")
+
+ // A small clockwise loop in the northern & eastern hemisperes.
+ smallNeCw *Loop = makeLoop("35:20, 45:20, 40:25")
+
+ // Loop around the north pole at 80 degrees.
+ arctic80 *Loop = makeLoop("80:-150, 80:-30, 80:90")
+
+ // Loop around the south pole at 80 degrees.
+ antarctic80 *Loop = makeLoop("-80:120, -80:0, -80:-120")
+
+ // The northern hemisphere, defined using two pairs of antipodal points.
+ northHemi *Loop = makeLoop("0:-180, 0:-90, 0:0, 0:90")
+
+ // The northern hemisphere, defined using three points 120 degrees apart.
+ northHemi3 *Loop = makeLoop("0:-180, 0:-60, 0:60")
+
+ // The western hemisphere, defined using two pairs of antipodal points.
+ westHemi *Loop = makeLoop("0:-180, -90:0, 0:0, 90:0")
+
+ // The "near" hemisphere, defined using two pairs of antipodal points.
+ nearHemi *Loop = makeLoop("0:-90, -90:0, 0:90, 90:0")
+
+ // A diamond-shaped loop around the point 0:180.
+ loopA *Loop = makeLoop("0:178, -1:180, 0:-179, 1:-180")
+
+ // Another diamond-shaped loop around the point 0:180.
+ loopB *Loop = makeLoop("0:179, -1:180, 0:-178, 1:-180")
+
+ // The intersection of A and B.
+ aIntersectB *Loop = makeLoop("0:179, -1:180, 0:-179, 1:-180")
+
+ // The union of A and B.
+ aUnionB *Loop = makeLoop("0:178, -1:180, 0:-178, 1:-180")
+
+ // A minus B (concave)
+ aMinusB *Loop = makeLoop("0:178, -1:180, 0:179, 1:-180")
+
+ // B minus A (concave)
+ bMinusA *Loop = makeLoop("0:-179, -1:180, 0:-178, 1:-180")
+
+ // A self-crossing loop with a duplicated vertex
+ bowtie *Loop = makeLoop("0:0, 2:0, 1:1, 0:2, 2:2, 1:1")
+
+ // Initialized below.
+ southHemi *Loop
+ eastHemi *Loop
+ farHemi *Loop
+)
+
+func init() {
+ southHemi = makeLoop("0:-180, 0:-90, 0:0, 0:90")
+ southHemi.Invert()
+ eastHemi = makeLoop("0:-180, -90:0, 0:0, 90:0")
+ eastHemi.Invert()
+ farHemi = makeLoop("0:-90, -90:0, 0:90, 90:0")
+ farHemi.Invert()
+}
+
+func makeLoop(s string) *Loop {
+ points := parsePoints(s)
+ return LoopFromPoints(points)
+}
+
+func TestLoopBounds(t *testing.T) {
+ if !(candyCane.RectBound().Lng.IsFull()) {
+ t.Fatal("ttttt")
+ }
+ if !(s1.Angle(candyCane.RectBound().Lat.Lo).Degrees() < -20) {
+ t.Fatal("")
+ }
+ if !(s1.Angle(candyCane.RectBound().Lat.Hi).Degrees() > 10) {
+ t.Fatal("")
+ }
+ if !(smallNeCw.RectBound().IsFull()) {
+ t.Fatal("")
+ }
+ if arctic80.RectBound() != RectFromLatLngLoHi(LatLngFromDegrees(80, -180), LatLngFromDegrees(90, 180)) {
+ t.Fatal("")
+ }
+ if antarctic80.RectBound() != RectFromLatLngLoHi(LatLngFromDegrees(-90, -180), LatLngFromDegrees(-80, 180)) {
+ t.Fatal("")
+ }
+
+ arctic80.Invert()
+ // The highest latitude of each edge is attained at its midpoint.
+ mid := Point{arctic80.Vertex(0).Add(arctic80.Vertex(1).Vector).Mul(0.5)}
+ if math.Abs(s1.Angle(arctic80.RectBound().Lat.Hi).Radians()-LatLngFromPoint(mid).Lat.Radians()) > EPSILON {
+ t.Fatal("")
+ }
+ arctic80.Invert()
+
+ if !(southHemi.RectBound().Lng.IsFull()) {
+ t.Fatal("")
+ }
+ if !(southHemi.RectBound().Lat == r1.IntervalFromPointPair(-math.Pi/2, 0)) {
+ t.Fatal("")
+ }
+}
+
+func TestLoopAreaCentroid(t *testing.T) {
+ if !(northHemi.GetArea() == 2*math.Pi) {
+ t.Fatal("")
+ }
+ if !(eastHemi.GetArea() == 2*math.Pi) {
+ t.Fatal("")
+ }
+
+ // Construct spherical caps of random height, and approximate their boundary
+ // with closely spaces vertices. Then check that the area and centroid are
+ // correct.
+
+ for i := 0; i < 100; i++ {
+ // Choose a coordinate frame for the spherical cap.
+ x := randomPoint()
+ y := x.Cross(randomPoint().Vector).Normalize()
+ z := x.Cross(y).Normalize()
+
+ // Given two points at latitude phi and whose longitudes differ by dtheta,
+ // the geodesic between the two points has a maximum latitude of
+ // atan(tan(phi) / cos(dtheta/2)). This can be derived by positioning
+ // the two points at (-dtheta/2, phi) and (dtheta/2, phi).
+ //
+ // We want to position the vertices close enough together so that their
+ // maximum distance from the boundary of the spherical cap is kMaxDist.
+ // Thus we want fabs(atan(tan(phi) / cos(dtheta/2)) - phi) <= kMaxDist.
+ kMaxDist := 1e-6
+ height := 2 * rand.Float64()
+ phi := math.Asin(1 - height)
+ maxDtheta := 2 * math.Acos(math.Tan(math.Abs(phi))/math.Tan(math.Abs(phi)+kMaxDist))
+ maxDtheta = math.Min(math.Pi, maxDtheta) // At least 3 vertices.
+
+ vertices := []Point{}
+ for theta := 0.0; theta < 2*math.Pi; theta += rand.Float64() * maxDtheta {
+
+ xCosThetaCosPhi := x.Mul((math.Cos(theta) * math.Cos(phi)))
+ ySinThetaCosPhi := y.Mul((math.Sin(theta) * math.Cos(phi)))
+ zSinPhi := z.Mul(math.Sin(phi))
+
+ sum := Point{xCosThetaCosPhi.Add(ySinThetaCosPhi.Add(zSinPhi))}
+
+ vertices = append(vertices, sum)
+ }
+
+ loop := LoopFromPoints(vertices)
+ areaCentroid := loop.GetAreaAndCentroid()
+
+ area := loop.GetArea()
+ centroid := loop.GetCentroid()
+ expectedArea := 2 * math.Pi * height
+ if areaCentroid.GetArea() != area {
+ t.Fatal("")
+ }
+ if !centroid.Equals(areaCentroid.GetCentroid()) {
+ t.Fatal("")
+ }
+ if !(math.Abs(area-expectedArea) <= 2*math.Pi*kMaxDist) {
+ t.Fatal("")
+ }
+
+ // high probability
+ if !(math.Abs(area-expectedArea) >= 0.01*kMaxDist) {
+ t.Fatal("")
+ }
+
+ expectedCentroid := z.Mul(expectedArea * (1 - 0.5*height))
+
+ if !(centroid.Sub(expectedCentroid).Norm() <= 2*kMaxDist) {
+ t.Fatal("")
+ }
+ }
+}
+
+func rotate(loop *Loop) *Loop {
+ vertices := make([]Point, 0, loop.NumVertices())
+ vertices = append(vertices, loop.vertices[1:loop.NumVertices()]...)
+ vertices = append(vertices, loop.vertices[0])
+ return LoopFromPoints(vertices)
+}
+
+func TestLoopContains(t *testing.T) {
+ if !(candyCane.ContainsPoint(PointFromLatLng(LatLngFromDegrees(5, 71)))) {
+ t.FailNow()
+ }
+ for i := 0; i < 4; i++ {
+ if !(northHemi.ContainsPoint(PointFromCoordsRaw(0, 0, 1))) {
+ t.FailNow()
+ }
+ if !(!northHemi.ContainsPoint(PointFromCoordsRaw(0, 0, -1))) {
+ t.FailNow()
+ }
+ if !(!southHemi.ContainsPoint(PointFromCoordsRaw(0, 0, 1))) {
+ t.FailNow()
+ }
+ if !(southHemi.ContainsPoint(PointFromCoordsRaw(0, 0, -1))) {
+ t.FailNow()
+ }
+ if !(!westHemi.ContainsPoint(PointFromCoordsRaw(0, 1, 0))) {
+ t.FailNow()
+ }
+ if !(westHemi.ContainsPoint(PointFromCoordsRaw(0, -1, 0))) {
+ t.FailNow()
+ }
+ if !(eastHemi.ContainsPoint(PointFromCoordsRaw(0, 1, 0))) {
+ t.FailNow()
+ }
+ if !(!eastHemi.ContainsPoint(PointFromCoordsRaw(0, -1, 0))) {
+ t.FailNow()
+ }
+ northHemi = rotate(northHemi)
+ southHemi = rotate(southHemi)
+ eastHemi = rotate(eastHemi)
+ westHemi = rotate(westHemi)
+ }
+
+ // This code checks each cell vertex is contained by exactly one of
+ // the adjacent cells.
+ for level := 0; level < 3; level++ {
+ loops := []*Loop{}
+ loopVertices := []Point{}
+ points := make(map[Point]bool)
+
+ for id := CellIDBegin(level); id != CellIDEnd(level); id = id.Next() {
+ cell := CellFromCellID(id)
+ points[cell.Id().Point()] = true
+ for k := 0; k < 4; k++ {
+ loopVertices = append(loopVertices, cell.Vertex(k))
+ points[cell.Vertex(k)] = true
+ }
+ loops = append(loops, LoopFromPoints(loopVertices))
+ loopVertices = []Point{}
+ }
+ for point := range points {
+ count := 0
+ for _, loop := range loops {
+ if loop.ContainsPoint(point) {
+ count++
+ }
+ }
+ if count != 1 {
+ t.Fatalf("Failed at level %d with count %d, loops:%d, points:%d", level, count, len(loops), len(points))
+ }
+ }
+ }
+}
+
+func testLoopRelation(t *testing.T, a, b *Loop, containsOrCrosses int, intersects, nestable bool) {
+ if a.ContainsLoop(b) != (containsOrCrosses == 1) {
+ if containsOrCrosses == 1 {
+ t.Fatalf("loop should be contained or crossing")
+ } else {
+ t.Fatalf("loop should not be contained or crossing")
+ }
+ }
+ if a.IntersectsLoop(b) != intersects {
+ if intersects {
+ t.Fatalf("loops should intersect")
+ } else {
+ t.Fatalf("loops should not intersect")
+ }
+ }
+ if nestable {
+ if a.ContainsNested(b) != a.ContainsLoop(b) {
+ t.Fatalf("loops should be nested")
+ }
+ }
+ if containsOrCrosses >= -1 {
+ if a.ContainsOrCrosses(b) != containsOrCrosses {
+ t.Fatalf("loops should contain or cross %d", containsOrCrosses)
+ }
+ }
+}
+
+func TestLoopRelations(t *testing.T) {
+ testLoopRelation(t, northHemi, northHemi, 1, true, false)
+ testLoopRelation(t, northHemi, southHemi, 0, false, false)
+ testLoopRelation(t, northHemi, eastHemi, -1, true, false)
+ testLoopRelation(t, northHemi, arctic80, 1, true, true)
+ testLoopRelation(t, northHemi, antarctic80, 0, false, true)
+ testLoopRelation(t, northHemi, candyCane, -1, true, false)
+
+ // // We can't compare northHemi3 vs. northHemi or southHemi.
+ testLoopRelation(t, northHemi3, northHemi3, 1, true, false)
+ testLoopRelation(t, northHemi3, eastHemi, -1, true, false)
+ testLoopRelation(t, northHemi3, arctic80, 1, true, true)
+ testLoopRelation(t, northHemi3, antarctic80, 0, false, true)
+ testLoopRelation(t, northHemi3, candyCane, -1, true, false)
+
+ testLoopRelation(t, southHemi, northHemi, 0, false, false)
+ testLoopRelation(t, southHemi, southHemi, 1, true, false)
+ testLoopRelation(t, southHemi, farHemi, -1, true, false)
+ testLoopRelation(t, southHemi, arctic80, 0, false, true)
+ testLoopRelation(t, southHemi, antarctic80, 1, true, true)
+ testLoopRelation(t, southHemi, candyCane, -1, true, false)
+
+ testLoopRelation(t, candyCane, northHemi, -1, true, false)
+ testLoopRelation(t, candyCane, southHemi, -1, true, false)
+ testLoopRelation(t, candyCane, arctic80, 0, false, true)
+ testLoopRelation(t, candyCane, antarctic80, 0, false, true)
+ testLoopRelation(t, candyCane, candyCane, 1, true, false)
+
+ testLoopRelation(t, nearHemi, westHemi, -1, true, false)
+
+ testLoopRelation(t, smallNeCw, southHemi, 1, true, false)
+ testLoopRelation(t, smallNeCw, westHemi, 1, true, false)
+ testLoopRelation(t, smallNeCw, northHemi, -2, true, false)
+ testLoopRelation(t, smallNeCw, eastHemi, -2, true, false)
+
+ testLoopRelation(t, loopA, loopA, 1, true, false)
+ testLoopRelation(t, loopA, loopB, -1, true, false)
+ testLoopRelation(t, loopA, aIntersectB, 1, true, false)
+ testLoopRelation(t, loopA, aUnionB, 0, true, false)
+ testLoopRelation(t, loopA, aMinusB, 1, true, false)
+ testLoopRelation(t, loopA, bMinusA, 0, false, false)
+
+ testLoopRelation(t, loopB, loopA, -1, true, false)
+ testLoopRelation(t, loopB, loopB, 1, true, false)
+ testLoopRelation(t, loopB, aIntersectB, 1, true, false)
+ testLoopRelation(t, loopB, aUnionB, 0, true, false)
+ testLoopRelation(t, loopB, aMinusB, 0, false, false)
+ testLoopRelation(t, loopB, bMinusA, 1, true, false)
+
+ testLoopRelation(t, aIntersectB, loopA, 0, true, false)
+ testLoopRelation(t, aIntersectB, loopB, 0, true, false)
+ testLoopRelation(t, aIntersectB, aIntersectB, 1, true, false)
+ testLoopRelation(t, aIntersectB, aUnionB, 0, true, true)
+ testLoopRelation(t, aIntersectB, aMinusB, 0, false, false)
+ testLoopRelation(t, aIntersectB, bMinusA, 0, false, false)
+
+ testLoopRelation(t, aUnionB, loopA, 1, true, false)
+ testLoopRelation(t, aUnionB, loopB, 1, true, false)
+ testLoopRelation(t, aUnionB, aIntersectB, 1, true, true)
+ testLoopRelation(t, aUnionB, aUnionB, 1, true, false)
+ testLoopRelation(t, aUnionB, aMinusB, 1, true, false)
+ testLoopRelation(t, aUnionB, bMinusA, 1, true, false)
+
+ testLoopRelation(t, aMinusB, loopA, 0, true, false)
+ testLoopRelation(t, aMinusB, loopB, 0, false, false)
+ testLoopRelation(t, aMinusB, aIntersectB, 0, false, false)
+ testLoopRelation(t, aMinusB, aUnionB, 0, true, false)
+ testLoopRelation(t, aMinusB, aMinusB, 1, true, false)
+ testLoopRelation(t, aMinusB, bMinusA, 0, false, true)
+
+ testLoopRelation(t, bMinusA, loopA, 0, false, false)
+ testLoopRelation(t, bMinusA, loopB, 0, true, false)
+ testLoopRelation(t, bMinusA, aIntersectB, 0, false, false)
+ testLoopRelation(t, bMinusA, aUnionB, 0, true, false)
+ testLoopRelation(t, bMinusA, aMinusB, 0, false, true)
+ testLoopRelation(t, bMinusA, bMinusA, 1, true, false)
+}
+
+/**
+ * Tests that nearly colinear points pass S2Loop.isValid()
+ */
+func TestLoopRoundingError(t *testing.T) {
+ points := []Point{
+ PointFromCoordsRaw(-0.9190364081111774, 0.17231932652084575, 0.35451111445694833),
+ PointFromCoordsRaw(-0.92130667053206, 0.17274500072476123, 0.3483578383756171),
+ PointFromCoordsRaw(-0.9257244057938284, 0.17357332608634282, 0.3360158106235289),
+ PointFromCoordsRaw(-0.9278712595449962, 0.17397586116468677, 0.32982923679138537),
+ }
+ loop := LoopFromPoints(points)
+ if !loop.IsValid() {
+ t.Errorf("loop should be valid")
+ }
+}
+
+func TestLoopIsValid(t *testing.T) {
+ if !loopA.IsValid() {
+ t.Errorf("loopA should be valid")
+ }
+ if !loopB.IsValid() {
+ t.Errorf("loopB should be valid")
+ }
+ if bowtie.IsValid() {
+ t.Errorf("bowtie should not be valid")
+ }
+}
+
+/**
+ * Tests {@link S2Loop#compareTo(S2Loop)}.
+ */
+func TestLoopComparisons(t *testing.T) {
+ abc := makeLoop("0:1, 0:2, 1:2")
+ abcd := makeLoop("0:1, 0:2, 1:2, 1:1")
+ abcde := makeLoop("0:1, 0:2, 1:2, 1:1, 1:0")
+ if !(abc.CompareTo(abcd) < 0) {
+ t.FailNow()
+ }
+ if !(abc.CompareTo(abcde) < 0) {
+ t.FailNow()
+ }
+ if !(abcd.CompareTo(abcde) < 0) {
+ t.FailNow()
+ }
+ if !(abcd.CompareTo(abc) > 0) {
+ t.FailNow()
+ }
+ if !(abcde.CompareTo(abc) > 0) {
+ t.FailNow()
+ }
+ if !(abcde.CompareTo(abcd) > 0) {
+ t.FailNow()
+ }
+
+ bcda := makeLoop("0:2, 1:2, 1:1, 0:1")
+ if 0 != abcd.CompareTo(bcda) {
+ t.FailNow()
+ }
+ if 0 != bcda.CompareTo(abcd) {
+ t.FailNow()
+ }
+
+ wxyz := makeLoop("10:11, 10:12, 11:12, 11:11")
+ if !(abcd.CompareTo(wxyz) > 0) {
+ t.FailNow()
+ }
+ if !(wxyz.CompareTo(abcd) < 0) {
+ t.FailNow()
+ }
+}
+func TestLoopGetDistance(t *testing.T) {
+ // Error margin since we're doing numerical computations
+ epsilon := 1e-15
+
+ // A square with (lat,lng) vertices (0,1), (1,1), (1,2) and (0,2)
+ // Tests the case where the shortest distance is along a normal to an edge,
+ // onto a vertex
+ s1 := makeLoop("0:1, 1:1, 1:2, 0:2")
+
+ // A square with (lat,lng) vertices (-1,1), (1,1), (1,2) and (-1,2)
+ // Tests the case where the shortest distance is along a normal to an edge,
+ // not onto a vertex
+ s2 := makeLoop("-1:1, 1:1, 1:2, -1:2")
+
+ // A diamond with (lat,lng) vertices (1,0), (2,1), (3,0) and (2,-1)
+ // Test the case where the shortest distance is NOT along a normal to an
+ // edge
+ s3 := makeLoop("1:0, 2:1, 3:0, 2:-1")
+
+ // All the vertices should be distance 0
+ for i := 0; i < s1.NumVertices(); i++ {
+ if math.Abs(s1.GetDistance(s1.Vertex(i)).Radians()) > epsilon {
+ t.FailNow()
+ }
+ }
+
+ // A point on one of the edges should be distance 0
+ if math.Abs(s1.GetDistance(PointFromLatLng(LatLngFromDegrees(0.5, 1))).Radians()) > epsilon {
+ t.FailNow()
+ }
+
+ // In all three cases, the closest point to the origin is (0,1), which is at
+ // a distance of 1 degree.
+ // Note: all of these are intentionally distances measured along the
+ // equator, since that makes the math significantly simpler. Otherwise, the
+ // distance wouldn't actually be 1 degree.
+ origin := PointFromLatLng(LatLngFromDegrees(0, 0))
+ if math.Abs(1-s1.GetDistance(origin).Degrees()) > epsilon {
+ t.FailNow()
+ }
+ if math.Abs(1-s2.GetDistance(origin).Degrees()) > epsilon {
+ t.FailNow()
+ }
+ if math.Abs(1-s3.GetDistance(origin).Degrees()) > epsilon {
+ t.FailNow()
+ }
+}
diff --git a/s2/point.go b/s2/point.go
index 9f90db92..44f4a23f 100644
--- a/s2/point.go
+++ b/s2/point.go
@@ -17,6 +17,7 @@ limitations under the License.
package s2
import (
+ "fmt"
"math"
"github.com/golang/geo/r3"
@@ -83,6 +84,18 @@ func PointFromCoords(x, y, z float64) Point {
return Point{r3.Vector{x, y, z}.Normalize()}
}
+func PointFromCoordsRaw(x, y, z float64) Point {
+ return Point{r3.Vector{x, y, z}}
+}
+
+// PointFromLatLng returns an Point for the given LatLng.
+func PointFromLatLng(ll LatLng) Point {
+ phi := ll.Lat.Radians()
+ theta := ll.Lng.Radians()
+ cosphi := math.Cos(phi)
+ return PointFromCoords(math.Cos(theta)*cosphi, math.Sin(theta)*cosphi, math.Sin(phi))
+}
+
// OriginPoint returns a unique "origin" on the sphere for operations that need a fixed
// reference point. In particular, this is the "point at infinity" used for
// point-in-polygon testing (by counting the number of edge crossings).
@@ -232,42 +245,53 @@ func RobustSign(a, b, c Point) Direction {
return Indeterminate
}
-// OrderedCCW returns 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 the
-// CCW ordering around O that starts at A, or equivalently, whether B is
-// contained in the range of angles (inclusive) that starts at A and extends
-// CCW to C. Properties:
-//
-// (1) If OrderedCCW(a,b,c,o) && OrderedCCW(b,a,c,o), then a == b
-// (2) If OrderedCCW(a,b,c,o) && OrderedCCW(a,c,b,o), then b == c
-// (3) If OrderedCCW(a,b,c,o) && OrderedCCW(c,b,a,o), then a == b == c
-// (4) If a == b or b == c, then OrderedCCW(a,b,c,o) is true
-// (5) Otherwise if a == c, then OrderedCCW(a,b,c,o) is false
-func OrderedCCW(a, b, c, o Point) bool {
- sum := 0
- if RobustSign(b, o, a) != Clockwise {
- sum++
+// Distance returns the angle between two points.
+func (p Point) Distance(b Point) s1.Angle {
+ return p.Vector.Angle(b.Vector)
+}
+
+// ApproxEqual reports if the two points are similar enough to be equal.
+func (p Point) ApproxEquals(other Point, maxError float64) bool {
+ return p.Vector.Angle(other.Vector).Radians() <= maxError
+}
+
+func (p Point) Equals(other Point) bool {
+ return p.X == other.X && p.Y == other.Y && p.Z == other.Z
+}
+
+func (p Point) LessThan(vb Point) bool {
+ if p.X < vb.X {
+ return true
}
- if RobustSign(c, o, b) != Clockwise {
- sum++
+ if vb.X < p.X {
+ return false
}
- if RobustSign(a, o, c) == CounterClockwise {
- sum++
+ if p.Y < vb.Y {
+ return true
}
- return sum >= 2
+ if vb.Y < p.Y {
+ return false
+ }
+ if p.Z < vb.Z {
+ return true
+ }
+ return false
}
-// Distance returns the angle between two points.
-func (p Point) Distance(b Point) s1.Angle {
- return p.Vector.Angle(b.Vector)
+func (p Point) CompareTo(other Point) int {
+ if p.LessThan(other) {
+ return -1
+ } else {
+ if p.Equals(other) {
+ return 0
+ }
+ return 1
+ }
}
-// ApproxEqual reports if the two points are similar enough to be equal.
-func (p Point) ApproxEqual(other Point) bool {
- const epsilon = 1e-14
- return p.Vector.Angle(other.Vector) <= epsilon
+func (p Point) DegreesString() string {
+ latLng := LatLngFromPoint(p)
+ return fmt.Sprintf("(%f, %f)", latLng.Lat.Degrees(), latLng.Lng.Degrees())
}
// PointArea returns the area on the unit sphere for the triangle defined by the
diff --git a/s2/point_test.go b/s2/point_test.go
index 20a4db23..27703b23 100644
--- a/s2/point_test.go
+++ b/s2/point_test.go
@@ -263,7 +263,7 @@ func TestPointDistance(t *testing.T) {
}
}
-func TestApproxEqual(t *testing.T) {
+func TestApproxEquals(t *testing.T) {
epsilon := 1e-14
tests := []struct {
x1, y1, z1 float64
@@ -286,8 +286,8 @@ func TestApproxEqual(t *testing.T) {
for _, test := range tests {
p1 := PointFromCoords(test.x1, test.y1, test.z1)
p2 := PointFromCoords(test.x2, test.y2, test.z2)
- if got := p1.ApproxEqual(p2); got != test.want {
- t.Errorf("%v.ApproxEqual(%v), got %v want %v", p1, p2, got, test.want)
+ if got := p1.ApproxEquals(p2, EPSILON); got != test.want {
+ t.Errorf("%v.ApproxEquals(%v), got %v want %v", p1, p2, got, test.want)
}
}
}
diff --git a/s2/polyline.go b/s2/polyline.go
new file mode 100644
index 00000000..b9f4e06d
--- /dev/null
+++ b/s2/polyline.go
@@ -0,0 +1,127 @@
+/*
+Copyright 2015 Google Inc. All rights reserved.
+
+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 s2
+
+import (
+ "fmt"
+
+ "github.com/golang/geo/s1"
+)
+
+/**
+ * A Polyline represents a sequence of zero or more vertices connected by
+ * straight edges (geodesics). Edges of length 0 and 180 degrees are not
+ * allowed, i.e. adjacent vertices should not be identical or antipodal.
+ *
+ * Note: Polylines do not have a Contains(S2Point) method, because
+ * "containment" is not numerically well-defined except at the polyline
+ * vertices.
+ */
+type Polyline struct {
+ Vertices []Point
+}
+
+func PolylineFromPoints(points []Point) Polyline {
+ return Polyline{points}
+}
+
+// Return true if the given vertices form a valid polyline.
+func (p Polyline) IsValid() bool {
+ // All vertices must be unit length.
+ n := len(p.Vertices)
+ for i := 0; i < n; i++ {
+ if !p.Vertices[i].IsUnit() {
+ fmt.Printf("Vertex %d is not unit length", i)
+ return false
+ }
+ }
+ // Adjacent vertices must not be identical or antipodal.
+ for i := 0; i < n; i++ {
+ if p.Vertices[i-1].ApproxEquals(p.Vertices[i], EPSILON) || p.Vertices[i-1].ApproxEquals(Point{p.Vertices[i].Neg()}, EPSILON) {
+ fmt.Printf("Vertices %d and %d are identical or antipodal", (i - 1), i)
+ return false
+ }
+ }
+ return true
+}
+
+func (p Polyline) NumVertices() int { return len(p.Vertices) }
+
+func (p Polyline) Vertex(k int) Point { return p.Vertices[k] }
+
+// Return the angle corresponding to the total arclength of the polyline on a unit sphere.
+func (p Polyline) GetArclengthAngle() s1.Angle {
+ var lengthSum s1.Angle = 0
+ for i := 1; i < p.NumVertices(); i++ {
+ lengthSum += p.Vertices[i-1].Angle(p.Vertices[i].Vector)
+ }
+ return lengthSum
+}
+
+// CapBound returns a bounding spherical cap. This is not guaranteed to be exact.
+func (p Polyline) CapBound() Cap {
+ return p.RectBound().CapBound()
+}
+
+// RectBound returns a bounding latitude-longitude rectangle that contains
+// the region. The bounds are not guaranteed to be tight.
+func (p Polyline) RectBound() Rect {
+ rb := NewRectBounder()
+ for i := 0; i < p.NumVertices(); i++ {
+ rb.AddPoint(p.Vertex(i))
+ }
+ return rb.GetBound()
+}
+
+// ContainsCell reports whether the region completely contains the given region.
+// It returns false if containment could not be determined.
+func (p Polyline) ContainsCell(cell Cell) bool {
+ return false
+}
+
+// IntersectsCell reports whether the region intersects the given cell or
+// if intersection could not be determined. It returns false if the region
+// does not intersect.
+func (p Polyline) IntersectsCell(cell Cell) bool {
+ if p.NumVertices() == 0 {
+ return false
+ }
+
+ // We only need to check whether the cell contains vertex 0 for correctness,
+ // but these tests are cheap compared to edge crossings so we might as well
+ // check all the vertices.
+ for i := 0; i < p.NumVertices(); i++ {
+ if cell.ContainsPoint(p.Vertex(i)) {
+ return true
+ }
+ }
+
+ cellVertices := make([]Point, 4)
+ for i := 0; i < 4; i++ {
+ cellVertices[i] = cell.Vertex(i)
+ }
+ for j := 0; j < 4; j++ {
+ ec := NewEdgeCrosser(cellVertices[j], cellVertices[(j+1)&3], p.Vertex(0))
+ for i := 1; i < p.NumVertices(); i++ {
+ if ec.RobustCrossing(p.Vertex(i)) >= 0 {
+ // There is a proper crossing, or two vertices were the same.
+ return true
+ }
+ }
+ }
+ return false
+}
diff --git a/s2/projections.go b/s2/projections.go
new file mode 100644
index 00000000..45d6e29b
--- /dev/null
+++ b/s2/projections.go
@@ -0,0 +1,154 @@
+package s2
+
+import (
+ "math"
+)
+
+var (
+ // Note that other parts of this GO port have been hardcoded to QUADRATIC projection
+
+ // Uncomment the desirect projection type
+ // S2_PROJECTION Projections = S2_LINEAR_PROJECTION{}
+ // S2_PROJECTION Projections = S2_TAN_PROJECTION{}
+ S2_PROJECTION Projections = S2_QUADRATIC_PROJECTION{}
+)
+
+type Projections interface {
+ MIN_AREA() Metric
+ MAX_AREA() Metric
+ AVG_AREA() Metric
+
+ MIN_ANGLE_SPAN() Metric
+ MAX_ANGLE_SPAN() Metric
+ AVG_ANGLE_SPAN() Metric
+
+ MIN_WIDTH() Metric
+ MAX_WIDTH() Metric
+ AVG_WIDTH() Metric
+
+ MIN_EDGE() Metric
+ MAX_EDGE() Metric
+ AVG_EDGE() Metric
+
+ MIN_DIAG() Metric
+ MAX_DIAG() Metric
+ AVG_DIAG() Metric
+
+ MAX_EDGE_ASPECT() float64
+ MAX_DIAG_ASPECT() float64
+}
+
+type S2_PROJECTION_COMMON struct{}
+
+func (m S2_PROJECTION_COMMON) AVG_AREA() Metric {
+ return Metric{
+ math.Pi / 6, // 0.524
+ 2,
+ }
+}
+
+func (m S2_PROJECTION_COMMON) AVG_ANGLE_SPAN() Metric {
+ return Metric{
+ math.Pi / 4, // 0.785
+ 1,
+ }
+}
+
+func (m S2_PROJECTION_COMMON) MAX_DIAG_ASPECT() float64 {
+ return math.Sqrt(3) // 1.732
+}
+
+type S2_QUADRATIC_PROJECTION struct{ S2_PROJECTION_COMMON }
+
+func (m S2_QUADRATIC_PROJECTION) MIN_AREA() Metric {
+ return Metric{
+ 2 * math.Sqrt2 / 9, // 0.314
+ 2,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MAX_AREA() Metric {
+ return Metric{
+ 0.65894981424079037, // 0.659
+ 2,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MIN_ANGLE_SPAN() Metric {
+ return Metric{
+ 2.0 / 3.0, // 0.667
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MAX_ANGLE_SPAN() Metric {
+ return Metric{
+ 0.85244858959960922, // 0.852
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MIN_WIDTH() Metric {
+ return Metric{
+ math.Sqrt2 / 3, // 0.471
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MAX_WIDTH() Metric {
+ return m.MAX_ANGLE_SPAN()
+}
+
+func (m S2_QUADRATIC_PROJECTION) AVG_WIDTH() Metric {
+ return Metric{
+ 0.71726183644304969, // 0.717
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MIN_EDGE() Metric {
+ return Metric{
+ math.Sqrt2 / 3, // 0.471
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MAX_EDGE() Metric {
+ return m.MAX_ANGLE_SPAN()
+}
+
+func (m S2_QUADRATIC_PROJECTION) AVG_EDGE() Metric {
+ return Metric{
+ 0.72960687319305303, // 0.730
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MIN_DIAG() Metric {
+ return Metric{
+ 4 * math.Sqrt2 / 9, // // 0.629
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MAX_DIAG() Metric {
+ return Metric{
+ 1.2193272972170106, // 1.219
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) AVG_DIAG() Metric {
+ return Metric{
+ 1.03021136949923584, // 1.030
+ 1,
+ }
+}
+
+func (m S2_QUADRATIC_PROJECTION) MAX_EDGE_ASPECT() float64 {
+ return 1.44261527445268292 // 1.443
+}
+
+// TODO: Add S2_LINEAR_PROJECTION, S2_TAN_PROJECTION
+// type S2_LINEAR_PROJECTION struct{ S2_PROJECTION_COMMON }
+// type S2_TAN_PROJECTION struct{ S2_PROJECTION_COMMON }
diff --git a/s2/rect.go b/s2/rect.go
index 56af72c3..5d04744b 100644
--- a/s2/rect.go
+++ b/s2/rect.go
@@ -49,6 +49,29 @@ func RectFromLatLng(p LatLng) Rect {
}
}
+func RectFromLatLngLoHi(lo, hi LatLng) Rect {
+ // assert (p1.isValid() && p2.isValid());
+ return Rect{
+ Lat: r1.IntervalFromEndpoints(lo.Lat.Radians(), hi.Lat.Radians()),
+ Lng: s1.IntervalFromEndpoints(lo.Lng.Radians(), hi.Lng.Radians()),
+ }
+}
+
+/**
+ * Convenience method to construct the minimal bounding rectangle containing
+ * the two given points. This is equivalent to starting with an empty
+ * rectangle and calling AddPoint() twice. Note that it is different than the
+ * S2LatLngRect(lo, hi) constructor, where the first point is always used as
+ * the lower-left corner of the resulting rectangle.
+ */
+func RectFromLatLngPointPair(p1, p2 LatLng) Rect {
+ // assert (p1.isValid() && p2.isValid());
+ return Rect{
+ Lat: r1.IntervalFromPointPair(p1.Lat.Radians(), p2.Lat.Radians()),
+ Lng: s1.IntervalFromPointPair(p1.Lng.Radians(), p2.Lng.Radians()),
+ }
+}
+
// RectFromCenterSize constructs a rectangle with the given size and center.
// center needs to be normalized, but size does not. The latitude
// interval of the result is clamped to [-90,90] degrees, and the longitude
@@ -92,6 +115,23 @@ func (r Rect) Hi() LatLng {
return LatLng{s1.Angle(r.Lat.Hi) * s1.Radian, s1.Angle(r.Lng.Hi) * s1.Radian}
}
+// Return the k-th vertex of the rectangle (k = 0,1,2,3) in CCW order.
+func (r Rect) Vertex(k int) LatLng {
+ // Return the points in CCW order (SW, SE, NE, NW).
+ switch k {
+ case 0:
+ return LatLng{s1.Angle(r.Lat.Lo), s1.Angle(r.Lng.Lo)}
+ case 1:
+ return LatLng{s1.Angle(r.Lat.Lo), s1.Angle(r.Lng.Hi)}
+ case 2:
+ return LatLng{s1.Angle(r.Lat.Hi), s1.Angle(r.Lng.Hi)}
+ case 3:
+ return LatLng{s1.Angle(r.Lat.Hi), s1.Angle(r.Lng.Lo)}
+ default:
+ panic("Invalid vertex index.")
+ }
+}
+
// Center returns the center of the rectangle.
func (r Rect) Center() LatLng {
return LatLng{s1.Angle(r.Lat.Center()) * s1.Radian, s1.Angle(r.Lng.Center()) * s1.Radian}
@@ -119,6 +159,16 @@ func (r Rect) ContainsLatLng(ll LatLng) bool {
return r.Lat.Contains(ll.Lat.Radians()) && r.Lng.Contains(ll.Lng.Radians())
}
+// Return true if and only if the rectangle contains the given other rectangle.
+func (r Rect) ContainsRect(other Rect) bool {
+ return r.Lat.ContainsInterval(other.Lat) && r.Lng.ContainsInterval(other.Lng)
+}
+
+// Return true if this rectangle and the given other rectangle have any points in common.
+func (r Rect) IntersectsRect(other Rect) bool {
+ return r.Lat.Intersects(other.Lat) && r.Lng.Intersects(other.Lng)
+}
+
// AddPoint increases the size of the rectangle to include the given point.
func (r Rect) AddPoint(ll LatLng) Rect {
if !ll.IsValid() {
@@ -151,7 +201,80 @@ func (r Rect) expanded(margin LatLng) Rect {
}
}
+// Return the smallest rectangle containing the union of this rectangle and the given rectangle.
+func (r Rect) Union(other Rect) Rect {
+ return Rect{
+ Lat: r.Lat.Union(other.Lat),
+ Lng: r.Lng.Union(other.Lng),
+ }
+}
+
func (r Rect) String() string { return fmt.Sprintf("[Lo%v, Hi%v]", r.Lo(), r.Hi()) }
+// CapBound returns a bounding spherical cap. This is not guaranteed to be exact.
+func (r Rect) CapBound() Cap {
+ // We consider two possible bounding caps, one whose axis passes
+ // through the center of the lat-long rectangle and one whose axis
+ // is the north or south pole. We return the smaller of the two caps.
+
+ if r.IsEmpty() {
+ return EmptyCap()
+ }
+
+ var poleZ, poleAngle float64
+ if r.Lat.Lo+r.Lat.Hi < 0 {
+ // South pole axis yields smaller cap.
+ poleZ = -1
+ poleAngle = math.Pi/2 + r.Lat.Hi
+ } else {
+ poleZ = 1
+ poleAngle = math.Pi/2 - r.Lat.Lo
+ }
+
+ poleCap := CapFromCenterAngle(PointFromCoordsRaw(0, 0, poleZ), s1.Angle(poleAngle))
+
+ // For bounding rectangles that span 180 degrees or less in longitude, the
+ // maximum cap size is achieved at one of the rectangle vertices. For
+ // rectangles that are larger than 180 degrees, we punt and always return a
+ // bounding cap centered at one of the two poles.
+ lngSpan := r.Lng.Hi - r.Lng.Lo
+ if math.Remainder(lngSpan, 2*math.Pi) >= 0 {
+ if lngSpan < 2*math.Pi {
+ midCap := CapFromCenterAngle(PointFromLatLng(r.Center()), s1.Angle(0))
+ for k := 0; k < 4; k++ {
+ midCap = midCap.AddPoint(PointFromLatLng(r.Vertex(k)))
+ }
+ if midCap.height < poleCap.height {
+ return midCap
+ }
+ }
+ }
+ return poleCap
+}
+
+// RectBound returns a bounding latitude-longitude rectangle that contains
+// the region. The bounds are not guaranteed to be tight.
+func (r Rect) RectBound() Rect {
+ return r
+}
+
+func (r Rect) ContainsCell(c Cell) bool {
+ // A latitude-longitude rectangle contains a cell if and only if it contains
+ // the cell's bounding rectangle. (This is an exact test.)
+ return r.ContainsRect(c.RectBound())
+}
+
+/**
+ * This test is cheap but is NOT exact. Use Intersects() if you want a more
+ * accurate and more expensive test. Note that when this method is used by an
+ * S2RegionCoverer, the accuracy isn't all that important since if a cell may
+ * intersect the region then it is subdivided, and the accuracy of this method
+ * goes up as the cells get smaller.
+ */
+func (r Rect) IntersectsCell(c Cell) bool {
+ // This test is cheap but is NOT exact (see s2latlngrect.h).
+ return r.IntersectsRect(c.RectBound())
+}
+
// BUG(dsymonds): The major differences from the C++ version are:
// - almost everything
diff --git a/s2/region.go b/s2/region.go
index 5edb8f5a..f0763d4e 100644
--- a/s2/region.go
+++ b/s2/region.go
@@ -40,10 +40,11 @@ type Region interface {
}
// Enforce interface satisfaction for a few types once they satisfy the interface.
-/*
var (
_ Region = Cap{}
- _ Region = CellUnion(nil)
+ _ Region = Cell{}
+ // _ Region = CellUnion(nil)
_ Region = Rect{}
+ _ Region = Polyline{}
+ _ Region = &Loop{}
)
-*/
diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go
new file mode 100644
index 00000000..447cda78
--- /dev/null
+++ b/s2/regioncoverer.go
@@ -0,0 +1,418 @@
+package s2
+
+var (
+ DEFAULT_MAX_CELLS int = 8
+
+ FACE_CELLS = []Cell{
+ CellFromCellID(CellIDFromFace(0)),
+ CellFromCellID(CellIDFromFace(1)),
+ CellFromCellID(CellIDFromFace(2)),
+ CellFromCellID(CellIDFromFace(3)),
+ CellFromCellID(CellIDFromFace(4)),
+ CellFromCellID(CellIDFromFace(5)),
+ }
+)
+
+type RegionCoverer struct {
+ minLevel int
+ maxLevel int
+ levelMod int
+ maxCells int
+
+ interiorCovering bool
+
+ candidatesCreatedCounter int
+
+ region Region
+
+ result []CellID
+
+ candidateQueue PriorityQueue
+}
+
+type candidate struct {
+ cell Cell
+ isTerminal bool
+ numChildren int
+ children []*candidate
+}
+
+type queueEntry struct {
+ priority int
+ candidate *candidate
+}
+
+func newQueueEntry(id int, candidate_ *candidate) *queueEntry {
+ return &queueEntry{priority: id, candidate: candidate_}
+}
+
+// A PriorityQueue implements heap.Interface and holds Items.
+type PriorityQueue []*queueEntry
+
+func newPriorityQueue(space int) PriorityQueue {
+ pq := PriorityQueue(make([]*queueEntry, 0, space))
+ return pq
+}
+
+func (pq PriorityQueue) Len() int { return len(pq) }
+
+func (pq PriorityQueue) compare(i, j int) int {
+ if pq[i].priority < pq[j].priority {
+ return 1
+ }
+ if pq[i].priority > pq[j].priority {
+ return -1
+ }
+ return 0
+}
+
+func (pq PriorityQueue) swap(i, j int) {
+ pq[i], pq[j] = pq[j], pq[i]
+}
+
+func (pq *PriorityQueue) Push(x interface{}) {
+ item := x.(*queueEntry)
+ *pq = append(*pq, item)
+
+ pq.up(pq.Len() - 1)
+}
+
+func (pq *PriorityQueue) Pop() interface{} {
+ n := pq.Len() - 1
+ pq.swap(0, n)
+ pq.down(0, n)
+
+ old := *pq
+ n = len(old)
+ item := old[n-1]
+ *pq = old[0 : n-1]
+ return item
+}
+
+func (pq *PriorityQueue) up(j int) {
+ for {
+ i := (j - 1) / 2 // parent
+ if i == j || pq.compare(j, i) >= 0 {
+ break
+ }
+ pq.swap(i, j)
+ j = i
+ }
+}
+
+func (pq *PriorityQueue) down(i, n int) {
+ for {
+ j1 := 2*i + 1
+ if j1 >= n || j1 < 0 { // j1 < 0 after int overflow
+ break
+ }
+ j := j1 // left child
+ if j2 := j1 + 1; j2 < n && pq.compare(j1, j2) > 0 {
+ j = j2 // = 2*i + 2 // right child
+ }
+ if pq.compare(i, j) <= 0 {
+ break
+ }
+ pq.swap(i, j)
+ i = j
+ }
+}
+
+func NewRegionCoverer() *RegionCoverer {
+ return &RegionCoverer{
+ minLevel: 0,
+ maxLevel: MAX_LEVEL,
+ levelMod: 1,
+ maxCells: DEFAULT_MAX_CELLS,
+ candidateQueue: newPriorityQueue(10),
+ }
+}
+
+func (rc *RegionCoverer) SetMinLevel(minLevel int) {
+ rc.minLevel = max(0, min(MAX_LEVEL, minLevel))
+}
+
+func (rc *RegionCoverer) SetMaxLevel(maxLevel int) {
+ rc.maxLevel = max(0, min(MAX_LEVEL, maxLevel))
+}
+
+func (rc *RegionCoverer) SetLevelMod(levelMod int) {
+ rc.levelMod = max(1, min(3, levelMod))
+}
+
+func (rc *RegionCoverer) SetMaxCells(maxCells int) {
+ rc.maxCells = maxCells
+}
+
+func (rc *RegionCoverer) MinLevel() int {
+ return rc.minLevel
+}
+
+func (rc *RegionCoverer) MaxLevel() int {
+ return rc.maxLevel
+}
+
+func (rc *RegionCoverer) LevelMod() int {
+ return rc.levelMod
+}
+
+func (rc *RegionCoverer) MaxCells() int {
+ return rc.maxCells
+}
+
+/**
+ * Computes a list of cell ids that covers the given region and satisfies the
+ * various restrictions specified above.
+ *
+ * @param region The region to cover
+ * @param covering The list filled in by this method
+ */
+func (rc *RegionCoverer) GetCovering(region Region, covering *[]CellID) {
+ // Rather than just returning the raw list of cell ids generated by
+ // GetCoveringInternal(), we construct a cell union and then denormalize it.
+ // This has the effect of replacing four child cells with their parent
+ // whenever this does not violate the covering parameters specified
+ // (min_level, level_mod, etc). This strategy significantly reduces the
+ // number of cells returned in many cases, and it is cheap compared to
+ // computing the covering in the first place.
+ tmp := rc.GetCoveringAsUnion(region)
+ tmp.DeNormalize(rc.minLevel, rc.levelMod, covering)
+}
+
+/**
+ * Computes a list of cell ids that is contained within the given region and
+ * satisfies the various restrictions specified above.
+ *
+ * @param region The region to fill
+ * @param interior The list filled in by this method
+ */
+func (rc *RegionCoverer) GetInteriorCovering(region Region, interior *[]CellID) {
+ tmp := rc.GetInteriorCoveringAsUnion(region)
+ tmp.DeNormalize(rc.minLevel, rc.levelMod, interior)
+}
+
+/**
+ * Return a normalized cell union that covers the given region and satisfies
+ * the restrictions *EXCEPT* for min_level() and level_mod(). These criteria
+ * cannot be satisfied using a cell union because cell unions are
+ * automatically normalized by replacing four child cells with their parent
+ * whenever possible. (Note that the list of cell ids passed to the cell union
+ * constructor does in fact satisfy all the given restrictions.)
+ */
+func (rc *RegionCoverer) GetCoveringAsUnion(region Region) *CellUnion {
+ rc.interiorCovering = false
+ rc.GetCoveringInternal(region)
+ union := CellUnionFromArrayAndSwap(&rc.result)
+ return union
+}
+
+/**
+ * Return a normalized cell union that is contained within the given region
+ * and satisfies the restrictions *EXCEPT* for min_level() and level_mod().
+ */
+func (rc *RegionCoverer) GetInteriorCoveringAsUnion(region Region) *CellUnion {
+ rc.interiorCovering = true
+ rc.GetCoveringInternal(region)
+ union := CellUnionFromArrayAndSwap(&rc.result)
+ return union
+}
+
+/**
+ * Given a connected region and a starting point, return a set of cells at the
+ * given level that cover the region.
+ */
+func GetSimpleCovering(region Region, start Point, level int, output *[]CellID) {
+ floodFill(region, CellFromPoint(start).Id().Parent(level), output)
+}
+
+/** Generates a covering and stores it in result. */
+func (rc *RegionCoverer) GetCoveringInternal(region Region) {
+ if len(rc.result) > 0 || rc.candidateQueue.Len() > 0 {
+ panic("Preconditions not met")
+ }
+
+ rc.region = region
+ rc.candidatesCreatedCounter = 0
+
+ rc.getInitialCandidates()
+ for rc.candidateQueue.Len() > 0 && (!rc.interiorCovering || len(rc.result) < rc.maxCells) {
+ qEntry := rc.candidateQueue.Pop().(*queueEntry)
+ candidate := qEntry.candidate
+ sz := len(rc.result) + candidate.numChildren
+ if !rc.interiorCovering {
+ sz = sz + rc.candidateQueue.Len()
+ }
+ if int(candidate.cell.Level()) < rc.minLevel || candidate.numChildren == 1 || sz <= rc.maxCells {
+ for i := 0; i < candidate.numChildren; i++ {
+ rc.addCandidate(candidate.children[i])
+ }
+
+ } else if rc.interiorCovering {
+ // do nothing
+ } else {
+ candidate.isTerminal = true
+ rc.addCandidate(candidate)
+ }
+
+ }
+
+ rc.candidateQueue = nil
+ rc.region = nil
+}
+
+func (rc *RegionCoverer) getInitialCandidates() {
+ // Optimization: if at least 4 cells are desired (the normal case),
+ // start with a 4-cell covering of the region's bounding cap. This
+ // lets us skip quite a few levels of refinement when the region to
+ // be covered is relatively small.
+ if rc.maxCells >= 4 {
+ // Find the maximum level such that the bounding cap contains at most one
+ // cell vertex at that level.
+ cap := rc.region.CapBound()
+ level := min(S2_PROJECTION.MIN_WIDTH().getMaxLevel(2*cap.Radius().Radians()), min(rc.maxLevel, MAX_LEVEL-1))
+ if rc.levelMod > 1 && level > rc.minLevel {
+ level -= (level - rc.minLevel) % rc.levelMod
+ }
+ // We don't bother trying to optimize the level == 0 case, since more than
+ // four face cells may be required.
+ if level > 0 {
+ // Find the leaf cell containing the cap axis, and determine which
+ // subcell of the parent cell contains it.
+ id := CellIDFromPoint(cap.Center())
+ base := id.VertexNeighbors(level)
+ for i := 0; i < len(base); i++ {
+ rc.addCandidate(rc.newCandidate(CellFromCellID(base[i])))
+ }
+ return
+ }
+ }
+ // Default: start with all six cube faces.
+ for face := 0; face < 6; face++ {
+ rc.addCandidate(rc.newCandidate(FACE_CELLS[face]))
+ }
+}
+
+func (rc *RegionCoverer) addCandidate(candidate_ *candidate) {
+ if candidate_ == nil {
+ return
+ }
+
+ if candidate_.isTerminal {
+ rc.result = append(rc.result, candidate_.cell.Id())
+ return
+ }
+
+ numLevels := rc.levelMod
+ if int(candidate_.cell.Level()) < rc.minLevel {
+ numLevels = 1
+ }
+ numTerminals := rc.expandChildren(candidate_, candidate_.cell, numLevels)
+
+ if candidate_.numChildren == 0 {
+ // do nothing
+ } else if !rc.interiorCovering && numTerminals == 1<
+ * 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.
+ */
+func RobustCCW(a, b, c Point) int {
+ return RobustCCWWithCross(a, b, c, Point{a.Cross(b.Vector)})
+}
+
+/**
+ * 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
+ */
+func RobustCCWWithCross(a, b, c, aCrossB Point) int {
+ // 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).
+ kMinAbsValue := 1.6e-15 // 2 * 14 * 2**-54
+
+ det := aCrossB.Dot(c.Vector)
+
+ // 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.
+ */
+func ExpensiveCCW(a, b, c Point) int {
+ // 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.
+ var sab float64 = 1
+ var sbc float64 = 1
+ var sca float64 = 1
+ if a.Dot(b.Vector) > 0 {
+ sab = -1
+ }
+ if b.Dot(c.Vector) > 0 {
+ sbc = -1
+ }
+ if c.Dot(a.Vector) > 0 {
+ sca = -1
+ }
+ vab := a.Add(b.Mul(sab))
+ vbc := b.Add(c.Mul(sbc))
+ vca := c.Add(a.Mul(sca))
+ dab := vab.Norm2()
+ dbc := vbc.Norm2()
+ 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.
+ var sign float64
+ 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 = vab.Cross(vca).Dot(a.Vector) * sab // BC is longest
+ // edge
+ } else {
+ sign = vca.Cross(vbc).Dot(c.Vector) * sca // AB is longest
+ // edge
+ }
+ } else {
+ if dab < dca || (dab == dca && b.LessThan(c)) {
+ sign = vbc.Cross(vab).Dot(b.Vector) * sbc // CA is longest
+ // edge
+ } else {
+ sign = vca.Cross(vbc).Dot(c.Vector) * 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.
+
+ ccw := PlanarOrderedCCW(r2.Vector{a.Y, a.Z}, r2.Vector{b.Y, b.Z}, r2.Vector{c.Y, c.Z})
+ if ccw == 0 {
+ ccw = PlanarOrderedCCW(r2.Vector{a.Z, a.X}, r2.Vector{b.Z, b.X}, r2.Vector{c.Z, c.X})
+ if ccw == 0 {
+ ccw = PlanarOrderedCCW(
+ r2.Vector{a.X, a.Y}, r2.Vector{b.X, b.Y}, r2.Vector{c.X, c.Y})
+ // assert (ccw != 0);
+ }
+ }
+ return ccw
+}
+
+func PlanarCCW(a, b r2.Vector) int {
+ // Return +1 if the edge AB is CCW around the origin, etc.
+ var sab float64 = 1
+ if a.Dot(b) > 0 {
+ sab = -1
+ }
+ vab := a.Add(b.Mul(sab))
+ da := a.Norm2()
+ db := b.Norm2()
+ var sign float64
+ if da < db || (da == db && a.LessThan(b)) {
+ sign = a.Cross(vab) * sab
+ } else {
+ sign = vab.Cross(b)
+ }
+ if sign > 0 {
+ return 1
+ }
+ if sign < 0 {
+ return -1
+ }
+ return 0
+}
+
+func PlanarOrderedCCW(a, b, c r2.Vector) 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:
+ *
+ *
+ */
+func OrderedCCW(a, b, c, o Point) bool {
+ // 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.
+
+ 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
+}
+
+// Defines an area or a length cell metric.
+type Metric struct {
+ deriv float64
+ dim uint
+}
+
+// Defines a cell metric of the given dimension (1 == length, 2 == area).
+func NewMetric(dim uint, deriv float64) Metric {
+ return Metric{deriv, 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.
+func (m Metric) Deriv() float64 { return m.deriv }
+
+// Return the value of a metric for cells at the given level.
+func (m Metric) GetValue(level int) float64 {
+ return m.deriv * math.Pow(2, float64(int(m.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.
+ */
+func (m Metric) getClosestLevel(value float64) int {
+ if m.dim == 1 {
+ return m.getMinLevel(math.Sqrt2 * value)
+ }
+ return m.getMinLevel(2 * 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.
+ */
+func (m Metric) getMinLevel(value float64) int {
+ if value <= 0 {
+ return MAX_LEVEL
+ }
+
+ // This code is equivalent to computing a floating-point "level"
+ // value and rounding up.
+ exponent := exp(value / (float64(int(1)<