From 99492ef7b4fa38c50d2ead421a677b49399445ab Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Thu, 9 Apr 2015 18:17:40 -0700 Subject: [PATCH 01/23] Added initial regioncoverer which for now always returns a fixed cube face --- s2/regioncoverer.go | 20 ++++++++++++++++++++ s2/regioncoverer_test.go | 19 +++++++++++++++++++ 2 files changed, 39 insertions(+) create mode 100644 s2/regioncoverer.go create mode 100644 s2/regioncoverer_test.go diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go new file mode 100644 index 00000000..ea8f02f1 --- /dev/null +++ b/s2/regioncoverer.go @@ -0,0 +1,20 @@ +package s2 + +type RegionCoverer struct { + maxCells int +} + +func NewRegionCoverer() *RegionCoverer { + return &RegionCoverer{} +} + +func (rc *RegionCoverer) SetMaxCells(maxCells int) { + rc.maxCells = maxCells + +} + +func (rc *RegionCoverer) GetCovering(region Region) *CellUnion { + return &CellUnion{ + 0x9000000000000000, + } +} diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go new file mode 100644 index 00000000..216a0951 --- /dev/null +++ b/s2/regioncoverer_test.go @@ -0,0 +1,19 @@ +package s2 + +import ( + "testing" +) + +func TestCovering(t *testing.T) { + coverer := NewRegionCoverer() + coverer.SetMaxCells(8) + cells := coverer.GetCovering(nil) + for _, cell := range *cells { + if cell.Level() != 1 { + t.Errorf("Level not as expected %s", cell.ToToken()) + } + if cell.Face() != 4 { + t.Error("Face not as expected") + } + } +} From 6a7bb06ffbc8772a127509e2810f144b7d18ec1b Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Thu, 9 Apr 2015 18:17:56 -0700 Subject: [PATCH 02/23] Added missing accessors --- s2/cell.go | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/s2/cell.go b/s2/cell.go index 49a26fda..7c7717a3 100644 --- a/s2/cell.go +++ b/s2/cell.go @@ -53,6 +53,22 @@ func CellFromLatLng(ll LatLng) Cell { return CellFromCellID(CellIDFromLatLng(ll)) } +func (c Cell) Id() CellID { + return c.id +} + +func (c Cell) Face() int8 { + return c.face +} + +func (c Cell) Level() int8 { + return c.level +} + +func (c Cell) Orientation() int8 { + return c.orientation +} + // IsLeaf returns whether this Cell is a leaf or not. func (c Cell) IsLeaf() bool { return c.level == maxLevel From 3a8834c7fb52a4e3704e28bf4965da08c817998e Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Thu, 23 Apr 2015 19:00:29 -0700 Subject: [PATCH 03/23] Work in progress on RegionCoverer --- s2/cellid.go | 2 + s2/cellunion.go | 36 ++++++++ s2/regioncoverer.go | 215 +++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 251 insertions(+), 2 deletions(-) diff --git a/s2/cellid.go b/s2/cellid.go index 0d1a2aee..ea29c8a9 100644 --- a/s2/cellid.go +++ b/s2/cellid.go @@ -38,6 +38,8 @@ type CellID uint64 // TODO(dsymonds): Some of these constants should probably be exported. const ( + MAX_LEVEL = 30 + faceBits = 3 numFaces = 6 maxLevel = 30 diff --git a/s2/cellunion.go b/s2/cellunion.go index 25b623b6..275cbe99 100644 --- a/s2/cellunion.go +++ b/s2/cellunion.go @@ -25,6 +25,42 @@ import "sort" // nor the four sibling CellIDs that are children of a single higher level CellID. type CellUnion []CellID +func min(a, b int) int { + if a > b { + return b + } + return a +} + +func max(a, b int) int { + if a < b { + return b + } + return a +} + +func (cu *CellUnion) DeNormalize(minLevel, levelMod int, output *CellUnion) { + *output = make([]CellID, len(*cu)) + for _, ci := range *cu { + level := ci.Level() + newLevel := max(minLevel, level) + if levelMod > 1 { + // Round up so that (new_level - min_level) is a multiple of level_mod. + // (Note that S2CellId::kMaxLevel is a multiple of 1, 2, and 3.) + newLevel += (MAX_LEVEL - (newLevel - minLevel)) % levelMod + newLevel = min(MAX_LEVEL, newLevel) + } + if newLevel == level { + *output = append(*output, ci) + } else { + end := ci.ChildEndAtLevel(newLevel) + for id := ci.ChildBeginAtLevel(newLevel); id != end; id = id.Next() { + *output = append(*output, id) + } + } + } +} + // Normalize normalizes the CellUnion. func (cu *CellUnion) Normalize() { sort.Sort(byID(*cu)) diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index ea8f02f1..4ce7559c 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -1,20 +1,231 @@ 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 []queueEntry +} + +type candidate struct { + cell Cell + isTerminal bool + numChildren int + children []*candidate +} + +type queueEntry struct { + id int + candidate *candidate +} + +func newQueueEntry(id int, candidate *candidate) *queueEntry { + return &queueEntry{id, candidate} } func NewRegionCoverer() *RegionCoverer { - return &RegionCoverer{} + return &RegionCoverer{ + minLevel: 0, + maxLevel: MAX_LEVEL, + levelMod: 1, + maxCells: DEFAULT_MAX_CELLS, + } +} + +func (rc *RegionCoverer) SetMinLevel(minLevel int) { + rc.minLevel = minLevel +} + +func (rc *RegionCoverer) SetMaxLevel(maxLevel int) { + rc.maxLevel = maxLevel +} + +func (rc *RegionCoverer) SetLevelMod(levelMod int) { + rc.levelMod = 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 } func (rc *RegionCoverer) GetCovering(region Region) *CellUnion { - return &CellUnion{ + rc.interiorCovering = false + rc.GetCoveringInternal(region) + return (*CellUnion)(&rc.result) +} + +func (rc *RegionCoverer) GetCoveringInternal(region Region) { + if len(rc.result) > 0 || len(rc.candidateQueue) > 0 { + panic("Preconditions not met") + } + + rc.region = region + rc.candidatesCreatedCounter = 0 + + rc.getInitialCandidates() + for len(rc.candidateQueue) > 0 && (!rc.interiorCovering || len(rc.result) < rc.maxCells) { + candidate := rc.candidateQueue[0].candidate // TODO: rc.candidateQueue.poll().candidate + sz := len(rc.result) + candidate.numChildren + if rc.interiorCovering { + sz = sz + len(rc.candidateQueue) + } + 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() { + // TODO + rc.result = []CellID{ 0x9000000000000000, } + +} + +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<= rc.minLevel { + // Optimization: add the parent cell rather than all of its children. + // We can't do this for interior coverings, since the children just + // intersect the region, but may not be contained by it - we need to + // subdivide them further. + candidate.isTerminal = true + rc.addCandidate(candidate) + } else { + // We negate the priority so that smaller absolute priorities are returned + // first. The heuristic is designed to refine the largest cells first, + // since those are where we have the largest potential gain. Among cells + // at the same level, we prefer the cells with the smallest number of + // intersecting children. Finally, we prefer cells that have the smallest + // number of children that cannot be refined any further. + // priority := -((((int(candidate.cell.Level()) << rc.maxChildrenShift()) + candidate.numChildren) << rc.maxChildrenShift()) + numTerminals) + // TODO: candidateQueue.add(new QueueEntry(priority, candidate)); + } +} + +func (rc *RegionCoverer) maxChildrenShift() uint { + return 2 * uint(rc.levelMod) +} + +func (rc *RegionCoverer) expandChildren(candidate *candidate, cell Cell, numLevels int) int { + numLevels-- + childCellIds := cell.Id().Children() + numTerminals := 0 + for i := 0; i < 4; i++ { + childCell := CellFromCellID(childCellIds[i]) + if numLevels > 0 { + if rc.region.IntersectsCell(childCell) { + numTerminals = numTerminals + rc.expandChildren(candidate, childCell, numLevels) + } + continue + } + child := rc.newCandidate(childCell) + if child != nil { + candidate.children[candidate.numChildren] = child + candidate.numChildren++ + if child.isTerminal { + numTerminals++ + } + } + } + return numTerminals +} + +func (rc *RegionCoverer) newCandidate(cell Cell) *candidate { + if rc.region.IntersectsCell(cell) { + return nil + } + isTerminal := false + if int(cell.Level()) >= rc.minLevel { + if rc.interiorCovering { + if rc.region.ContainsCell(cell) { + isTerminal = true + } else if int(cell.Level())+rc.levelMod > rc.maxLevel { + return nil + } + } else { + if int(cell.Level())+rc.levelMod > rc.maxLevel || rc.region.ContainsCell(cell) { + isTerminal = true + } + } + + } + candidate_ := &candidate{ + cell: cell, + isTerminal: isTerminal, + } + if !isTerminal { + candidate_.children = make([]*candidate, 1< Date: Thu, 23 Apr 2015 23:09:08 -0700 Subject: [PATCH 04/23] Added more RegionCoverer functionality --- s2/cap.go | 5 ++ s2/cell.go | 41 +++++++++++++-- s2/cellid.go | 6 +-- s2/cellid_test.go | 2 +- s2/cellunion.go | 14 ------ s2/projections.go | 53 ++++++++++++++++++++ s2/regioncoverer.go | 104 ++++++++++++++++++++++++++++++++------ s2/regioncoverer_test.go | 3 +- s2/s2.go | 105 +++++++++++++++++++++++++++++++++++++++ 9 files changed, 296 insertions(+), 37 deletions(-) create mode 100644 s2/projections.go create mode 100644 s2/s2.go diff --git a/s2/cap.go b/s2/cap.go index 54d932a8..344c121e 100644 --- a/s2/cap.go +++ b/s2/cap.go @@ -118,6 +118,11 @@ func (c Cap) IsFull() bool { return c.height == fullHeight } +// Center returns the cap's center or axis. +func (c Cap) Center() Point { + return c.center +} + // Radius returns the cap's radius. func (c Cap) Radius() s1.Angle { if c.IsEmpty() { diff --git a/s2/cell.go b/s2/cell.go index 7c7717a3..2a5d3dab 100644 --- a/s2/cell.go +++ b/s2/cell.go @@ -45,7 +45,7 @@ func CellFromCellID(id CellID) Cell { // CellFromPoint constructs a cell for the given Point. func CellFromPoint(p Point) Cell { - return CellFromCellID(cellIDFromPoint(p)) + return CellFromCellID(CellIDFromPoint(p)) } // CellFromLatLng constructs a cell for the given LatLng. @@ -106,5 +106,40 @@ func (c Cell) ExactArea() float64 { return PointArea(v0, v1, v2) + PointArea(v0, v2, v3) } -// TODO(roberts, or $SOMEONE): Differences from C++, almost everything else still. -// Implement the accessor methods on the internal fields. +// CapBound returns a bounding spherical cap. This is not guaranteed to be exact. +func (c Cell) CapBound() Cap { + // Use the cell center in (u,v)-space as the cap axis. This vector is + // very close to GetCenter() and faster to compute. Neither one of these + // vectors yields the bounding cap with minimal surface area, but they + // are both pretty close. + // + // It's possible to show that the two vertices that are furthest from + // the (u,v)-origin never determine the maximum cap size (this is a + // possible future optimization). + u := c.uv.Center().X + v := c.uv.Center().Y + cap := CapFromCenterHeight(Point{faceUVToXYZ(int(c.face), u, v).Normalize()}, 0) + for k := 0; k < 4; k++ { + cap = cap.AddPoint(c.Vertex(k)) + } + return cap +} + +// RectBound returns a bounding latitude-longitude rectangle that contains +// the region. The bounds are not guaranteed to be tight. +func (c Cell) RectBound() Rect { + return EmptyRect() +} + +// ContainsCell reports whether the region completely contains the given region. +// It returns false if containment could not be determined. +func (c Cell) ContainsCell(other Cell) bool { + return c.Id().Contains(other.Id()) +} + +// 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 (c Cell) IntersectsCell(other Cell) bool { + return c.Id().Intersects(other.Id()) +} diff --git a/s2/cellid.go b/s2/cellid.go index ea29c8a9..011450bd 100644 --- a/s2/cellid.go +++ b/s2/cellid.go @@ -63,7 +63,7 @@ func CellIDFromFace(face int) CellID { // CellIDFromLatLng returns the leaf cell containing ll. func CellIDFromLatLng(ll LatLng) CellID { - return cellIDFromPoint(PointFromLatLng(ll)) + return CellIDFromPoint(PointFromLatLng(ll)) } // CellIDFromToken returns a cell given a hex-encoded string of its uint64 ID. @@ -449,8 +449,8 @@ func stToIJ(s float64) int { return clamp(int(math.Floor(maxSize*s)), 0, maxSize-1) } -// cellIDFromPoint returns the leaf cell containing point p. -func cellIDFromPoint(p Point) CellID { +// CellIDFromPoint returns the leaf cell containing point p. +func CellIDFromPoint(p Point) CellID { f, u, v := xyzToFaceUV(r3.Vector{p.X, p.Y, p.Z}) i := stToIJ(uvToST(u)) j := stToIJ(uvToST(v)) diff --git a/s2/cellid_test.go b/s2/cellid_test.go index bbed769f..53c8c0b2 100644 --- a/s2/cellid_test.go +++ b/s2/cellid_test.go @@ -210,7 +210,7 @@ func (v byCellID) Less(i, j int) bool { return uint64(v[i]) < uint64(v[j]) } func TestVertexNeighbors(t *testing.T) { // Check the vertex neighbors of the center of face 2 at level 5. - id := cellIDFromPoint(PointFromCoords(0, 0, 1)) + id := CellIDFromPoint(PointFromCoords(0, 0, 1)) neighbors := id.VertexNeighbors(5) sort.Sort(byCellID(neighbors)) diff --git a/s2/cellunion.go b/s2/cellunion.go index 275cbe99..c2fc0f75 100644 --- a/s2/cellunion.go +++ b/s2/cellunion.go @@ -25,20 +25,6 @@ import "sort" // nor the four sibling CellIDs that are children of a single higher level CellID. type CellUnion []CellID -func min(a, b int) int { - if a > b { - return b - } - return a -} - -func max(a, b int) int { - if a < b { - return b - } - return a -} - func (cu *CellUnion) DeNormalize(minLevel, levelMod int, output *CellUnion) { *output = make([]CellID, len(*cu)) for _, ci := range *cu { diff --git a/s2/projections.go b/s2/projections.go new file mode 100644 index 00000000..db779a18 --- /dev/null +++ b/s2/projections.go @@ -0,0 +1,53 @@ +package s2 + +import ( + "math" +) + +// type Projections int + +// const ( +// S2_LINEAR_PROJECTION Projections = iota +// S2_TAN_PROJECTION +// S2_QUADRATIC_PROJECTION +// ) + +// const S2_PROJECTION Projections = S2_QUADRATIC_PROJECTION + +var ( + S2_PROJECTION Projections = S2_QUADRATIC_PROJECTION{} +) + +type Projections interface { + MAX_ANGLE_SPAN() Metric + MIN_WIDTH() Metric + MAX_WIDTH() Metric + AVG_WIDTH() Metric +} + +type S2_QUADRATIC_PROJECTION struct{} + +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, + } +} diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index 4ce7559c..0f03fa37 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -1,5 +1,9 @@ package s2 +import ( + "container/heap" +) + var ( DEFAULT_MAX_CELLS int = 8 @@ -27,7 +31,7 @@ type RegionCoverer struct { result []CellID - candidateQueue []queueEntry + candidateQueue PriorityQueue } type candidate struct { @@ -38,20 +42,66 @@ type candidate struct { } type queueEntry struct { - id int + priority int candidate *candidate + // The index is needed by update and is maintained by the heap.Interface methods. + index int // The index of the item in the heap. } func newQueueEntry(id int, candidate *candidate) *queueEntry { - return &queueEntry{id, candidate} + return &queueEntry{priority: id, candidate: candidate} +} + +// A PriorityQueue implements heap.Interface and holds Items. +type PriorityQueue []*queueEntry + +func newPriorityQueue(space int) PriorityQueue { + return make([]*queueEntry, 0, space) +} + +func (pq PriorityQueue) Len() int { return len(pq) } + +func (pq PriorityQueue) Less(i, j int) bool { + // We want Pop to give us the highest, not lowest, priority so we use greater than here. + return pq[i].priority > pq[j].priority +} + +func (pq PriorityQueue) Swap(i, j int) { + pq[i], pq[j] = pq[j], pq[i] + pq[i].index = i + pq[j].index = j +} + +func (pq *PriorityQueue) Push(x interface{}) { + n := len(*pq) + item := x.(*queueEntry) + item.index = n + *pq = append(*pq, item) +} + +func (pq *PriorityQueue) Pop() interface{} { + old := *pq + n := len(old) + item := old[n-1] + item.index = -1 // for safety + *pq = old[0 : n-1] + return item +} + +// update modifies the priority and value of an Item in the queue. +func (pq *PriorityQueue) update(item *queueEntry, candidate *candidate, priority int) { + item.candidate = candidate + item.priority = priority + heap.Fix(pq, item.index) } func NewRegionCoverer() *RegionCoverer { return &RegionCoverer{ - minLevel: 0, - maxLevel: MAX_LEVEL, - levelMod: 1, - maxCells: DEFAULT_MAX_CELLS, + minLevel: 0, + maxLevel: MAX_LEVEL, + levelMod: 1, + maxCells: DEFAULT_MAX_CELLS, + candidateQueue: newPriorityQueue(10), } } @@ -94,7 +144,7 @@ func (rc *RegionCoverer) GetCovering(region Region) *CellUnion { } func (rc *RegionCoverer) GetCoveringInternal(region Region) { - if len(rc.result) > 0 || len(rc.candidateQueue) > 0 { + if len(rc.result) > 0 || rc.candidateQueue.Len() > 0 { panic("Preconditions not met") } @@ -103,7 +153,7 @@ func (rc *RegionCoverer) GetCoveringInternal(region Region) { rc.getInitialCandidates() for len(rc.candidateQueue) > 0 && (!rc.interiorCovering || len(rc.result) < rc.maxCells) { - candidate := rc.candidateQueue[0].candidate // TODO: rc.candidateQueue.poll().candidate + candidate := rc.candidateQueue.Pop().(*queueEntry).candidate sz := len(rc.result) + candidate.numChildren if rc.interiorCovering { sz = sz + len(rc.candidateQueue) @@ -127,11 +177,35 @@ func (rc *RegionCoverer) GetCoveringInternal(region Region) { } func (rc *RegionCoverer) getInitialCandidates() { - // TODO - rc.result = []CellID{ - 0x9000000000000000, + // 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) { @@ -167,8 +241,8 @@ func (rc *RegionCoverer) addCandidate(candidate *candidate) { // at the same level, we prefer the cells with the smallest number of // intersecting children. Finally, we prefer cells that have the smallest // number of children that cannot be refined any further. - // priority := -((((int(candidate.cell.Level()) << rc.maxChildrenShift()) + candidate.numChildren) << rc.maxChildrenShift()) + numTerminals) - // TODO: candidateQueue.add(new QueueEntry(priority, candidate)); + priority := -((((int(candidate.cell.Level()) << rc.maxChildrenShift()) + candidate.numChildren) << rc.maxChildrenShift()) + numTerminals) + rc.candidateQueue.Push(newQueueEntry(priority, candidate)) } } diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index 216a0951..9e6de8a1 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -7,7 +7,8 @@ import ( func TestCovering(t *testing.T) { coverer := NewRegionCoverer() coverer.SetMaxCells(8) - cells := coverer.GetCovering(nil) + region := CellFromCellID(CellIDFromToken("80c297c53c")) + cells := coverer.GetCovering(region) for _, cell := range *cells { if cell.Level() != 1 { t.Errorf("Level not as expected %s", cell.ToToken()) diff --git a/s2/s2.go b/s2/s2.go new file mode 100644 index 00000000..9dbd6c79 --- /dev/null +++ b/s2/s2.go @@ -0,0 +1,105 @@ +package s2 + +import ( + "math" +) + +// Number of bits in the mantissa of a double. +const EXPONENT_SHIFT uint = 52 + +// Mask to extract the exponent from a double. +const EXPONENT_MASK uint64 = 0x7ff0000000000000 + +func min(a, b int) int { + if a > b { + return b + } + return a +} + +func max(a, b int) int { + if a < b { + return b + } + return a +} + +func exp(v float64) int { + if v == 0 { + return 0 + } + bits := math.Float64bits(v) + return (int)((EXPONENT_MASK&bits)>>EXPONENT_SHIFT) - 1022 +} + +// 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 math.Pow(m.deriv, 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 { + return m.getMinLevel(math.Sqrt2 * value) +} + +/** + * Return the minimum level such that the metric is at most the given value, + * or S2CellId::kMaxLevel if there is no such level. For example, + * S2::kMaxDiag.GetMinLevel(0.1) returns the minimum level such that all + * cell diagonal lengths are 0.1 or smaller. The return value is always a + * valid level. + */ +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)<>(m.dim-1)))) + // assert (level == S2CellId.MAX_LEVEL || getValue(level) <= value); + // assert (level == 0 || getValue(level - 1) > value); + return level +} + +/** + * Return the maximum level such that the metric is at least the given + * value, or zero if there is no such level. For example, + * S2.kMinWidth.GetMaxLevel(0.1) returns the maximum level such that all + * cells have a minimum width of 0.1 or larger. The return value is always a + * valid level. + */ +func (m Metric) getMaxLevel(value float64) int { + if value <= 0 { + return MAX_LEVEL + } + + // This code is equivalent to computing a floating-point "level" + // value and rounding down. + exponent := exp(float64(int(1)<>(m.dim-1)))) + // assert (level == 0 || getValue(level) >= value); + // assert (level == S2CellId.MAX_LEVEL || getValue(level + 1) < value); + return level +} From 17eb84dc06043db4cbb6ec33a98ca534fae698cc Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Sat, 25 Apr 2015 14:07:50 -0700 Subject: [PATCH 05/23] More work in progress --- s2/cellunion.go | 11 ++++++++- s2/regioncoverer.go | 48 ++++++++++++++++++++++++++++++++-------- s2/regioncoverer_test.go | 9 +++++--- 3 files changed, 55 insertions(+), 13 deletions(-) diff --git a/s2/cellunion.go b/s2/cellunion.go index c2fc0f75..935ae609 100644 --- a/s2/cellunion.go +++ b/s2/cellunion.go @@ -25,7 +25,16 @@ import "sort" // nor the four sibling CellIDs that are children of a single higher level CellID. type CellUnion []CellID -func (cu *CellUnion) DeNormalize(minLevel, levelMod int, output *CellUnion) { +func CellUnionFromArrayAndSwap(ids *[]CellID) *CellUnion { + union := &CellUnion{} + *union = append(*union, *ids...) + union.Normalize() + *ids = []CellID{} + return union +} + + +func (cu *CellUnion) DeNormalize(minLevel, levelMod int, output *[]CellID) { *output = make([]CellID, len(*cu)) for _, ci := range *cu { level := ci.Level() diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index 0f03fa37..bac8a00f 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -63,7 +63,7 @@ func (pq PriorityQueue) Len() int { return len(pq) } func (pq PriorityQueue) Less(i, j int) bool { // We want Pop to give us the highest, not lowest, priority so we use greater than here. - return pq[i].priority > pq[j].priority + return pq[i].priority < pq[j].priority } func (pq PriorityQueue) Swap(i, j int) { @@ -106,15 +106,15 @@ func NewRegionCoverer() *RegionCoverer { } func (rc *RegionCoverer) SetMinLevel(minLevel int) { - rc.minLevel = minLevel + rc.minLevel = max(0, min(MAX_LEVEL, minLevel)) } func (rc *RegionCoverer) SetMaxLevel(maxLevel int) { - rc.maxLevel = maxLevel + rc.maxLevel = max(0, min(MAX_LEVEL, maxLevel)) } func (rc *RegionCoverer) SetLevelMod(levelMod int) { - rc.levelMod = levelMod + rc.levelMod = max(1, min(3, levelMod)) } func (rc *RegionCoverer) SetMaxCells(maxCells int) { @@ -137,12 +137,42 @@ func (rc *RegionCoverer) MaxCells() int { return rc.maxCells } -func (rc *RegionCoverer) GetCovering(region Region) *CellUnion { +/** + * 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) +} + + +/** + * 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) - return (*CellUnion)(&rc.result) + union := CellUnionFromArrayAndSwap(&rc.result) + return union } + /** 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") @@ -152,11 +182,11 @@ func (rc *RegionCoverer) GetCoveringInternal(region Region) { rc.candidatesCreatedCounter = 0 rc.getInitialCandidates() - for len(rc.candidateQueue) > 0 && (!rc.interiorCovering || len(rc.result) < rc.maxCells) { + for rc.candidateQueue.Len() > 0 && (!rc.interiorCovering || len(rc.result) < rc.maxCells) { candidate := rc.candidateQueue.Pop().(*queueEntry).candidate sz := len(rc.result) + candidate.numChildren - if rc.interiorCovering { - sz = sz + len(rc.candidateQueue) + 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++ { diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index 9e6de8a1..c44fbf20 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -7,9 +7,12 @@ import ( func TestCovering(t *testing.T) { coverer := NewRegionCoverer() coverer.SetMaxCells(8) - region := CellFromCellID(CellIDFromToken("80c297c53c")) - cells := coverer.GetCovering(region) - for _, cell := range *cells { + region := CellFromCellID(CellIDFromToken("80c297c53")) + cells := []CellID{} + // cells := coverer.GetCoveringAsUnion(region) + coverer.GetCovering(region, &cells) + for _, cell := range cells { + t.Errorf("cell: %s", cell.ToToken()) if cell.Level() != 1 { t.Errorf("Level not as expected %s", cell.ToToken()) } From 862dbdfa584a306b3a5d9413410438430e35dea2 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Sat, 25 Apr 2015 23:26:32 -0700 Subject: [PATCH 06/23] Fixed two bugs, added Cap Region interface --- s2/cap.go | 83 ++++++++++++++++++++++++++++++++++++++++ s2/cell.go | 24 ++++++++++-- s2/cellunion.go | 3 +- s2/regioncoverer.go | 51 ++++++++++++------------ s2/regioncoverer_test.go | 14 +++---- 5 files changed, 134 insertions(+), 41 deletions(-) diff --git a/s2/cap.go b/s2/cap.go index 344c121e..3018f00e 100644 --- a/s2/cap.go +++ b/s2/cap.go @@ -317,6 +317,89 @@ func radiusToHeight(r s1.Angle) float64 { } +// ContainsCell reports whether the region completely contains the given region. +// It returns false if containment could not be determined. +func (c Cap) ContainsCell(cell Cell) bool { + vertices := make([]Point, 4) + for k := 0; k < 4; k++ { + vertices[k] = cell.Vertex(k) + if !c.ContainsPoint(vertices[k]) { + return false + } + } + return !c.Complement().intersects(cell, vertices) +} + +// 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 (c Cap) IntersectsCell(cell Cell) bool { + vertices := make([]Point, 4) + for k := 0; k < 4; k++ { + vertices[k] = cell.Vertex(k) + if c.ContainsPoint(vertices[k]) { + return true + } + } + return c.intersects(cell, vertices) +} + +/** + * Return true if the cap intersects 'cell', given that the cap vertices have + * alrady been checked. + */ +func (c Cap) intersects(cell Cell, vertices []Point) bool { + // Return true if this cap intersects any point of 'cell' excluding its + // vertices (which are assumed to already have been checked). + + // If the cap is a hemisphere or larger, the cell and the complement of the + // cap are both convex. Therefore since no vertex of the cell is contained, + // no other interior point of the cell is contained either. + if c.height >= 1 { + return false + } + + // We need to check for empty caps due to the axis check just below. + if c.IsEmpty() { + return false + } + + // Optimization: return true if the cell contains the cap axis. (This + // allows half of the edge checks below to be skipped.) + if cell.ContainsPoint(c.center) { + return true + } + + // At this point we know that the cell does not contain the cap axis, + // and the cap does not contain any cell vertex. The only way that they + // can intersect is if the cap intersects the interior of some edge. + + sin2Angle := c.height * (2 - c.height) // sin^2(capAngle) + for k := 0; k < 4; k++ { + edge := cell.EdgeRaw(k) + dot := c.center.Dot(edge.Vector) + if dot > 0 { + // The axis is in the interior half-space defined by the edge. We don't + // need to consider these edges, since if the cap intersects this edge + // then it also intersects the edge on the opposite side of the cell + // (because we know the axis is not contained with the cell). + continue + } + // The Norm2() factor is necessary because "edge" is not normalized. + if dot*dot > sin2Angle*edge.Norm2() { + return false // Entire cap is on the exterior side of this edge. + } + // Otherwise, the great circle containing this edge intersects + // the interior of the cap. We just need to check whether the point + // of closest approach occurs between the two edge endpoints. + dir := edge.Cross(c.center.Vector) + if dir.Dot(vertices[k].Vector) < 0 && dir.Dot(vertices[(k+1)&3].Vector) > 0 { + return true + } + } + return false +} + // TODO(roberts): Differences from C++ // Intersects(S2Cell), Contains(S2Cell), MayIntersect(S2Cell) // Centroid, Union diff --git a/s2/cell.go b/s2/cell.go index 2a5d3dab..0a218c9d 100644 --- a/s2/cell.go +++ b/s2/cell.go @@ -88,15 +88,19 @@ func (c Cell) Vertex(k int) Point { // Edge returns the inward-facing normal of the great circle passing through // the CCW ordered edge from vertex k to vertex k+1 (mod 4). func (c Cell) Edge(k int) Point { + return Point{c.EdgeRaw(k).Normalize()} +} + +func (c Cell) EdgeRaw(k int) Point { switch k { case 0: - return Point{vNorm(int(c.face), c.uv.Y.Lo).Normalize()} // Bottom + return Point{vNorm(int(c.face), c.uv.Y.Lo)} // Bottom case 1: - return Point{uNorm(int(c.face), c.uv.X.Hi).Normalize()} // Right + return Point{uNorm(int(c.face), c.uv.X.Hi)} // Right case 2: - return Point{vNorm(int(c.face), c.uv.Y.Hi).Mul(-1.0).Normalize()} // Top + return Point{vNorm(int(c.face), c.uv.Y.Hi).Mul(-1.0)} // Top default: - return Point{uNorm(int(c.face), c.uv.X.Lo).Mul(-1.0).Normalize()} // Left + return Point{uNorm(int(c.face), c.uv.X.Lo).Mul(-1.0)} // Left } } @@ -125,9 +129,21 @@ func (c Cell) CapBound() Cap { return cap } +func (c Cell) ContainsPoint(point Point) bool { + // We can't just call XYZtoFaceUV, because for points that lie on the + // boundary between two faces (i.e. u or v is +1/-1) we need to return + // true for both adjacent cells. + u, v, ok := faceXYZToUV(int(c.face), point) + if !ok { + return false + } + return u >= c.uv.X.Lo && u <= c.uv.X.Hi && v >= c.uv.Y.Lo && v <= c.uv.Y.Hi +} + // RectBound returns a bounding latitude-longitude rectangle that contains // the region. The bounds are not guaranteed to be tight. func (c Cell) RectBound() Rect { + // TODO: Implement return EmptyRect() } diff --git a/s2/cellunion.go b/s2/cellunion.go index 935ae609..f770ce1b 100644 --- a/s2/cellunion.go +++ b/s2/cellunion.go @@ -33,9 +33,8 @@ func CellUnionFromArrayAndSwap(ids *[]CellID) *CellUnion { return union } - func (cu *CellUnion) DeNormalize(minLevel, levelMod int, output *[]CellID) { - *output = make([]CellID, len(*cu)) + *output = make([]CellID, 0, len(*cu)) for _, ci := range *cu { level := ci.Level() newLevel := max(minLevel, level) diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index bac8a00f..93ce52f3 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -138,33 +138,32 @@ func (rc *RegionCoverer) MaxCells() int { } /** - * 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 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) } - /** - * 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.) - */ + * 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) @@ -172,7 +171,7 @@ func (rc *RegionCoverer) GetCoveringAsUnion(region Region) *CellUnion { return union } - /** Generates a covering and stores it in result. */ +/** 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") @@ -305,7 +304,7 @@ func (rc *RegionCoverer) expandChildren(candidate *candidate, cell Cell, numLeve } func (rc *RegionCoverer) newCandidate(cell Cell) *candidate { - if rc.region.IntersectsCell(cell) { + if !rc.region.IntersectsCell(cell) { return nil } isTerminal := false diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index c44fbf20..35d7ad39 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -7,17 +7,13 @@ import ( func TestCovering(t *testing.T) { coverer := NewRegionCoverer() coverer.SetMaxCells(8) - region := CellFromCellID(CellIDFromToken("80c297c53")) + region1 := CellFromCellID(CellIDFromToken("80c297c574")) + region := region1.CapBound() + cells := []CellID{} // cells := coverer.GetCoveringAsUnion(region) coverer.GetCovering(region, &cells) - for _, cell := range cells { - t.Errorf("cell: %s", cell.ToToken()) - if cell.Level() != 1 { - t.Errorf("Level not as expected %s", cell.ToToken()) - } - if cell.Face() != 4 { - t.Error("Face not as expected") - } + for i, cell := range cells { + t.Errorf("cell %d: %x - %s", i, uint64(cell), cell.ToToken()) } } From 00dcd9bca0e464bee1c09b6af1fa5810296b079a Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Sun, 26 Apr 2015 22:37:24 -0700 Subject: [PATCH 07/23] Bug fix for Metrics --- s2/s2.go | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/s2/s2.go b/s2/s2.go index 9dbd6c79..00a86cac 100644 --- a/s2/s2.go +++ b/s2/s2.go @@ -59,7 +59,10 @@ func (m Metric) GetValue(level int) float64 { * always a valid level. */ func (m Metric) getClosestLevel(value float64) int { - return m.getMinLevel(math.Sqrt2 * value) + if m.dim == 1 { + return m.getMinLevel(math.Sqrt2 * value) + } + return m.getMinLevel(2 * value) } /** From ae63777e9f5e0b145e70dddac54942f2a70b518f Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Wed, 29 Apr 2015 09:32:22 -0700 Subject: [PATCH 08/23] Fixed PriorityQueue in RegionCoverer --- r3/vector.go | 3 +++ s2/cell.go | 4 ++-- s2/cellid.go | 4 ++++ s2/regioncoverer.go | 36 +++++++++++++++++++----------------- s2/regioncoverer_test.go | 18 ++++++++++++++++++ 5 files changed, 46 insertions(+), 19 deletions(-) diff --git a/r3/vector.go b/r3/vector.go index 68aaf4f2..f5b3d81a 100644 --- a/r3/vector.go +++ b/r3/vector.go @@ -59,6 +59,9 @@ func (v Vector) IsUnit() bool { // Abs returns the vector with nonnegative components. func (v Vector) Abs() Vector { return Vector{math.Abs(v.X), math.Abs(v.Y), math.Abs(v.Z)} } +// Neg returns the negated vector +func (v Vector) Neg() Vector { return Vector{-v.X, -v.Y, -v.Z} } + // Add returns the standard vector sum of v and ov. func (v Vector) Add(ov Vector) Vector { return Vector{v.X + ov.X, v.Y + ov.Y, v.Z + ov.Z} } diff --git a/s2/cell.go b/s2/cell.go index 0a218c9d..04afb8b5 100644 --- a/s2/cell.go +++ b/s2/cell.go @@ -98,9 +98,9 @@ func (c Cell) EdgeRaw(k int) Point { case 1: return Point{uNorm(int(c.face), c.uv.X.Hi)} // Right case 2: - return Point{vNorm(int(c.face), c.uv.Y.Hi).Mul(-1.0)} // Top + return Point{vNorm(int(c.face), c.uv.Y.Hi).Neg()} // Top default: - return Point{uNorm(int(c.face), c.uv.X.Lo).Mul(-1.0)} // Left + return Point{uNorm(int(c.face), c.uv.X.Lo).Neg()} // Left } } diff --git a/s2/cellid.go b/s2/cellid.go index 011450bd..138bd378 100644 --- a/s2/cellid.go +++ b/s2/cellid.go @@ -267,6 +267,10 @@ func (ci CellID) String() string { return b.String() } +func (ci CellID) ToString() string { + return fmt.Sprintf("(face=%d, pos=%x, level=%d))", ci.Face(), ci.Pos(), ci.Level()) +} + // Point returns the center of the s2 cell on the sphere as a Point. func (ci CellID) Point() Point { return Point{ci.rawPoint().Normalize()} } diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index 93ce52f3..10925b76 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -48,22 +48,24 @@ type queueEntry struct { index int // The index of the item in the heap. } -func newQueueEntry(id int, candidate *candidate) *queueEntry { - return &queueEntry{priority: id, 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 { - return make([]*queueEntry, 0, space) + pq := PriorityQueue(make([]*queueEntry, 0, space)) + heap.Init(&pq) + return pq } func (pq PriorityQueue) Len() int { return len(pq) } func (pq PriorityQueue) Less(i, j int) bool { // We want Pop to give us the highest, not lowest, priority so we use greater than here. - return pq[i].priority < pq[j].priority + return pq[i].priority > pq[j].priority } func (pq PriorityQueue) Swap(i, j int) { @@ -182,7 +184,7 @@ func (rc *RegionCoverer) GetCoveringInternal(region Region) { rc.getInitialCandidates() for rc.candidateQueue.Len() > 0 && (!rc.interiorCovering || len(rc.result) < rc.maxCells) { - candidate := rc.candidateQueue.Pop().(*queueEntry).candidate + candidate := heap.Pop(&rc.candidateQueue).(*queueEntry).candidate sz := len(rc.result) + candidate.numChildren if !rc.interiorCovering { sz = sz + rc.candidateQueue.Len() @@ -237,32 +239,32 @@ func (rc *RegionCoverer) getInitialCandidates() { } } -func (rc *RegionCoverer) addCandidate(candidate *candidate) { - if candidate == nil { +func (rc *RegionCoverer) addCandidate(candidate_ *candidate) { + if candidate_ == nil { return } - if candidate.isTerminal { - rc.result = append(rc.result, candidate.cell.Id()) + if candidate_.isTerminal { + rc.result = append(rc.result, candidate_.cell.Id()) return } numLevels := rc.levelMod - if int(candidate.cell.Level()) < rc.minLevel { + if int(candidate_.cell.Level()) < rc.minLevel { numLevels = 1 } - numTerminals := rc.expandChildren(candidate, candidate.cell, numLevels) + numTerminals := rc.expandChildren(candidate_, candidate_.cell, numLevels) - if candidate.numChildren == 0 { + if candidate_.numChildren == 0 { // do nothing } else if !rc.interiorCovering && numTerminals == 1<= rc.minLevel { + int(candidate_.cell.Level()) >= rc.minLevel { // Optimization: add the parent cell rather than all of its children. // We can't do this for interior coverings, since the children just // intersect the region, but may not be contained by it - we need to // subdivide them further. - candidate.isTerminal = true - rc.addCandidate(candidate) + candidate_.isTerminal = true + rc.addCandidate(candidate_) } else { // We negate the priority so that smaller absolute priorities are returned // first. The heuristic is designed to refine the largest cells first, @@ -270,8 +272,8 @@ func (rc *RegionCoverer) addCandidate(candidate *candidate) { // at the same level, we prefer the cells with the smallest number of // intersecting children. Finally, we prefer cells that have the smallest // number of children that cannot be refined any further. - priority := -((((int(candidate.cell.Level()) << rc.maxChildrenShift()) + candidate.numChildren) << rc.maxChildrenShift()) + numTerminals) - rc.candidateQueue.Push(newQueueEntry(priority, candidate)) + priority := -((((int(candidate_.cell.Level()) << rc.maxChildrenShift()) + candidate_.numChildren) << rc.maxChildrenShift()) + numTerminals) + heap.Push(&rc.candidateQueue, newQueueEntry(priority, candidate_)) } } diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index 35d7ad39..3841fb1f 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -1,15 +1,33 @@ package s2 import ( + "container/heap" + "fmt" "testing" ) +func TestPriorityQueue(t *testing.T) { + pq := newPriorityQueue(4) + heap.Init(&pq) + for i := -4; i < 0; i++ { + heap.Push(&pq, newQueueEntry(i, nil)) + } + for i := -1; i >= -4; i-- { + entry := heap.Pop(&pq).(*queueEntry) + if entry.priority != i { + t.Errorf("expected %d, got %d", i, entry.priority) + } + } +} + func TestCovering(t *testing.T) { coverer := NewRegionCoverer() coverer.SetMaxCells(8) region1 := CellFromCellID(CellIDFromToken("80c297c574")) region := region1.CapBound() + fmt.Printf("%s\n", region.String()) + cells := []CellID{} // cells := coverer.GetCoveringAsUnion(region) coverer.GetCovering(region, &cells) From 8d24c78e7c32d2573ae514454bc98d170b3e22e5 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Wed, 29 Apr 2015 11:22:04 -0700 Subject: [PATCH 09/23] Added Cell.RectBound support --- r1/interval.go | 11 +++++ s1/interval.go | 18 +++++++ s2/cap.go | 1 - s2/cell.go | 102 ++++++++++++++++++++++++++++++++++++++- s2/regioncoverer_test.go | 5 ++ 5 files changed, 134 insertions(+), 3 deletions(-) diff --git a/r1/interval.go b/r1/interval.go index 41208af4..256effbe 100644 --- a/r1/interval.go +++ b/r1/interval.go @@ -34,6 +34,17 @@ func EmptyInterval() Interval { return Interval{1, 0} } // IntervalFromPoint returns an interval representing a single point. func IntervalFromPoint(p float64) Interval { return Interval{p, p} } +// Convenience method to construct the minimal interval containing the two +// given points. This is equivalent to starting with an empty interval and +// calling AddPoint() twice, but it is more efficient. +func IntervalFromPointPair(p1, p2 float64) Interval { + if p1 <= p2 { + return Interval{p1, p2} + } else { + return Interval{p2, p1} + } +} + // IsEmpty reports whether the interval is empty. func (i Interval) IsEmpty() bool { return i.Lo > i.Hi } diff --git a/s1/interval.go b/s1/interval.go index f6e2b296..f4dcb6cf 100644 --- a/s1/interval.go +++ b/s1/interval.go @@ -47,6 +47,24 @@ func IntervalFromEndpoints(lo, hi float64) Interval { return i } +// Convenience method to construct the minimal interval containing the two +// given points. This is equivalent to starting with an empty interval and +// calling AddPoint() twice, but it is more efficient. +func IntervalFromPointPair(p1, p2 float64) Interval { + // assert (Math.abs(p1) <= S2.M_PI && Math.abs(p2) <= S2.M_PI); + if p1 == -math.Pi { + p1 = math.Pi + } + if p2 == -math.Pi { + p2 = math.Pi + } + if positiveDistance(p1, p2) <= math.Pi { + return Interval{p1, p2} + } else { + return Interval{p2, p1} + } +} + // EmptyInterval returns an empty interval. func EmptyInterval() Interval { return Interval{math.Pi, -math.Pi} } diff --git a/s2/cap.go b/s2/cap.go index 3018f00e..6d4ded4f 100644 --- a/s2/cap.go +++ b/s2/cap.go @@ -401,5 +401,4 @@ func (c Cap) intersects(cell Cell, vertices []Point) bool { } // TODO(roberts): Differences from C++ -// Intersects(S2Cell), Contains(S2Cell), MayIntersect(S2Cell) // Centroid, Union diff --git a/s2/cell.go b/s2/cell.go index 04afb8b5..e35b65ce 100644 --- a/s2/cell.go +++ b/s2/cell.go @@ -17,7 +17,11 @@ limitations under the License. package s2 import ( + "math" + + "github.com/golang/geo/r1" "github.com/golang/geo/r2" + "github.com/golang/geo/s1" ) // Cell is an S2 region object that represents a cell. Unlike CellIDs, @@ -140,11 +144,79 @@ func (c Cell) ContainsPoint(point Point) bool { return u >= c.uv.X.Lo && u <= c.uv.X.Hi && v >= c.uv.Y.Lo && v <= c.uv.Y.Hi } +// We grow the bounds slightly to make sure that the bounding rectangle +// also contains the normalized versions of the vertices. Note that the +// maximum result magnitude is Pi, with a floating-point exponent of 1. +// Therefore adding or subtracting 2**-51 will always change the result. +var MAX_ERROR float64 = 1.0 / (1 << 51) + +// The 4 cells around the equator extend to +/-45 degrees latitude at the +// midpoints of their top and bottom edges. The two cells covering the +// poles extend down to +/-35.26 degrees at their vertices. +// adding kMaxError (as opposed to the C version) because of asin and atan2 +// roundoff errors +var POLE_MIN_LAT float64 = math.Asin(math.Sqrt(1.0/3.0)) - MAX_ERROR // 35.26 degrees + // RectBound returns a bounding latitude-longitude rectangle that contains // the region. The bounds are not guaranteed to be tight. func (c Cell) RectBound() Rect { - // TODO: Implement - return EmptyRect() + if c.level > 0 { + // Except for cells at level 0, the latitude and longitude extremes are + // attained at the vertices. Furthermore, the latitude range is + // determined by one pair of diagonally opposite vertices and the + // longitude range is determined by the other pair. + // + // We first determine which corner (i,j) of the cell has the largest + // absolute latitude. To maximize latitude, we want to find the point in + // the cell that has the largest absolute z-coordinate and the smallest + // absolute x- and y-coordinates. To do this we look at each coordinate + // (u and v), and determine whether we want to minimize or maximize that + // coordinate based on the axis direction and the cell's (u,v) quadrant. + u := c.uv.X.Lo + c.uv.X.Hi + v := c.uv.Y.Lo + c.uv.Y.Hi + var i, j int + if uAxis(int(c.face)).Z == 0 { + if u < 0 { + i = 1 + } + } else { + if u > 0 { + i = 1 + } + } + if vAxis(int(c.face)).Z == 0 { + if v < 0 { + j = 1 + } + } else { + if v > 0 { + j = 1 + } + } + + lat := r1.IntervalFromPointPair(c.latitude(i, j), c.latitude(1-i, 1-j)) + lat = lat.Expanded(MAX_ERROR).Intersection(validRectLatRange) + if lat.Lo == validRectLatRange.Lo || lat.Hi == validRectLatRange.Hi { + return Rect{lat, s1.FullInterval()} + } + lng := s1.IntervalFromPointPair(c.longitude(i, 1-j), c.longitude(1-i, j)) + return Rect{lat, lng.Expanded(MAX_ERROR)} + } + + switch c.face { + case 0: + return Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{-math.Pi / 4, math.Pi / 4}} + case 1: + return Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{math.Pi / 4, 3 * math.Pi / 4}} + case 2: + return Rect{r1.Interval{POLE_MIN_LAT, math.Pi / 2}, s1.Interval{-math.Pi, math.Pi}} + case 3: + return Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{3 * math.Pi / 4, -3 * math.Pi / 4}} + case 4: + return Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{-3 * math.Pi / 4, -math.Pi / 4}} + default: + return Rect{r1.Interval{-math.Pi / 2, -POLE_MIN_LAT}, s1.Interval{-math.Pi, math.Pi}} + } } // ContainsCell reports whether the region completely contains the given region. @@ -159,3 +231,29 @@ func (c Cell) ContainsCell(other Cell) bool { func (c Cell) IntersectsCell(other Cell) bool { return c.Id().Intersects(other.Id()) } + +func (c Cell) latitude(i, j int) float64 { + u := c.uv.X.Lo + if i == 1 { + u = c.uv.X.Hi + } + v := c.uv.Y.Lo + if j == 1 { + v = c.uv.Y.Hi + } + p := Point{faceUVToXYZ(int(c.face), u, v)} + return math.Atan2(p.Z, math.Sqrt(p.X*p.X+p.Y*p.Y)) +} + +func (c Cell) longitude(i, j int) float64 { + u := c.uv.X.Lo + if i == 1 { + u = c.uv.X.Hi + } + v := c.uv.Y.Lo + if j == 1 { + v = c.uv.Y.Hi + } + p := Point{faceUVToXYZ(int(c.face), u, v)} + return math.Atan2(p.Y, p.X) +} diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index 3841fb1f..addec1c2 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -24,7 +24,12 @@ func TestCovering(t *testing.T) { coverer := NewRegionCoverer() coverer.SetMaxCells(8) region1 := CellFromCellID(CellIDFromToken("80c297c574")) + rect1 := region1.RectBound() + fmt.Printf("%s\n", rect1.String()) + region := region1.CapBound() + rect := region.RectBound() + fmt.Printf("%s\n", rect.String()) fmt.Printf("%s\n", region.String()) From d80b92c8eca928e979ebec95c60d608b7e0d4497 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Wed, 29 Apr 2015 18:03:08 -0700 Subject: [PATCH 10/23] Added more stuff --- s2/cell.go | 4 +- s2/edgeutil.go | 146 +++++++++++++++++++++++++++++++++++++++ s2/latlng.go | 25 ++----- s2/point.go | 8 +++ s2/polyline.go | 128 ++++++++++++++++++++++++++++++++++ s2/rect.go | 117 ++++++++++++++++++++++++++++++- s2/region.go | 6 +- s2/regioncoverer_test.go | 27 ++++++++ 8 files changed, 437 insertions(+), 24 deletions(-) create mode 100644 s2/edgeutil.go create mode 100644 s2/polyline.go diff --git a/s2/cell.go b/s2/cell.go index e35b65ce..3bb98067 100644 --- a/s2/cell.go +++ b/s2/cell.go @@ -242,7 +242,7 @@ func (c Cell) latitude(i, j int) float64 { v = c.uv.Y.Hi } p := Point{faceUVToXYZ(int(c.face), u, v)} - return math.Atan2(p.Z, math.Sqrt(p.X*p.X+p.Y*p.Y)) + return latitude(p).Radians() } func (c Cell) longitude(i, j int) float64 { @@ -255,5 +255,5 @@ func (c Cell) longitude(i, j int) float64 { v = c.uv.Y.Hi } p := Point{faceUVToXYZ(int(c.face), u, v)} - return math.Atan2(p.Y, p.X) + return longitude(p).Radians() } diff --git a/s2/edgeutil.go b/s2/edgeutil.go new file mode 100644 index 00000000..0532351f --- /dev/null +++ b/s2/edgeutil.go @@ -0,0 +1,146 @@ +package s2 + +import ( + "fmt" + "math" + + "github.com/golang/geo/r1" +) + +type EdgeCrosser struct { + // The fields below are all constant. + a Point + b 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, + } + ec.RestartAt(c) + return ec +} + +func (ec EdgeCrosser) RestartAt(c Point) { + ec.c = c + ec.acb = -int(RobustSign(ec.a, ec.b, c)) +} + +/** + * 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(RobustSign(ec.a, ec.b, d)) + 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 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. + cbd := -int(RobustSign(ec.c, d, ec.b)) + if cbd != ec.acb { + return -1 + } + + dac := int(RobustSign(ec.c, d, ec.a)) + 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() { + fmt.Printf("AddPoint %s\n", rb.bound.String()) + rb.bound = rb.bound.AddPoint(bLatLng) + fmt.Printf("AddPoint %s\n", rb.bound.String()) + } 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(RectFromLatLngPair(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(PointFromCoords(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 } diff --git a/s2/latlng.go b/s2/latlng.go index 6111b9b4..8dea3867 100644 --- a/s2/latlng.go +++ b/s2/latlng.go @@ -33,6 +33,11 @@ func LatLngFromDegrees(lat, lng float64) LatLng { return LatLng{s1.Angle(lat) * s1.Degree, s1.Angle(lng) * s1.Degree} } +// LatLngFromPoint returns an LatLng for a given Point. +func LatLngFromPoint(p Point) LatLng { + return LatLng{latitude(p), longitude(p)} +} + // IsValid returns true iff the LatLng is normalized, with Lat ∈ [-π/2,π/2] and Lng ∈ [-π,π]. func (ll LatLng) IsValid() bool { return math.Abs(ll.Lat.Radians()) <= math.Pi/2 && math.Abs(ll.Lng.Radians()) <= math.Pi @@ -51,28 +56,12 @@ func (ll LatLng) Distance(ll2 LatLng) s1.Angle { return s1.Angle(2*math.Atan2(math.Sqrt(x), math.Sqrt(math.Max(0, 1-x)))) * s1.Radian } -// 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)) -} - -// LatLngFromPoint returns an LatLng for a given Point. -func LatLngFromPoint(p Point) LatLng { - return LatLng{latitude(p), longitude(p)} + return s1.Angle(math.Atan2(p.Y, p.X)) } // BUG(dsymonds): The major differences from the C++ version are: diff --git a/s2/point.go b/s2/point.go index 9f90db92..a06331d0 100644 --- a/s2/point.go +++ b/s2/point.go @@ -83,6 +83,14 @@ func PointFromCoords(x, y, z float64) Point { return Point{r3.Vector{x, y, z}.Normalize()} } +// 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). diff --git a/s2/polyline.go b/s2/polyline.go new file mode 100644 index 00000000..b686fcc7 --- /dev/null +++ b/s2/polyline.go @@ -0,0 +1,128 @@ +/* +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].ApproxEqual(p.Vertices[i]) || p.Vertices[i-1].ApproxEqual(Point{p.Vertices[i].Neg()}) { + 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()-1; i++ { + fmt.Printf("polyline Rectbound %d\n", 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/rect.go b/s2/rect.go index 56af72c3..541a0b8b 100644 --- a/s2/rect.go +++ b/s2/rect.go @@ -49,6 +49,21 @@ func RectFromLatLng(p LatLng) Rect { } } +/** + * 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 RectFromLatLngPair(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 @@ -58,7 +73,7 @@ func RectFromLatLng(p LatLng) Rect { // Examples of clamping (in degrees): // center=(80,170), size=(40,60) -> lat=[60,90], lng=[140,-160] // center=(10,40), size=(210,400) -> lat=[-90,90], lng=[-180,180] -// center=(-90,180), size=(20,50) -> lat=[-90,-80], lng=[155,-155] +// center=(-90,180), size=(20,50) -> lat=[-90,-80], lng=[1 func RectFromCenterSize(center, size LatLng) Rect { half := LatLng{size.Lat / 2, size.Lng / 2} return RectFromLatLng(center).expanded(half) @@ -92,6 +107,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 +151,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 +193,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(PointFromCoords(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..431bed69 100644 --- a/s2/region.go +++ b/s2/region.go @@ -40,10 +40,10 @@ 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{} ) -*/ diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index addec1c2..ac28c47e 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -40,3 +40,30 @@ func TestCovering(t *testing.T) { t.Errorf("cell %d: %x - %s", i, uint64(cell), cell.ToToken()) } } + +func TestCoveringPolyline(t *testing.T) { + coverer := NewRegionCoverer() + coverer.SetMaxCells(8) + + points := []Point{ + PointFromLatLng(LatLngFromDegrees(34.0909533022671600, -118.3914214745164100)), + PointFromLatLng(LatLngFromDegrees(34.0906409358360560, -118.3911871165037200)), + } + + for _, point := range points { + fmt.Printf("point: %v\n", point.String()) + } + + // region := RectFromLatLng(LatLngFromDegrees(34.0909533022671600, -118.3914214745164100)) + // region = region.AddPoint(LatLngFromDegrees(34.0906409358360560, -118.3911871165037200)) + + polyline := PolylineFromPoints(points) + region := polyline.RectBound() + fmt.Printf("%s\n", region.String()) + + cells := []CellID{} + coverer.GetCovering(region, &cells) + for i, cell := range cells { + t.Errorf("cell %d: %x - %s", i, uint64(cell), cell.ToToken()) + } +} From 7bf21c9363355050d9a32bdd2197bc1720bca365 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Wed, 29 Apr 2015 22:24:56 -0700 Subject: [PATCH 11/23] Fixed bugs in edgeutil --- s2/edgeutil.go | 11 ++++------- s2/polyline.go | 3 +-- s2/regioncoverer_test.go | 6 +++--- 3 files changed, 8 insertions(+), 12 deletions(-) diff --git a/s2/edgeutil.go b/s2/edgeutil.go index 0532351f..79a10179 100644 --- a/s2/edgeutil.go +++ b/s2/edgeutil.go @@ -1,7 +1,6 @@ package s2 import ( - "fmt" "math" "github.com/golang/geo/r1" @@ -29,7 +28,7 @@ func NewEdgeCrosser(a, b, c Point) EdgeCrosser { return ec } -func (ec EdgeCrosser) RestartAt(c Point) { +func (ec *EdgeCrosser) RestartAt(c Point) { ec.c = c ec.acb = -int(RobustSign(ec.a, ec.b, c)) } @@ -42,7 +41,7 @@ func (ec EdgeCrosser) RestartAt(c Point) { * degenerate. As a side effect, it saves vertex D to be used as the next * vertex C. */ -func (ec EdgeCrosser) RobustCrossing(d Point) int { +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 @@ -76,7 +75,7 @@ func (ec EdgeCrosser) RobustCrossing(d Point) int { /** * This function handles the "slow path" of robustCrossing(). */ -func (ec EdgeCrosser) robustCrossingInternal(d Point) int { +func (ec *EdgeCrosser) robustCrossingInternal(d Point) int { // ACB and BDA have the appropriate orientations, so now we check the // triangles CBD and DAC. cbd := -int(RobustSign(ec.c, d, ec.b)) @@ -102,13 +101,11 @@ func NewRectBounder() RectBounder { return RectBounder{bound: EmptyRect()} } -func (rb RectBounder) AddPoint(b Point) { +func (rb *RectBounder) AddPoint(b Point) { bLatLng := LatLngFromPoint(b) if rb.bound.IsEmpty() { - fmt.Printf("AddPoint %s\n", rb.bound.String()) rb.bound = rb.bound.AddPoint(bLatLng) - fmt.Printf("AddPoint %s\n", rb.bound.String()) } 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. diff --git a/s2/polyline.go b/s2/polyline.go index b686fcc7..3b9f262b 100644 --- a/s2/polyline.go +++ b/s2/polyline.go @@ -81,8 +81,7 @@ func (p Polyline) CapBound() Cap { // the region. The bounds are not guaranteed to be tight. func (p Polyline) RectBound() Rect { rb := NewRectBounder() - for i := 0; i < p.NumVertices()-1; i++ { - fmt.Printf("polyline Rectbound %d\n", i) + for i := 0; i < p.NumVertices(); i++ { rb.AddPoint(p.Vertex(i)) } return rb.GetBound() diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index ac28c47e..d6f02338 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -58,11 +58,11 @@ func TestCoveringPolyline(t *testing.T) { // region = region.AddPoint(LatLngFromDegrees(34.0906409358360560, -118.3911871165037200)) polyline := PolylineFromPoints(points) - region := polyline.RectBound() - fmt.Printf("%s\n", region.String()) + // region := polyline.RectBound() + // fmt.Printf("%s\n", region.String()) cells := []CellID{} - coverer.GetCovering(region, &cells) + coverer.GetCovering(polyline, &cells) for i, cell := range cells { t.Errorf("cell %d: %x - %s", i, uint64(cell), cell.ToToken()) } From 9598df8b1feea38ccb37b4eb69274c9171fe2063 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Thu, 30 Apr 2015 19:55:52 -0700 Subject: [PATCH 12/23] Added S2Loop and some more functionality and fixes, still work in progress --- s2/edge.go | 25 +++ s2/edgeutil.go | 220 +++++++++++++++++++- s2/loop.go | 429 +++++++++++++++++++++++++++++++++++++++ s2/point.go | 40 ++++ s2/region.go | 1 + s2/regioncoverer.go | 2 +- s2/regioncoverer_test.go | 2 +- 7 files changed, 712 insertions(+), 7 deletions(-) create mode 100644 s2/edge.go create mode 100644 s2/loop.go diff --git a/s2/edge.go b/s2/edge.go new file mode 100644 index 00000000..63938e76 --- /dev/null +++ b/s2/edge.go @@ -0,0 +1,25 @@ +package s2 + +import ( + "fmt" +) + +type Edge struct { + start Point + end Point +} + +func NewEdgeFromStartEnd(start, end Point) Edge { + return Edge{start, end} +} + +func (e Edge) Start() Point { return e.start } +func (e Edge) End() Point { return e.end } + +func (e Edge) String() string { + return fmt.Sprintf("Edge: (%s -> %s)\n or [%s -> %s]", e.start.DegreesString(), e.end.DegreesString(), e.start.String(), e.end.String()) +} + +func (e Edge) Equals(other Edge) bool { + return e.start.ApproxEqual(other.start) && e.end.ApproxEqual(other.end) +} diff --git a/s2/edgeutil.go b/s2/edgeutil.go index 79a10179..92a22b76 100644 --- a/s2/edgeutil.go +++ b/s2/edgeutil.go @@ -19,8 +19,8 @@ type EdgeCrosser struct { acb int } -func NewEdgeCrosser(a, b, c Point) EdgeCrosser { - ec := EdgeCrosser{ +func NewEdgeCrosser(a, b, c Point) *EdgeCrosser { + ec := &EdgeCrosser{ a: a, b: b, } @@ -72,6 +72,27 @@ func (ec *EdgeCrosser) RobustCrossing(d Point) int { 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 := Point{ec.c.Vector} + + 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(). */ @@ -97,8 +118,8 @@ type RectBounder struct { bound Rect } -func NewRectBounder() RectBounder { - return RectBounder{bound: EmptyRect()} +func NewRectBounder() *RectBounder { + return &RectBounder{bound: EmptyRect()} } func (rb *RectBounder) AddPoint(b Point) { @@ -140,4 +161,193 @@ func (rb *RectBounder) AddPoint(b Point) { rb.aLatLng = bLatLng } -func (rb RectBounder) GetBound() Rect { return rb.bound } +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. +} diff --git a/s2/loop.go b/s2/loop.go new file mode 100644 index 00000000..ddcdd16b --- /dev/null +++ b/s2/loop.go @@ -0,0 +1,429 @@ +package s2 + +import ( + "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, 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) 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(PointFromCoords(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(PointFromCoords(0, 0, -1)) { + b = Rect{r1.IntervalFromPointPair(-math.Pi/2, b.Lat.Hi), b.Lng} + } + l.bound = b +} + +/** + * 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 +} + +// CapBound returns a bounding spherical cap. This is not guaranteed to be exact. +func (l *Loop) CapBound() Cap { + // TODO: Implement + return EmptyCap() +} + +// 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 := LoopFromCell(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 LoopFromCell(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 := PointFromCoords(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 { + panic("TODO") + // DataEdgeIterator it = getEdgeIterator(numVertices) + // int previousIndex = -2 + // for (it.getCandidates(origin, p); it.hasNext(); it.next()) { + // int ai = it.index() + // if (previousIndex != ai - 1) { + // crosser.restartAt(vertices[ai]) + // } + // previousIndex = ai + // inside ^= crosser.EdgeOrVertexCrossing(vertex(ai + 1)) + // } + } + + return inside +} + +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 +} + +/** + * 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 { + // DataEdgeIterator it = getEdgeIterator(b.numVertices); + result := 1 + // since 'this' usually has many more vertices than 'b', use the index on + // 'this' and loop over 'b' + // for (int j = 0; j < b.numVertices(); ++j) { + // S2EdgeUtil.EdgeCrosser crosser = + // new S2EdgeUtil.EdgeCrosser(b.vertex(j), b.vertex(j + 1), vertex(0)); + // int previousIndex = -2; + // for (it.getCandidates(b.vertex(j), b.vertex(j + 1)); it.hasNext(); it.next()) { + // int i = it.index(); + // if (previousIndex != i - 1) { + // crosser.restartAt(vertex(i)); + // } + // previousIndex = i; + // int crossing = crosser.robustCrossing(vertex(i + 1)); + // if (crossing < 0) { + // continue; + // } + // if (crossing > 0) { + // return -1; // There is a proper edge crossing. + // } + // if (vertex(i + 1).equals(b.vertex(j + 1))) { + // result = Math.min(result, relation.test( + // vertex(i), vertex(i + 1), vertex(i + 2), b.vertex(j), b.vertex(j + 2))); + // if (result < 0) { + // return result; + // } + // } + // } + // } + return result +} diff --git a/s2/point.go b/s2/point.go index a06331d0..d44f362e 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" @@ -278,6 +279,45 @@ func (p Point) ApproxEqual(other Point) bool { return p.Vector.Angle(other.Vector) <= epsilon } +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 vb.X < p.X { + return false + } + if p.Y < vb.Y { + return true + } + if vb.Y < p.Y { + return false + } + if p.Z < vb.Z { + return true + } + return false +} + +func (p Point) CompareTo(other Point) int { + if p.LessThan(other) { + return -1 + } else { + if p.Equals(other) { + return 0 + } + return 1 + } +} + +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 // given points. // diff --git a/s2/region.go b/s2/region.go index 431bed69..f0763d4e 100644 --- a/s2/region.go +++ b/s2/region.go @@ -46,4 +46,5 @@ var ( // _ Region = CellUnion(nil) _ Region = Rect{} _ Region = Polyline{} + _ Region = &Loop{} ) diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index 10925b76..2b24907d 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -65,7 +65,7 @@ func (pq PriorityQueue) Len() int { return len(pq) } func (pq PriorityQueue) Less(i, j int) bool { // We want Pop to give us the highest, not lowest, priority so we use greater than here. - return pq[i].priority > pq[j].priority + return pq[i].priority >= pq[j].priority } func (pq PriorityQueue) Swap(i, j int) { diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index d6f02338..437afcfa 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -43,7 +43,7 @@ func TestCovering(t *testing.T) { func TestCoveringPolyline(t *testing.T) { coverer := NewRegionCoverer() - coverer.SetMaxCells(8) + coverer.SetMaxCells(4) points := []Point{ PointFromLatLng(LatLngFromDegrees(34.0909533022671600, -118.3914214745164100)), From 0397a454db1d5de98f449997d56339d56aae3372 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Sat, 2 May 2015 00:11:33 -0700 Subject: [PATCH 13/23] Added edgeindex --- r3/vector.go | 5 +- s2/cap.go | 13 +- s2/cellid.go | 5 + s2/edge.go | 2 +- s2/edgeindex.go | 590 +++++++++++++++++++++++++++++++++++++++++++++++ s2/edgeutil.go | 61 +++++ s2/latlng.go | 4 + s2/loop.go | 540 +++++++++++++++++++++++++++++++++---------- s2/point.go | 3 +- s2/point_test.go | 6 +- s2/polyline.go | 2 +- s2/s2.go | 45 ++++ 12 files changed, 1142 insertions(+), 134 deletions(-) create mode 100644 s2/edgeindex.go diff --git a/r3/vector.go b/r3/vector.go index f5b3d81a..515f8e65 100644 --- a/r3/vector.go +++ b/r3/vector.go @@ -69,7 +69,10 @@ func (v Vector) Add(ov Vector) Vector { return Vector{v.X + ov.X, v.Y + ov.Y, v. func (v Vector) Sub(ov Vector) Vector { return Vector{v.X - ov.X, v.Y - ov.Y, v.Z - ov.Z} } // Mul returns the standard scalar product of v and m. -func (v Vector) Mul(m float64) Vector { return Vector{m * v.X, m * v.Y, m * v.Z} } +func (v Vector) Mul(m float64) Vector { return Vector{v.X * m, v.Y * m, v.Z * m} } + +// Mul returns the standard scalar product of v and m. +func (v Vector) Div(m float64) Vector { return Vector{v.X / m, v.Y / m, v.Z / m} } // Dot returns the standard dot product of v and ov. func (v Vector) Dot(ov Vector) float64 { return v.X*ov.X + v.Y*ov.Y + v.Z*ov.Z } diff --git a/s2/cap.go b/s2/cap.go index 6d4ded4f..019162b1 100644 --- a/s2/cap.go +++ b/s2/cap.go @@ -249,13 +249,12 @@ func (c Cap) RectBound() Rect { // ApproxEqual reports if this caps' center and height are within // a reasonable epsilon from the other cap. func (c Cap) ApproxEqual(other Cap) bool { - const epsilon = 1e-14 - return c.center.ApproxEqual(other.center) && - math.Abs(c.height-other.height) <= epsilon || - c.IsEmpty() && other.height <= epsilon || - other.IsEmpty() && c.height <= epsilon || - c.IsFull() && other.height >= 2-epsilon || - other.IsFull() && c.height >= 2-epsilon + return c.center.ApproxEquals(other.center, EPSILON) && + math.Abs(c.height-other.height) <= EPSILON || + c.IsEmpty() && other.height <= EPSILON || + other.IsEmpty() && c.height <= EPSILON || + c.IsFull() && other.height >= 2-EPSILON || + other.IsFull() && c.height >= 2-EPSILON } // AddPoint increases the cap if necessary to include the given point. If this cap is empty, diff --git a/s2/cellid.go b/s2/cellid.go index 138bd378..d8aad72f 100644 --- a/s2/cellid.go +++ b/s2/cellid.go @@ -47,6 +47,11 @@ const ( maxSize = 1 << maxLevel ) +func CellIDNone() CellID { return CellID(0) } +func CellIDSentinel() CellID { return CellID(^uint64(0)) } +func CellIDBegin(level int) CellID { return CellIDFromFacePosLevel(0, 0, 0).ChildBeginAtLevel(level) } +func CellIDEnd(level int) CellID { return CellIDFromFacePosLevel(5, 0, 0).ChildBeginAtLevel(level) } + // CellIDFromFacePosLevel returns a cell given its face in the range // [0,5], the 61-bit Hilbert curve position pos within that face, and // the level in the range [0,maxLevel]. The position in the cell ID diff --git a/s2/edge.go b/s2/edge.go index 63938e76..e2fb788c 100644 --- a/s2/edge.go +++ b/s2/edge.go @@ -21,5 +21,5 @@ func (e Edge) String() string { } func (e Edge) Equals(other Edge) bool { - return e.start.ApproxEqual(other.start) && e.end.ApproxEqual(other.end) + return e.start.ApproxEquals(other.start, EPSILON) && e.end.ApproxEquals(other.end, EPSILON) } diff --git a/s2/edgeindex.go b/s2/edgeindex.go new file mode 100644 index 00000000..e07ed190 --- /dev/null +++ b/s2/edgeindex.go @@ -0,0 +1,590 @@ +package s2 + +import ( + "math" + "sort" +) + +var THICKENING float64 = 0.01 +var MAX_DET_ERROR float64 = 1e-14 + +type EdgeIndex struct { + cells []uint64 + edges []int + minimumS2LevelUsed int + indexComputed bool + queryCount int + + getNumEdges func() int + edgeFrom func(int) Point + edgeTo func(int) Point +} + +func NewEdgeIndex(getNumEdges func() int, edgeFrom func(int) Point, edgeTo func(int) Point) *EdgeIndex { + return &EdgeIndex{ + getNumEdges: getNumEdges, + edgeFrom: edgeFrom, + edgeTo: edgeTo, + } +} + +/** + * Empties the index in case it already contained something. + */ +func (e *EdgeIndex) Reset() { + e.minimumS2LevelUsed = MAX_LEVEL + e.indexComputed = false + e.queryCount = 0 + e.cells = nil + e.edges = nil +} + +/** + * Compares [cell1, edge1] to [cell2, edge2], by cell first and edge second. + * + * @return -1 if [cell1, edge1] is less than [cell2, edge2], 1 if [cell1, + * edge1] is greater than [cell2, edge2], 0 otherwise. + */ +func compare(cell1 uint64, edge1 int, cell2 uint64, edge2 int) int { + if cell1 < cell2 { + return -1 + } else if cell1 > cell2 { + return 1 + } else if edge1 < edge2 { + return -1 + } else if edge1 > edge2 { + return 1 + } else { + return 0 + } +} + +/** Computes the index (if it has not been previously done). */ +func (e *EdgeIndex) ComputeIndex() { + if e.indexComputed { + return + } + cellList := []uint64{} + edgeList := []int{} + for i := 0; i < e.getNumEdges(); i++ { + from := e.edgeFrom(i) + to := e.edgeTo(i) + cover := []CellID{} + level := e.getCovering(from, to, true, &cover) + e.minimumS2LevelUsed = min(e.minimumS2LevelUsed, level) + for _, cellId := range cover { + cellList = append(cellList, uint64(cellId)) + edgeList = append(edgeList, i) + } + } + e.cells = make([]uint64, len(cellList)) + e.edges = make([]int, len(edgeList)) + for i := 0; i < len(e.cells); i++ { + e.cells[i] = cellList[i] + e.edges[i] = edgeList[i] + } + e.sortIndex() + e.indexComputed = true +} + +type ByCellThenEdge struct { + e *EdgeIndex + a []int +} + +func (s ByCellThenEdge) Len() int { return len(s.a) } +func (s ByCellThenEdge) Swap(i, j int) { s.a[i], s.a[j] = s.a[j], s.a[i] } +func (s ByCellThenEdge) Less(i, j int) bool { + return compare(s.e.cells[i], s.e.edges[i], s.e.cells[j], s.e.edges[j]) < 0 +} + +/** Sorts the parallel 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, -(int(^uint(0)>>1)-1)), + -1 - e.binarySearch(cell2, int(^uint(0)>>1)), + } +} + +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 index 92a22b76..1e1d6eff 100644 --- a/s2/edgeutil.go +++ b/s2/edgeutil.go @@ -4,6 +4,7 @@ import ( "math" "github.com/golang/geo/r1" + "github.com/golang/geo/s1" ) type EdgeCrosser struct { @@ -351,3 +352,63 @@ func (w WedgeContainsOrCrosses) Test(a0, ab1, a2, b0, b2 Point) int { } 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 8dea3867..8ccb636c 100644 --- a/s2/latlng.go +++ b/s2/latlng.go @@ -45,6 +45,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. diff --git a/s2/loop.go b/s2/loop.go index ddcdd16b..4d65fd4a 100644 --- a/s2/loop.go +++ b/s2/loop.go @@ -1,6 +1,7 @@ package s2 import ( + "fmt" "math" "github.com/golang/geo/r1" @@ -33,7 +34,7 @@ type Loop struct { // 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 + index *EdgeIndex // Maps each S2Point to its order in the loop, from 1 to numVertices. vertexToIndex map[Point]int @@ -61,7 +62,11 @@ func LoopFromPoints(points []Point) *Loop { return l } -func LoopFromCell(cell Cell, bound Rect) *Loop { +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), @@ -120,59 +125,25 @@ func (l *Loop) Vertex(i int) Point { return l.vertices[i] } -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)") +func (l *Loop) CompareTo(other *Loop) int { + if l.NumVertices() != other.NumVertices() { + return l.NumVertices() - other.NumVertices() } - - // 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(PointFromCoords(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(PointFromCoords(0, 0, -1)) { - b = Rect{r1.IntervalFromPointPair(-math.Pi/2, b.Lat.Hi), b.Lng} + // 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++ } - l.bound = b + return 0 } /** @@ -189,86 +160,54 @@ func (l *Loop) initFirstLogicalVertex() { l.firstLogicalVertex = first } -// CapBound returns a bounding spherical cap. This is not guaranteed to be exact. -func (l *Loop) CapBound() Cap { - // TODO: Implement - return EmptyCap() -} - -// 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. + * Return true if the loop area is at most 2*Pi. */ -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 := LoopFromCell(cell, cellBound) - return l.ContainsLoop(cellLoop) +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 } /** - * 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. + * Invert the loop if necessary so that the area enclosed by the loop is at + * most 2*Pi. */ -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 +func (l *Loop) Normalize() { + if !l.IsNormalized() { + l.Invert() } - return LoopFromCell(cell, cellBound).IntersectsLoop(l) } /** - * The point 'p' does not need to be normalized. + * Reverse the order of the loop vertices, effectively complementing the + * region represented by the loop. */ -func (l *Loop) ContainsPoint(p Point) bool { - if !l.bound.ContainsLatLng(LatLngFromPoint(p)) { - return false +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 } - - inside := l.originInside - origin := PointFromCoords(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)) - } + 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 { - panic("TODO") - // DataEdgeIterator it = getEdgeIterator(numVertices) - // int previousIndex = -2 - // for (it.getCandidates(origin, p); it.hasNext(); it.next()) { - // int ai = it.index() - // if (previousIndex != ai - 1) { - // crosser.restartAt(vertices[ai]) - // } - // previousIndex = ai - // inside ^= crosser.EdgeOrVertexCrossing(vertex(ai + 1)) - // } + l.initBound() } + l.initFirstLogicalVertex() +} - return inside +/** + * 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 1 //l.GetAreaCentroid(false).getArea() } func (l *Loop) ContainsLoop(b *Loop) bool { @@ -346,7 +285,7 @@ func (l *Loop) IntersectsLoop(b *Loop) bool { // Now check whether there are any edge crossings, and also check the loop // relationship at any shared vertices. - if (l.checkEdgeCrossings(b, WedgeIntersects{}) < 0) { + if l.checkEdgeCrossings(b, WedgeIntersects{}) < 0 { return true } @@ -365,6 +304,369 @@ func (l *Loop) IntersectsLoop(b *Loop) bool { 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 := PointFromCoords(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 { + panic("TODO") + // DataEdgeIterator it = getEdgeIterator(numVertices) + // int previousIndex = -2 + // for (it.getCandidates(origin, p); it.hasNext(); it.next()) { + // int ai = it.index() + // if (previousIndex != ai - 1) { + // crosser.restartAt(vertices[ai]) + // } + // previousIndex = ai + // inside ^= crosser.EdgeOrVertexCrossing(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", 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", 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", a1, b1) + fmt.Printf("Edge locations in degrees: %s-%s and %s-%s", + 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(PointFromCoords(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(PointFromCoords(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. diff --git a/s2/point.go b/s2/point.go index d44f362e..d573eb73 100644 --- a/s2/point.go +++ b/s2/point.go @@ -274,8 +274,7 @@ func (p Point) Distance(b Point) s1.Angle { } // ApproxEqual reports if the two points are similar enough to be equal. -func (p Point) ApproxEqual(other Point) bool { - const epsilon = 1e-14 +func (p Point) ApproxEquals(other Point, maxError float64) bool { return p.Vector.Angle(other.Vector) <= epsilon } 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 index 3b9f262b..b9f4e06d 100644 --- a/s2/polyline.go +++ b/s2/polyline.go @@ -51,7 +51,7 @@ func (p Polyline) IsValid() bool { } // Adjacent vertices must not be identical or antipodal. for i := 0; i < n; i++ { - if p.Vertices[i-1].ApproxEqual(p.Vertices[i]) || p.Vertices[i-1].ApproxEqual(Point{p.Vertices[i].Neg()}) { + 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 } diff --git a/s2/s2.go b/s2/s2.go index 00a86cac..494b1ac2 100644 --- a/s2/s2.go +++ b/s2/s2.go @@ -4,6 +4,8 @@ import ( "math" ) +const EPSILON float64 = 1e-14 + // Number of bits in the mantissa of a double. const EXPONENT_SHIFT uint = 52 @@ -32,6 +34,49 @@ func exp(v float64) int { return (int)((EXPONENT_MASK&bits)>>EXPONENT_SHIFT) - 1022 } +/** + * Return the angle at the vertex B in the triangle ABC. The return value is + * always in the range [0, Pi]. The points do not need to be normalized. + * Ensures that Angle(a,b,c) == Angle(c,b,a) for all a,b,c. + * + * The angle is undefined if A or C is diametrically opposite from B, and + * becomes numerically unstable as the length of edge AB or BC approaches 180 + * degrees. + */ +func angle(a, b, c Point) float64 { + return a.Cross(b.Vector).Angle(c.Cross(b.Vector)).Radians() +} + +func approxEqualsNumber(a, b, maxError float64) bool { + return math.Abs(a-b) <= maxError +} + +/** + * Return true if the points A, B, C are strictly counterclockwise. Return + * false if the points are clockwise or colinear (i.e. if they are all + * contained on some great circle). + * + * Due to numerical errors, situations may arise that are mathematically + * impossible, e.g. ABC may be considered strictly CCW while BCA is not. + * However, the implementation guarantees the following: + * + * If SimpleCCW(a,b,c), then !SimpleCCW(c,b,a) for all a,b,c. + * + * In other words, ABC and CBA are guaranteed not to be both CCW + */ +func simpleCCW(a, b, c Point) bool { + // We compute the signed volume of the parallelepiped ABC. The usual + // formula for this is (AxB).C, but we compute it here using (CxA).B + // in order to ensure that ABC and CBA are not both CCW. This follows + // from the following identities (which are true numerically, not just + // mathematically): + // + // (1) x.CrossProd(y) == -(y.CrossProd(x)) + // (2) (-x).DotProd(y) == -(x.DotProd(y)) + + return c.Cross(a.Vector).Dot(b.Vector) > 0 +} + // Defines an area or a length cell metric. type Metric struct { deriv float64 From dcaf8e91c992c172aa033c9a3c766371c14eab2a Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Mon, 4 May 2015 13:02:59 -0700 Subject: [PATCH 14/23] Loop supports Region interface --- r2/vector.go | 104 +++++++++++ r3/vector.go | 31 ++++ s2/areacentroid.go | 33 ++++ s2/loop.go | 171 +++++++++++++----- s2/point.go | 27 --- s2/s2.go | 424 +++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 724 insertions(+), 66 deletions(-) create mode 100644 r2/vector.go create mode 100644 s2/areacentroid.go diff --git a/r2/vector.go b/r2/vector.go new file mode 100644 index 00000000..ab4a87be --- /dev/null +++ b/r2/vector.go @@ -0,0 +1,104 @@ +/* + * Copyright 2005 Google Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package r2 + +import ( + "fmt" + "math" +) + +/** + * r2.Vector represents a vector in the two-dimensional space. It defines the + * basic geometrical operations for 2D vectors, e.g. cross product, addition, + * norm, comparison etc. + * + */ +type Vector struct { + X, Y float64 +} + +func (v Vector) String() string { return fmt.Sprintf("(%v, %v)", v.X, v.Y) } + +// Norm returns the vector's norm. +func (v Vector) Norm() float64 { return math.Sqrt(v.Dot(v)) } + +// Norm2 returns the square of the norm. +func (v Vector) Norm2() float64 { return v.Dot(v) } + +// Normalize returns a unit vector in the same direction as v. +func (v Vector) Normalize() Vector { + if v == (Vector{0, 0}) { + return v + } + return v.Mul(1 / v.Norm()) +} + +// Abs returns the vector with nonnegative components. +func (v Vector) Abs() Vector { return Vector{math.Abs(v.X), math.Abs(v.Y)} } + +// Neg returns the negated vector +func (v Vector) Neg() Vector { return Vector{-v.X, -v.Y} } + +// Add returns the standard vector sum of v and ov. +func (v Vector) Add(ov Vector) Vector { return Vector{v.X + ov.X, v.Y + ov.Y} } + +// Sub returns the standard vector difference of v and ov. +func (v Vector) Sub(ov Vector) Vector { return Vector{v.X - ov.X, v.Y - ov.Y} } + +// Mul returns the standard scalar product of v and m. +func (v Vector) Mul(m float64) Vector { return Vector{v.X * m, v.Y * m} } + +// Mul returns the standard scalar product of v and m. +func (v Vector) Div(m float64) Vector { return Vector{v.X / m, v.Y / m} } + +// Dot returns the standard dot product of v and ov. +func (v Vector) Dot(ov Vector) float64 { return v.X*ov.X + v.Y*ov.Y } + +// Cross returns the standard cross product of v and ov. +func (v Vector) Cross(ov Vector) float64 { + return v.X*ov.Y - v.Y*ov.X +} + +func (v Vector) Equals(other Vector) bool { + return v.X == other.X && v.Y == other.Y +} + +func (v Vector) LessThan(vb Vector) bool { + if v.X < vb.X { + return true + } + if vb.X < v.X { + return false + } + if v.Y < vb.Y { + return true + } + if vb.Y < v.Y { + return false + } + return false +} + +func (v Vector) CompareTo(other Vector) int { + if v.LessThan(other) { + return -1 + } else { + if v.Equals(other) { + return 0 + } + return 1 + } +} diff --git a/r3/vector.go b/r3/vector.go index 515f8e65..1b2e87ac 100644 --- a/r3/vector.go +++ b/r3/vector.go @@ -28,6 +28,37 @@ type Vector struct { X, Y, Z float64 } +func (v Vector) GetAxis(axis int) float64 { + switch axis { + case 0: + return v.X + case 1: + return v.Y + case 2: + return v.Z + default: + return v.GetAxis((axis + 3) % 3) + } +} + +/** Return the index of the largest component fabs */ +func (v Vector) LargestAbsComponent() int { + temp := v.Abs() + if temp.X > temp.Y { + if temp.X > temp.Z { + return 0 + } else { + return 2 + } + } else { + if temp.Y > temp.Z { + return 1 + } else { + return 2 + } + } +} + // ApproxEqual reports whether v and ov are equal within a small epsilon. func (v Vector) ApproxEqual(ov Vector) bool { const epsilon = 1e-14 diff --git a/s2/areacentroid.go b/s2/areacentroid.go new file mode 100644 index 00000000..094d0bf9 --- /dev/null +++ b/s2/areacentroid.go @@ -0,0 +1,33 @@ +package s2 + +/** + * The area of an interior, i.e. the region on the left side of an odd + * number of loops and optionally a centroid. + * The area is between 0 and 4*Pi. If it has a centroid, it is + * the true centroid of the interiord multiplied by the area of the shape. + * Note that the centroid may not be contained by the shape. + * + * @author dbentley@google.com (Daniel Bentley) + */ +type AreaCentroid struct { + area float64 + centroid *Point +} + +func NewAreaCentroid(area float64, centroid *Point) AreaCentroid { + return AreaCentroid{ + area: area, + centroid: centroid, + } +} + +func (ac AreaCentroid) GetArea() float64 { + return ac.area +} + +func (ac AreaCentroid) GetCentroid() Point { + if ac.centroid == nil { + panic("no centroid") + } + return *ac.centroid +} diff --git a/s2/loop.go b/s2/loop.go index 4d65fd4a..a8f00ace 100644 --- a/s2/loop.go +++ b/s2/loop.go @@ -202,12 +202,102 @@ func (l *Loop) Invert() { l.initFirstLogicalVertex() } +/** + * Helper method to get area and optionally centroid. + */ +func (l *Loop) getAreaCentroid(doCentroid bool) AreaCentroid { + var centroid *Point + // Don't crash even if loop is not well-defined. + if l.NumVertices() < 3 { + return NewAreaCentroid(0, nil) + } + + // 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 = PointFromCoords(slightlyDisplaced, origin.Y, origin.Z) + case 1: + origin = PointFromCoords(origin.X, slightlyDisplaced, origin.Z) + case 2: + origin = PointFromCoords(origin.X, origin.Y, slightlyDisplaced) + } + origin = Point{origin.Normalize()} + + var areaSum float64 = 0 + centroidSum := PointFromCoords(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 = ¢roidSum + } + 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 1 //l.GetAreaCentroid(false).getArea() + 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 { @@ -457,17 +547,16 @@ func (l *Loop) ContainsPoint(p Point) bool { inside = inside != crosser.EdgeOrVertexCrossing(l.Vertex(i)) } } else { - panic("TODO") - // DataEdgeIterator it = getEdgeIterator(numVertices) - // int previousIndex = -2 - // for (it.getCandidates(origin, p); it.hasNext(); it.next()) { - // int ai = it.index() - // if (previousIndex != ai - 1) { - // crosser.restartAt(vertices[ai]) - // } - // previousIndex = ai - // inside ^= crosser.EdgeOrVertexCrossing(vertex(ai + 1)) - // } + 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 @@ -697,35 +786,39 @@ func (l *Loop) findVertex(p Point) int { * intersections and no shared vertices. */ func (l *Loop) checkEdgeCrossings(b *Loop, relation WedgeRelation) int { - // DataEdgeIterator it = getEdgeIterator(b.numVertices); + 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 (int j = 0; j < b.numVertices(); ++j) { - // S2EdgeUtil.EdgeCrosser crosser = - // new S2EdgeUtil.EdgeCrosser(b.vertex(j), b.vertex(j + 1), vertex(0)); - // int previousIndex = -2; - // for (it.getCandidates(b.vertex(j), b.vertex(j + 1)); it.hasNext(); it.next()) { - // int i = it.index(); - // if (previousIndex != i - 1) { - // crosser.restartAt(vertex(i)); - // } - // previousIndex = i; - // int crossing = crosser.robustCrossing(vertex(i + 1)); - // if (crossing < 0) { - // continue; - // } - // if (crossing > 0) { - // return -1; // There is a proper edge crossing. - // } - // if (vertex(i + 1).equals(b.vertex(j + 1))) { - // result = Math.min(result, relation.test( - // vertex(i), vertex(i + 1), vertex(i + 2), b.vertex(j), b.vertex(j + 2))); - // if (result < 0) { - // return result; - // } - // } - // } - // } + 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/point.go b/s2/point.go index d573eb73..9756fea1 100644 --- a/s2/point.go +++ b/s2/point.go @@ -241,33 +241,6 @@ 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++ - } - if RobustSign(c, o, b) != Clockwise { - sum++ - } - if RobustSign(a, o, c) == CounterClockwise { - sum++ - } - return sum >= 2 -} - // Distance returns the angle between two points. func (p Point) Distance(b Point) s1.Angle { return p.Vector.Angle(b.Vector) diff --git a/s2/s2.go b/s2/s2.go index 494b1ac2..7b6c27ca 100644 --- a/s2/s2.go +++ b/s2/s2.go @@ -2,6 +2,8 @@ package s2 import ( "math" + + "github.com/golang/geo/r2" ) const EPSILON float64 = 1e-14 @@ -51,6 +53,181 @@ func approxEqualsNumber(a, b, maxError float64) bool { return math.Abs(a-b) <= maxError } +/** + * Return the area of triangle ABC. The method used is about twice as + * expensive as Girard's formula, but it is numerically stable for both large + * and very small triangles. The points do not need to be normalized. The area + * is always positive. + * + * The triangle area is undefined if it contains two antipodal points, and + * becomes numerically unstable as the length of any edge approaches 180 + * degrees. + */ +func Area(a, b, c Point) float64 { + // This method is based on l'Huilier's theorem, + // + // tan(E/4) = sqrt(tan(s/2) tan((s-a)/2) tan((s-b)/2) tan((s-c)/2)) + // + // where E is the spherical excess of the triangle (i.e. its area), + // a, b, c, are the side lengths, and + // s is the semiperimeter (a + b + c) / 2 . + // + // The only significant source of error using l'Huilier's method is the + // cancellation error of the terms (s-a), (s-b), (s-c). This leads to a + // *relative* error of about 1e-16 * s / min(s-a, s-b, s-c). This compares + // to a relative error of about 1e-15 / E using Girard's formula, where E is + // the true area of the triangle. Girard's formula can be even worse than + // this for very small triangles, e.g. a triangle with a true area of 1e-30 + // might evaluate to 1e-5. + // + // So, we prefer l'Huilier's formula unless dmin < s * (0.1 * E), where + // dmin = min(s-a, s-b, s-c). This basically includes all triangles + // except for extremely long and skinny ones. + // + // Since we don't know E, we would like a conservative upper bound on + // the triangle area in terms of s and dmin. It's possible to show that + // E <= k1 * s * sqrt(s * dmin), where k1 = 2*sqrt(3)/Pi (about 1). + // Using this, it's easy to show that we should always use l'Huilier's + // method if dmin >= k2 * s^5, where k2 is about 1e-2. Furthermore, + // if dmin < k2 * s^5, the triangle area is at most k3 * s^4, where + // k3 is about 0.1. Since the best case error using Girard's formula + // is about 1e-15, this means that we shouldn't even consider it unless + // s >= 3e-4 or so. + + // We use volatile doubles to force the compiler to truncate all of these + // quantities to 64 bits. Otherwise it may compute a value of dmin > 0 + // simply because it chose to spill one of the intermediate values to + // memory but not one of the others. + sa := b.Angle(c.Vector).Radians() + sb := c.Angle(a.Vector).Radians() + sc := a.Angle(b.Vector).Radians() + s := 0.5 * (sa + sb + sc) + if s >= 3e-4 { + // Consider whether Girard's formula might be more accurate. + s2 := s * s + dmin := s - math.Max(sa, math.Max(sb, sc)) + if dmin < 1e-2*s*s2*s2 { + // This triangle is skinny enough to consider Girard's formula. + area := GirardArea(a, b, c) + if dmin < s*(0.1*area) { + return area + } + } + } + // Use l'Huilier's formula. + return 4 * math.Atan( + math.Sqrt( + math.Max(0.0, math.Tan(0.5*s)*math.Tan(0.5*(s-sa))*math.Tan(0.5*(s-sb))*math.Tan(0.5*(s-sc))))) +} + +/** + * Return the area of the triangle computed using Girard's formula. This is + * slightly faster than the Area() method above is not accurate for very small + * triangles. + */ +func GirardArea(a, b, c Point) float64 { + // This is equivalent to the usual Girard's formula but is slightly + // more accurate, faster to compute, and handles a == b == c without + // a special case. + + ab := a.Cross(b.Vector) + bc := b.Cross(c.Vector) + ac := a.Cross(c.Vector) + return math.Max(0.0, ab.Angle(ac).Radians()-ab.Angle(bc).Radians()+bc.Angle(ac).Radians()) +} + +/** + * Like Area(), but returns a positive value for counterclockwise triangles + * and a negative value otherwise. + */ +func SignedArea(a, b, c Point) float64 { + return Area(a, b, c) * float64(RobustCCW(a, b, c)) +} + +// About centroids: +// ---------------- +// +// There are several notions of the "centroid" of a triangle. First, there +// // is the planar centroid, which is simply the centroid of the ordinary +// (non-spherical) triangle defined by the three vertices. Second, there is +// the surface centroid, which is defined as the intersection of the three +// medians of the spherical triangle. It is possible to show that this +// point is simply the planar centroid projected to the surface of the +// sphere. Finally, there is the true centroid (mass centroid), which is +// defined as the area integral over the spherical triangle of (x,y,z) +// divided by the triangle area. This is the point that the triangle would +// rotate around if it was spinning in empty space. +// +// The best centroid for most purposes is the true centroid. Unlike the +// planar and surface centroids, the true centroid behaves linearly as +// regions are added or subtracted. That is, if you split a triangle into +// pieces and compute the average of their centroids (weighted by triangle +// area), the result equals the centroid of the original triangle. This is +// not true of the other centroids. +// +// Also note that the surface centroid may be nowhere near the intuitive +// "center" of a spherical triangle. For example, consider the triangle +// with vertices A=(1,eps,0), B=(0,0,1), C=(-1,eps,0) (a quarter-sphere). +// The surface centroid of this triangle is at S=(0, 2*eps, 1), which is +// within a distance of 2*eps of the vertex B. Note that the median from A +// (the segment connecting A to the midpoint of BC) passes through S, since +// this is the shortest path connecting the two endpoints. On the other +// hand, the true centroid is at M=(0, 0.5, 0.5), which when projected onto +// the surface is a much more reasonable interpretation of the "center" of +// this triangle. + +/** + * Return the centroid of the planar triangle ABC. This can be normalized to + * unit length to obtain the "surface centroid" of the corresponding spherical + * triangle, i.e. the intersection of the three medians. However, note that + * for large spherical triangles the surface centroid may be nowhere near the + * intuitive "center" (see example above). + */ +func PlanarCentroid(a, b, c Point) Point { + return PointFromCoords((a.X+b.X+c.X)/3.0, (a.Y+b.Y+c.Y)/3.0, (a.Z+b.Z+c.Z)/3.0) +} + +/** + * Returns the true centroid of the spherical triangle ABC multiplied by the + * signed area of spherical triangle ABC. The reasons for multiplying by the + * signed area are (1) this is the quantity that needs to be summed to compute + * the centroid of a union or difference of triangles, and (2) it's actually + * easier to calculate this way. + */ +func TrueCentroid(a, b, c Point) Point { + // I couldn't find any references for computing the true centroid of a + // spherical triangle... I have a truly marvellous demonstration of this + // formula which this margin is too narrow to contain :) + + // assert (isUnitLength(a) && isUnitLength(b) && isUnitLength(c)); + sina := b.Cross(c.Vector).Norm() + sinb := c.Cross(a.Vector).Norm() + sinc := a.Cross(b.Vector).Norm() + ra := math.Asin(sina) / sina + rb := math.Asin(sinb) / sinb + rc := math.Asin(sinc) / sinc + if sina == 0 { + ra = 1 + } + if sinb == 0 { + rb = 1 + } + if sinc == 0 { + rc = 1 + } + + // Now compute a point M such that M.X = rX * det(ABC) / 2 for X in A,B,C. + x := PointFromCoords(a.X, b.X, c.X) + y := PointFromCoords(a.Y, b.Y, c.Y) + z := PointFromCoords(a.Z, b.Z, c.Z) + r := PointFromCoords(ra, rb, rc) + return PointFromCoords( + 0.5*y.Cross(z.Vector).Dot(r.Vector), + 0.5*z.Cross(x.Vector).Dot(r.Vector), + 0.5*x.Cross(y.Vector).Dot(r.Vector), + ) +} + /** * Return true if the points A, B, C are strictly counterclockwise. Return * false if the points are clockwise or colinear (i.e. if they are all @@ -77,6 +254,253 @@ func simpleCCW(a, b, c Point) bool { return c.Cross(a.Vector).Dot(b.Vector) > 0 } +/** + * WARNING! This requires arbitrary precision arithmetic to be truly robust. + * This means that for nearly colinear AB and AC, this function may return the + * wrong answer. + * + *

+ * Like SimpleCCW(), but returns +1 if the points are counterclockwise and -1 + * if the points are clockwise. It satisfies the following conditions: + * + * (1) RobustCCW(a,b,c) == 0 if and only if a == b, b == c, or c == a (2) + * RobustCCW(b,c,a) == RobustCCW(a,b,c) for all a,b,c (3) RobustCCW(c,b,a) + * ==-RobustCCW(a,b,c) for all a,b,c + * + * In other words: + * + * (1) The result is zero if and only if two points are the same. (2) + * Rotating the order of the arguments does not affect the result. (3) + * Exchanging any two arguments inverts the result. + * + * This function is essentially like taking the sign of the determinant of + * a,b,c, except that it has additional logic to make sure that the above + * properties hold even when the three points are coplanar, and to deal with + * the limitations of floating-point arithmetic. + * + * Note: a, b and c are expected to be of unit length. Otherwise, the results + * are undefined. + */ +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: + *

    + *
  1. If orderedCCW(a,b,c,o) && orderedCCW(b,a,c,o), then a == b
  2. + *
  3. If orderedCCW(a,b,c,o) && orderedCCW(a,c,b,o), then b == c
  4. + *
  5. If orderedCCW(a,b,c,o) && orderedCCW(c,b,a,o), then a == b == c
  6. + *
  7. If a == b or b == c, then orderedCCW(a,b,c,o) is true
  8. + *
  9. Otherwise if a == c, then orderedCCW(a,b,c,o) is false
  10. + *
+ */ +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 From c727539a605c7fb600308ea402e96a31974b6209 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Mon, 4 May 2015 19:27:49 -0700 Subject: [PATCH 15/23] Added some of the Loop tests and fixed encountered bugs --- s2/edgeutil.go | 21 ++--- s2/loop.go | 22 +++--- s2/loop_test.go | 204 ++++++++++++++++++++++++++++++++++++++++++++++++ s2/point.go | 4 + 4 files changed, 231 insertions(+), 20 deletions(-) create mode 100644 s2/loop_test.go diff --git a/s2/edgeutil.go b/s2/edgeutil.go index 1e1d6eff..68e5542d 100644 --- a/s2/edgeutil.go +++ b/s2/edgeutil.go @@ -9,8 +9,9 @@ import ( type EdgeCrosser struct { // The fields below are all constant. - a Point - b Point + a Point + b Point + aCrossB Point // The fields below are updated for each vertex in the chain. @@ -22,8 +23,9 @@ type EdgeCrosser struct { func NewEdgeCrosser(a, b, c Point) *EdgeCrosser { ec := &EdgeCrosser{ - a: a, - b: b, + a: a, + b: b, + aCrossB: Point{a.Cross(b.Vector)}, } ec.RestartAt(c) return ec @@ -31,7 +33,7 @@ func NewEdgeCrosser(a, b, c Point) *EdgeCrosser { func (ec *EdgeCrosser) RestartAt(c Point) { ec.c = c - ec.acb = -int(RobustSign(ec.a, ec.b, c)) + ec.acb = -int(RobustCCWWithCross(ec.a, ec.b, ec.c, ec.aCrossB)) } /** @@ -52,7 +54,7 @@ func (ec *EdgeCrosser) RobustCrossing(d Point) int { // Recall that robustCCW is invariant with respect to rotating its // arguments, i.e. ABC has the same orientation as BDA. - bda := int(RobustSign(ec.a, ec.b, d)) + bda := int(RobustCCWWithCross(ec.a, ec.b, d, ec.aCrossB)) var result int if bda == -ec.acb && bda != 0 { @@ -81,7 +83,7 @@ func (ec *EdgeCrosser) RobustCrossing(d Point) int { */ func (ec *EdgeCrosser) EdgeOrVertexCrossing(d Point) bool { // We need to copy c since it is clobbered by robustCrossing(). - c2 := Point{ec.c.Vector} + c2 := PointFromCoordsRaw(ec.c.X, ec.c.Y, ec.c.Z) crossing := ec.RobustCrossing(d) if crossing < 0 { @@ -100,12 +102,13 @@ func (ec *EdgeCrosser) EdgeOrVertexCrossing(d Point) bool { func (ec *EdgeCrosser) robustCrossingInternal(d Point) int { // ACB and BDA have the appropriate orientations, so now we check the // triangles CBD and DAC. - cbd := -int(RobustSign(ec.c, d, ec.b)) + cCrossD := Point{ec.c.Cross(d.Vector)} + cbd := -int(RobustCCWWithCross(ec.c, d, ec.b, cCrossD)) if cbd != ec.acb { return -1 } - dac := int(RobustSign(ec.c, d, ec.a)) + dac := int(RobustCCWWithCross(ec.c, d, ec.a, cCrossD)) if dac == ec.acb { return 1 } else { diff --git a/s2/loop.go b/s2/loop.go index a8f00ace..a960382e 100644 --- a/s2/loop.go +++ b/s2/loop.go @@ -233,16 +233,16 @@ func (l *Loop) getAreaCentroid(doCentroid bool) AreaCentroid { slightlyDisplaced := origin.GetAxis(axis) + math.E*1e-10 switch axis { case 0: - origin = PointFromCoords(slightlyDisplaced, origin.Y, origin.Z) + origin = PointFromCoordsRaw(slightlyDisplaced, origin.Y, origin.Z) case 1: - origin = PointFromCoords(origin.X, slightlyDisplaced, origin.Z) + origin = PointFromCoordsRaw(origin.X, slightlyDisplaced, origin.Z) case 2: - origin = PointFromCoords(origin.X, origin.Y, slightlyDisplaced) + origin = PointFromCoordsRaw(origin.X, origin.Y, slightlyDisplaced) } origin = Point{origin.Normalize()} var areaSum float64 = 0 - centroidSum := PointFromCoords(0, 0, 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 { @@ -537,7 +537,7 @@ func (l *Loop) ContainsPoint(p Point) bool { } inside := l.originInside - origin := PointFromCoords(0, 1, 0) + origin := PointFromCoordsRaw(0, 1, 0) crosser := NewEdgeCrosser(origin, p, l.Vertex(l.NumVertices()-1)) // The s2edgeindex library is not optimized yet for long edges, @@ -608,7 +608,7 @@ func (l *Loop) IsValid() bool { // 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", i) + fmt.Printf("Vertex %d is not unit length\n", i) return false } } @@ -617,7 +617,7 @@ func (l *Loop) IsValid() bool { 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", previousVertexIndex, i) + fmt.Printf("Duplicate vertices: %d and %d\n", previousVertexIndex, i) return false } vmap[l.Vertex(i)] = i @@ -667,8 +667,8 @@ func (l *Loop) IsValid() bool { crosses = crosser.RobustCrossing(l.Vertex(b2)) > 0 previousIndex = b2 if crosses { - fmt.Printf("Edges %d and %d cross", a1, b1) - fmt.Printf("Edge locations in degrees: %s-%s and %s-%s", + 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(), @@ -743,14 +743,14 @@ func (l *Loop) initBound() { // 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(PointFromCoords(0, 0, 1)) { + 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(PointFromCoords(0, 0, -1)) { + 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 diff --git a/s2/loop_test.go b/s2/loop_test.go new file mode 100644 index 00000000..a161b42f --- /dev/null +++ b/s2/loop_test.go @@ -0,0 +1,204 @@ +package s2 + +import ( + // "runtime/debug" + "testing" +) + +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 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") + // } +} diff --git a/s2/point.go b/s2/point.go index 9756fea1..86afbbe7 100644 --- a/s2/point.go +++ b/s2/point.go @@ -84,6 +84,10 @@ 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() From ba3970ba26f7c177dad8a1828541474ca7967080 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Mon, 4 May 2015 23:01:11 -0700 Subject: [PATCH 16/23] Added remaining Loop tests and fixed more bugs --- s2/cellid.go | 2 +- s2/edgeutil.go | 4 +- s2/loop.go | 6 +- s2/loop_test.go | 297 ++++++++++++++++++++++++++++++++++++++++++++++-- s2/rect.go | 12 +- s2/s2.go | 30 ++--- 6 files changed, 318 insertions(+), 33 deletions(-) diff --git a/s2/cellid.go b/s2/cellid.go index d8aad72f..4abe276c 100644 --- a/s2/cellid.go +++ b/s2/cellid.go @@ -50,7 +50,7 @@ const ( func CellIDNone() CellID { return CellID(0) } func CellIDSentinel() CellID { return CellID(^uint64(0)) } func CellIDBegin(level int) CellID { return CellIDFromFacePosLevel(0, 0, 0).ChildBeginAtLevel(level) } -func CellIDEnd(level int) CellID { return CellIDFromFacePosLevel(5, 0, 0).ChildBeginAtLevel(level) } +func CellIDEnd(level int) CellID { return CellIDFromFacePosLevel(5, 0, 0).ChildEndAtLevel(level) } // CellIDFromFacePosLevel returns a cell given its face in the range // [0,5], the 61-bit Hilbert curve position pos within that face, and diff --git a/s2/edgeutil.go b/s2/edgeutil.go index 68e5542d..5cf19609 100644 --- a/s2/edgeutil.go +++ b/s2/edgeutil.go @@ -134,7 +134,7 @@ func (rb *RectBounder) AddPoint(b Point) { } 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(RectFromLatLngPair(rb.aLatLng, bLatLng)) + 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 @@ -142,7 +142,7 @@ func (rb *RectBounder) AddPoint(b Point) { // 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(PointFromCoords(0, 0, 1).Vector) + dir := aCrossB.Cross(PointFromCoordsRaw(0, 0, 1).Vector) da := dir.Dot(rb.a.Vector) db := dir.Dot(b.Vector) diff --git a/s2/loop.go b/s2/loop.go index a960382e..12ee4098 100644 --- a/s2/loop.go +++ b/s2/loop.go @@ -206,10 +206,10 @@ func (l *Loop) Invert() { * Helper method to get area and optionally centroid. */ func (l *Loop) getAreaCentroid(doCentroid bool) AreaCentroid { - var centroid *Point + var centroid *Point = nil // Don't crash even if loop is not well-defined. if l.NumVertices() < 3 { - return NewAreaCentroid(0, nil) + return NewAreaCentroid(0, centroid) } // The triangle area calculation becomes numerically unstable as the length @@ -268,7 +268,7 @@ func (l *Loop) getAreaCentroid(doCentroid bool) AreaCentroid { // The loop's sign() does not affect the return result and should be taken // into account by the caller. if doCentroid { - centroid = ¢roidSum + centroid = &Point{centroidSum.Vector} } return NewAreaCentroid(areaSum, centroid) } diff --git a/s2/loop_test.go b/s2/loop_test.go index a161b42f..f92ef8cb 100644 --- a/s2/loop_test.go +++ b/s2/loop_test.go @@ -1,8 +1,12 @@ package s2 import ( - // "runtime/debug" + "math" + "math/rand" "testing" + + "github.com/golang/geo/r1" + "github.com/golang/geo/s1" ) var ( @@ -70,6 +74,189 @@ 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 { @@ -192,13 +379,103 @@ func TestLoopRoundingError(t *testing.T) { } 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") - // } + 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/rect.go b/s2/rect.go index 541a0b8b..e3eef5a3 100644 --- a/s2/rect.go +++ b/s2/rect.go @@ -49,6 +49,14 @@ func RectFromLatLng(p LatLng) Rect { } } +func RectFromLatLngLoHi(lo, hi LatLng) Rect { + // assert (p1.isValid() && p2.isValid()); + return Rect{ + Lat: r1.IntervalFromPointPair(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 @@ -56,7 +64,7 @@ func RectFromLatLng(p LatLng) Rect { * S2LatLngRect(lo, hi) constructor, where the first point is always used as * the lower-left corner of the resulting rectangle. */ -func RectFromLatLngPair(p1, p2 LatLng) Rect { +func RectFromLatLngPointPair(p1, p2 LatLng) Rect { // assert (p1.isValid() && p2.isValid()); return Rect{ Lat: r1.IntervalFromPointPair(p1.Lat.Radians(), p2.Lat.Radians()), @@ -223,7 +231,7 @@ func (r Rect) CapBound() Cap { poleAngle = math.Pi/2 - r.Lat.Lo } - poleCap := CapFromCenterAngle(PointFromCoords(0, 0, poleZ), s1.Angle(poleAngle)) + 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 diff --git a/s2/s2.go b/s2/s2.go index 7b6c27ca..85069629 100644 --- a/s2/s2.go +++ b/s2/s2.go @@ -184,7 +184,7 @@ func SignedArea(a, b, c Point) float64 { * intuitive "center" (see example above). */ func PlanarCentroid(a, b, c Point) Point { - return PointFromCoords((a.X+b.X+c.X)/3.0, (a.Y+b.Y+c.Y)/3.0, (a.Z+b.Z+c.Z)/3.0) + return PointFromCoordsRaw((a.X+b.X+c.X)/3.0, (a.Y+b.Y+c.Y)/3.0, (a.Z+b.Z+c.Z)/3.0) } /** @@ -203,25 +203,25 @@ func TrueCentroid(a, b, c Point) Point { sina := b.Cross(c.Vector).Norm() sinb := c.Cross(a.Vector).Norm() sinc := a.Cross(b.Vector).Norm() - ra := math.Asin(sina) / sina - rb := math.Asin(sinb) / sinb - rc := math.Asin(sinc) / sinc - if sina == 0 { - ra = 1 + var ra float64 = 1 + var rb float64 = 1 + var rc float64 = 1 + if sina != 0 { + ra = math.Asin(sina) / sina } - if sinb == 0 { - rb = 1 + if sinb != 0 { + rb = math.Asin(sinb) / sinb } - if sinc == 0 { - rc = 1 + if sinc != 0 { + rc = math.Asin(sinc) / sinc } // Now compute a point M such that M.X = rX * det(ABC) / 2 for X in A,B,C. - x := PointFromCoords(a.X, b.X, c.X) - y := PointFromCoords(a.Y, b.Y, c.Y) - z := PointFromCoords(a.Z, b.Z, c.Z) - r := PointFromCoords(ra, rb, rc) - return PointFromCoords( + x := PointFromCoordsRaw(a.X, b.X, c.X) + y := PointFromCoordsRaw(a.Y, b.Y, c.Y) + z := PointFromCoordsRaw(a.Z, b.Z, c.Z) + r := PointFromCoordsRaw(ra, rb, rc) + return PointFromCoordsRaw( 0.5*y.Cross(z.Vector).Dot(r.Vector), 0.5*z.Cross(x.Vector).Dot(r.Vector), 0.5*x.Cross(y.Vector).Dot(r.Vector), From 4a2e47f687875995655592a81707b3934d92cfd7 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Tue, 5 May 2015 18:20:26 -0700 Subject: [PATCH 17/23] Made RegionCoverer PriorityQueue behave the same as the Java version of the s2 geometry library and added some simple tests --- s2/regioncoverer.go | 70 +++++++++++++++++---------- s2/regioncoverer_test.go | 101 ++++++++++++++++++++++++--------------- 2 files changed, 108 insertions(+), 63 deletions(-) diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index 2b24907d..acfc89bc 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -1,9 +1,5 @@ package s2 -import ( - "container/heap" -) - var ( DEFAULT_MAX_CELLS int = 8 @@ -44,8 +40,6 @@ type candidate struct { type queueEntry struct { priority int candidate *candidate - // The index is needed by update and is maintained by the heap.Interface methods. - index int // The index of the item in the heap. } func newQueueEntry(id int, candidate_ *candidate) *queueEntry { @@ -57,44 +51,71 @@ type PriorityQueue []*queueEntry func newPriorityQueue(space int) PriorityQueue { pq := PriorityQueue(make([]*queueEntry, 0, space)) - heap.Init(&pq) return pq } func (pq PriorityQueue) Len() int { return len(pq) } -func (pq PriorityQueue) Less(i, j int) bool { - // We want Pop to give us the highest, not lowest, priority so we use greater than here. - return pq[i].priority >= pq[j].priority +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) { +func (pq PriorityQueue) swap(i, j int) { pq[i], pq[j] = pq[j], pq[i] - pq[i].index = i - pq[j].index = j } func (pq *PriorityQueue) Push(x interface{}) { - n := len(*pq) item := x.(*queueEntry) - item.index = n *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) + n = len(old) item := old[n-1] - item.index = -1 // for safety *pq = old[0 : n-1] return item } -// update modifies the priority and value of an Item in the queue. -func (pq *PriorityQueue) update(item *queueEntry, candidate *candidate, priority int) { - item.candidate = candidate - item.priority = priority - heap.Fix(pq, item.index) +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 { @@ -184,7 +205,8 @@ func (rc *RegionCoverer) GetCoveringInternal(region Region) { rc.getInitialCandidates() for rc.candidateQueue.Len() > 0 && (!rc.interiorCovering || len(rc.result) < rc.maxCells) { - candidate := heap.Pop(&rc.candidateQueue).(*queueEntry).candidate + qEntry := rc.candidateQueue.Pop().(*queueEntry) + candidate := qEntry.candidate sz := len(rc.result) + candidate.numChildren if !rc.interiorCovering { sz = sz + rc.candidateQueue.Len() @@ -273,7 +295,7 @@ func (rc *RegionCoverer) addCandidate(candidate_ *candidate) { // intersecting children. Finally, we prefer cells that have the smallest // number of children that cannot be refined any further. priority := -((((int(candidate_.cell.Level()) << rc.maxChildrenShift()) + candidate_.numChildren) << rc.maxChildrenShift()) + numTerminals) - heap.Push(&rc.candidateQueue, newQueueEntry(priority, candidate_)) + rc.candidateQueue.Push(newQueueEntry(priority, candidate_)) } } diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index 437afcfa..4e20c5d0 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -1,69 +1,92 @@ package s2 import ( - "container/heap" - "fmt" "testing" ) -func TestPriorityQueue(t *testing.T) { - pq := newPriorityQueue(4) - heap.Init(&pq) - for i := -4; i < 0; i++ { - heap.Push(&pq, newQueueEntry(i, nil)) +func TestCoveringCap(t *testing.T) { + expectedResultsFromJavaLibrary := []string{ + "80c297c50b", + "80c297c50d", + "80c297c512aaaaab", + "80c297c57", + "80c297c582aaaaab", + "80c297c59d", + "80c297c59f", + "80c297c5a1", } - for i := -1; i >= -4; i-- { - entry := heap.Pop(&pq).(*queueEntry) - if entry.priority != i { - t.Errorf("expected %d, got %d", i, entry.priority) - } - } -} - -func TestCovering(t *testing.T) { coverer := NewRegionCoverer() coverer.SetMaxCells(8) - region1 := CellFromCellID(CellIDFromToken("80c297c574")) - rect1 := region1.RectBound() - fmt.Printf("%s\n", rect1.String()) - - region := region1.CapBound() - rect := region.RectBound() - fmt.Printf("%s\n", rect.String()) - - fmt.Printf("%s\n", region.String()) + cell := CellFromCellID(CellIDFromToken("80c297c574")) + cap := cell.CapBound() cells := []CellID{} - // cells := coverer.GetCoveringAsUnion(region) - coverer.GetCovering(region, &cells) + coverer.GetCovering(cap, &cells) for i, cell := range cells { - t.Errorf("cell %d: %x - %s", i, uint64(cell), cell.ToToken()) + if cell.ToToken() != expectedResultsFromJavaLibrary[i] { + t.Fatal("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) + } } } func TestCoveringPolyline(t *testing.T) { - coverer := NewRegionCoverer() - coverer.SetMaxCells(4) + expectedResultsFromJavaLibrary := []string{ + "80c2bea0418c", + "80c2bea04194", + "80c2bea041c4", + "80c2bea041eb", + "80c2bea0423", + "80c2bea0434", + "80c2bea043b", + "80c2bea043caac", + } points := []Point{ PointFromLatLng(LatLngFromDegrees(34.0909533022671600, -118.3914214745164100)), PointFromLatLng(LatLngFromDegrees(34.0906409358360560, -118.3911871165037200)), } + polyline := PolylineFromPoints(points) + + coverer := NewRegionCoverer() + coverer.SetMaxCells(8) - for _, point := range points { - fmt.Printf("point: %v\n", point.String()) + cells := []CellID{} + coverer.GetCovering(polyline, &cells) + for i, cell := range cells { + if cell.ToToken() != expectedResultsFromJavaLibrary[i] { + t.Fatal("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) + } } +} - // region := RectFromLatLng(LatLngFromDegrees(34.0909533022671600, -118.3914214745164100)) - // region = region.AddPoint(LatLngFromDegrees(34.0906409358360560, -118.3911871165037200)) +func TestCoveringLoop(t *testing.T) { + expectedResultsFromJavaLibrary := []string{ + "80c2c7b46204", + "80c2c7b4620c", + "80c2c7b4899", + "80c2c7b489c4", + "80c2c7b489dc", + "80c2c7b489f", + "80c2c7b48a04", + "80c2c7b48a1d", + } + vertices := []Point{ + PointFromLatLng(LatLngFromDegrees(34.0487325747361496, -118.2554703578353070)), + PointFromLatLng(LatLngFromDegrees(34.0486331233724258, -118.2555538415908671)), + PointFromLatLng(LatLngFromDegrees(34.0486803488948837, -118.2556309551000595)), + PointFromLatLng(LatLngFromDegrees(34.0487739664704705, -118.2555437833070613)), + } + loop := LoopFromPoints(vertices) + loop.Normalize() - polyline := PolylineFromPoints(points) - // region := polyline.RectBound() - // fmt.Printf("%s\n", region.String()) + coverer := NewRegionCoverer() + coverer.SetMaxCells(8) cells := []CellID{} - coverer.GetCovering(polyline, &cells) + coverer.GetCovering(loop, &cells) for i, cell := range cells { - t.Errorf("cell %d: %x - %s", i, uint64(cell), cell.ToToken()) + if cell.ToToken() != expectedResultsFromJavaLibrary[i] { + t.Fatal("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) + } } } From cb1064200c9713173d5e0830c4b5eaa908fecf45 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Tue, 5 May 2015 18:25:18 -0700 Subject: [PATCH 18/23] Typo in test reporting --- s2/regioncoverer_test.go | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index 4e20c5d0..fa755a95 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -24,7 +24,7 @@ func TestCoveringCap(t *testing.T) { coverer.GetCovering(cap, &cells) for i, cell := range cells { if cell.ToToken() != expectedResultsFromJavaLibrary[i] { - t.Fatal("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) + t.Fatalf("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) } } } @@ -54,7 +54,7 @@ func TestCoveringPolyline(t *testing.T) { coverer.GetCovering(polyline, &cells) for i, cell := range cells { if cell.ToToken() != expectedResultsFromJavaLibrary[i] { - t.Fatal("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) + t.Fatalf("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) } } } @@ -86,7 +86,7 @@ func TestCoveringLoop(t *testing.T) { coverer.GetCovering(loop, &cells) for i, cell := range cells { if cell.ToToken() != expectedResultsFromJavaLibrary[i] { - t.Fatal("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) + t.Fatalf("TestCoveringLoop result %d got %s expected %s", i, cell.ToToken(), expectedResultsFromJavaLibrary[i]) } } } From 9029e33d95aa2df863ce97593697ac457a8e33fb Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Tue, 5 May 2015 19:05:41 -0700 Subject: [PATCH 19/23] Minor adjustments in form --- s2/cap.go | 13 +++++++------ s2/latlng.go | 2 +- s2/loop_test.go | 2 +- s2/point.go | 2 +- s2/rect.go | 2 +- 5 files changed, 11 insertions(+), 10 deletions(-) diff --git a/s2/cap.go b/s2/cap.go index 019162b1..169caf88 100644 --- a/s2/cap.go +++ b/s2/cap.go @@ -249,12 +249,13 @@ func (c Cap) RectBound() Rect { // ApproxEqual reports if this caps' center and height are within // a reasonable epsilon from the other cap. func (c Cap) ApproxEqual(other Cap) bool { - return c.center.ApproxEquals(other.center, EPSILON) && - math.Abs(c.height-other.height) <= EPSILON || - c.IsEmpty() && other.height <= EPSILON || - other.IsEmpty() && c.height <= EPSILON || - c.IsFull() && other.height >= 2-EPSILON || - other.IsFull() && c.height >= 2-EPSILON + const epsilon = 1e-14 + return c.center.ApproxEquals(other.center, epsilon) && + math.Abs(c.height-other.height) <= epsilon || + c.IsEmpty() && other.height <= epsilon || + other.IsEmpty() && c.height <= epsilon || + c.IsFull() && other.height >= 2-epsilon || + other.IsFull() && c.height >= 2-epsilon } // AddPoint increases the cap if necessary to include the given point. If this cap is empty, diff --git a/s2/latlng.go b/s2/latlng.go index 8ccb636c..bf0c9021 100644 --- a/s2/latlng.go +++ b/s2/latlng.go @@ -57,7 +57,7 @@ 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)))) } func latitude(p Point) s1.Angle { diff --git a/s2/loop_test.go b/s2/loop_test.go index f92ef8cb..c2a7f372 100644 --- a/s2/loop_test.go +++ b/s2/loop_test.go @@ -98,7 +98,7 @@ func TestLoopBounds(t *testing.T) { 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 { + if math.Abs(s1.Angle(arctic80.RectBound().Lat.Hi).Radians()-LatLngFromPoint(mid).Lat.Radians()) > EPSILON { t.Fatal("") } arctic80.Invert() diff --git a/s2/point.go b/s2/point.go index 86afbbe7..44f4a23f 100644 --- a/s2/point.go +++ b/s2/point.go @@ -252,7 +252,7 @@ func (p Point) Distance(b Point) s1.Angle { // 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) <= epsilon + return p.Vector.Angle(other.Vector).Radians() <= maxError } func (p Point) Equals(other Point) bool { diff --git a/s2/rect.go b/s2/rect.go index e3eef5a3..bad52e2d 100644 --- a/s2/rect.go +++ b/s2/rect.go @@ -81,7 +81,7 @@ func RectFromLatLngPointPair(p1, p2 LatLng) Rect { // Examples of clamping (in degrees): // center=(80,170), size=(40,60) -> lat=[60,90], lng=[140,-160] // center=(10,40), size=(210,400) -> lat=[-90,90], lng=[-180,180] -// center=(-90,180), size=(20,50) -> lat=[-90,-80], lng=[1 +// center=(-90,180), size=(20,50) -> lat=[-90,-80], lng=[155,-155] func RectFromCenterSize(center, size LatLng) Rect { half := LatLng{size.Lat / 2, size.Lng / 2} return RectFromLatLng(center).expanded(half) From be3ae28d8ee096a9bc36ce36b78a54a480d66347 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Tue, 5 May 2015 19:12:33 -0700 Subject: [PATCH 20/23] Reverting a comment removal and code move --- s2/latlng.go | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/s2/latlng.go b/s2/latlng.go index bf0c9021..b34efae8 100644 --- a/s2/latlng.go +++ b/s2/latlng.go @@ -33,11 +33,6 @@ func LatLngFromDegrees(lat, lng float64) LatLng { return LatLng{s1.Angle(lat) * s1.Degree, s1.Angle(lng) * s1.Degree} } -// LatLngFromPoint returns an LatLng for a given Point. -func LatLngFromPoint(p Point) LatLng { - return LatLng{latitude(p), longitude(p)} -} - // IsValid returns true iff the LatLng is normalized, with Lat ∈ [-π/2,π/2] and Lng ∈ [-π,π]. func (ll LatLng) IsValid() bool { return math.Abs(ll.Lat.Radians()) <= math.Pi/2 && math.Abs(ll.Lng.Radians()) <= math.Pi @@ -60,6 +55,9 @@ func (ll LatLng) Distance(ll2 LatLng) s1.Angle { 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))) } @@ -68,5 +66,10 @@ func longitude(p Point) s1.Angle { return s1.Angle(math.Atan2(p.Y, p.X)) } +// LatLngFromPoint returns an LatLng for a given Point. +func LatLngFromPoint(p Point) LatLng { + return LatLng{latitude(p), longitude(p)} +} + // BUG(dsymonds): The major differences from the C++ version are: // - normalization From 415810e487912a46dda25490073a63445471225f Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Wed, 6 May 2015 17:51:00 -0700 Subject: [PATCH 21/23] Added missing RegionCoverer functions, but untested --- s2/regioncoverer.go | 60 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/s2/regioncoverer.go b/s2/regioncoverer.go index acfc89bc..45b48b77 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -179,6 +179,18 @@ func (rc *RegionCoverer) GetCovering(region Region, covering *[]CellID) { 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 @@ -194,6 +206,25 @@ func (rc *RegionCoverer) GetCoveringAsUnion(region Region) *CellUnion { 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 (rc *RegionCoverer) GetSimpleCovering(region Region, start Point, level int, output *[]CellID) { + rc.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 { @@ -356,3 +387,32 @@ func (rc *RegionCoverer) newCandidate(cell Cell) *candidate { rc.candidatesCreatedCounter++ return candidate_ } + +/** + * Given a region and a starting cell, return the set of all the + * edge-connected cells at the same level that intersect "region". The output + * cells are returned in arbitrary order. + */ +func (rc *RegionCoverer) floodFill(region Region, start CellID, output *[]CellID) { + all := make(map[CellID]bool) + frontier := []CellID{} + *output = []CellID{} + all[start] = true + frontier = append(frontier, start) + for len(frontier) > 0 { + id := frontier[len(frontier)-1] + frontier = frontier[0 : len(frontier)-1] + if !region.IntersectsCell(CellFromCellID(id)) { + continue + } + *output = append(*output, id) + + neighbors := id.EdgeNeighbors() + for _, nbr := range neighbors { + if _, hasNbr := all[nbr]; !hasNbr { + frontier = append(frontier, nbr) + all[nbr] = true + } + } + } +} From ba36ea30975cc322571895882af869294a828372 Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Wed, 6 May 2015 22:37:41 -0700 Subject: [PATCH 22/23] Added more RegionCoverer tests which brought in more code into other parts as well --- s2/cell.go | 16 +++++ s2/cellunion.go | 30 ++++++++++ s2/edgeindex.go | 4 +- s2/projections.go | 123 +++++++++++++++++++++++++++++++++++---- s2/regioncoverer.go | 6 +- s2/regioncoverer_test.go | 106 +++++++++++++++++++++++++++++++++ s2/s2.go | 2 +- s2/s2_test.go | 8 +++ 8 files changed, 278 insertions(+), 17 deletions(-) diff --git a/s2/cell.go b/s2/cell.go index 3bb98067..c46272e0 100644 --- a/s2/cell.go +++ b/s2/cell.go @@ -108,6 +108,22 @@ func (c Cell) EdgeRaw(k int) Point { } } +/** + * Return the average area for cells at the given level. + */ +func AverageArea(level int) float64 { + return S2_PROJECTION.AVG_AREA().GetValue(level) +} + +/** + * Return the average area of cells at this level. This is accurate to within + * a factor of 1.7 (for S2_QUADRATIC_PROJECTION) and is extremely cheap to + * compute. + */ +func (c Cell) AverageArea() float64 { + return AverageArea(int(c.level)) +} + // ExactArea return the area of this cell as accurately as possible. func (c Cell) ExactArea() float64 { v0, v1, v2, v3 := c.Vertex(0), c.Vertex(1), c.Vertex(2), c.Vertex(3) diff --git a/s2/cellunion.go b/s2/cellunion.go index f770ce1b..c8a093b9 100644 --- a/s2/cellunion.go +++ b/s2/cellunion.go @@ -25,6 +25,13 @@ import "sort" // nor the four sibling CellIDs that are children of a single higher level CellID. type CellUnion []CellID +func CellUnionFromCellIDs(ids []CellID) *CellUnion { + union := &CellUnion{} + *union = append(*union, ids...) + union.Normalize() + return union +} + func CellUnionFromArrayAndSwap(ids *[]CellID) *CellUnion { union := &CellUnion{} *union = append(*union, *ids...) @@ -116,6 +123,29 @@ func (cu *CellUnion) Normalize() { *cu = output } +/** + * Return true if the cell union contains the given cell id. Containment is + * defined with respect to regions, e.g. a cell contains its 4 children. This + * is a fast operation (logarithmic in the size of the cell union). + */ +func (cu *CellUnion) Contains(id CellID) bool { + // This function requires that Normalize has been called first. + // + // This is an exact test. Each cell occupies a linear span of the S2 + // space-filling curve, and the cell id is simply the position at the center + // of this span. The cell union ids are sorted in increasing order along + // the space-filling curve. So we simply find the pair of cell ids that + // surround the given cell id (using binary search). There is containment + // if and only if one of these two cell ids contains this cell. + + i := sort.Search(len(*cu), func(i int) bool { return id < (*cu)[i] }) + + if i != len(*cu) && (*cu)[i].RangeMin() <= id { + return true + } + return i != 0 && (*cu)[i-1].RangeMax() >= id +} + // Intersects reports whether this cell union intersects the given cell ID. // // This method assumes that the CellUnion has been normalized. diff --git a/s2/edgeindex.go b/s2/edgeindex.go index e07ed190..f8751c83 100644 --- a/s2/edgeindex.go +++ b/s2/edgeindex.go @@ -347,8 +347,8 @@ func (e *EdgeIndex) getEdges(cell1, cell2 uint64) []int { // not valid edge indices, we will always get -N-1, so we immediately // convert to N. return []int{ - -1 - e.binarySearch(cell1, -(int(^uint(0)>>1)-1)), - -1 - e.binarySearch(cell2, int(^uint(0)>>1)), + -1 - e.binarySearch(cell1, math.MinInt32), + -1 - e.binarySearch(cell2, math.MaxInt32), } } diff --git a/s2/projections.go b/s2/projections.go index db779a18..45d6e29b 100644 --- a/s2/projections.go +++ b/s2/projections.go @@ -4,28 +4,82 @@ import ( "math" ) -// type Projections int - -// const ( -// S2_LINEAR_PROJECTION Projections = iota -// S2_TAN_PROJECTION -// S2_QUADRATIC_PROJECTION -// ) - -// const S2_PROJECTION Projections = S2_QUADRATIC_PROJECTION - 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_QUADRATIC_PROJECTION struct{} +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{ @@ -51,3 +105,50 @@ func (m S2_QUADRATIC_PROJECTION) AVG_WIDTH() Metric { 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/regioncoverer.go b/s2/regioncoverer.go index 45b48b77..447cda78 100644 --- a/s2/regioncoverer.go +++ b/s2/regioncoverer.go @@ -221,8 +221,8 @@ func (rc *RegionCoverer) GetInteriorCoveringAsUnion(region Region) *CellUnion { * Given a connected region and a starting point, return a set of cells at the * given level that cover the region. */ -func (rc *RegionCoverer) GetSimpleCovering(region Region, start Point, level int, output *[]CellID) { - rc.floodFill(region, CellFromPoint(start).Id().Parent(level), output) +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. */ @@ -393,7 +393,7 @@ func (rc *RegionCoverer) newCandidate(cell Cell) *candidate { * edge-connected cells at the same level that intersect "region". The output * cells are returned in arbitrary order. */ -func (rc *RegionCoverer) floodFill(region Region, start CellID, output *[]CellID) { +func floodFill(region Region, start CellID, output *[]CellID) { all := make(map[CellID]bool) frontier := []CellID{} *output = []CellID{} diff --git a/s2/regioncoverer_test.go b/s2/regioncoverer_test.go index fa755a95..096a4d29 100644 --- a/s2/regioncoverer_test.go +++ b/s2/regioncoverer_test.go @@ -1,6 +1,8 @@ package s2 import ( + "math" + "math/rand" "testing" ) @@ -90,3 +92,107 @@ func TestCoveringLoop(t *testing.T) { } } } + +func random(n int32) int32 { + if n == 0 { + return 0 + } + return rand.Int31n(n) +} + +/** + * Checks that "covering" completely covers the given region. If "check_tight" + * is true, also checks that it does not contain any cells that do not + * intersect the given region. ("id" is only used internally.) + */ +func checkCoveringRegion(region Region, covering CellUnion, checkTight bool, id CellID) { + if !id.IsValid() { + for face := 0; face < 6; face++ { + checkCoveringRegion(region, covering, checkTight, CellIDFromFacePosLevel(face, 0, 0)) + } + return + } + + if !region.IntersectsCell(CellFromCellID(id)) { + // If region does not intersect id, then neither should the covering. + if checkTight { + if !(!covering.Intersects(id)) { + panic("") + } + } + + } else if !covering.Contains(id) { + // The region may intersect id, but we can't assert that the covering + // intersects id because we may discover that the region does not actually + // intersect upon further subdivision. (MayIntersect is not exact.) + if region.ContainsCell(CellFromCellID(id)) { + panic("") + } + if id.IsLeaf() { + panic("") + } + end := id.ChildEnd() + for child := id.ChildBegin(); child != end; child = child.Next() { + checkCoveringRegion(region, covering, checkTight, child) + } + } +} + +func checkCovering(coverer *RegionCoverer, region Region, covering []CellID, interior bool) { + // Keep track of how many cells have the same coverer.min_level() ancestor. + minLevelCells := make(map[CellID]int) + for i := 0; i < len(covering); i++ { + level := covering[i].Level() + if !(level >= coverer.MinLevel()) { + panic("") + } + if !(level <= coverer.MaxLevel()) { + panic("") + } + if (level-coverer.MinLevel())%coverer.LevelMod() != 0 { + panic("") + } + + key := covering[i].Parent(coverer.MinLevel()) + if i, ok := minLevelCells[key]; !ok { + minLevelCells[key] = 1 + } else { + minLevelCells[key] = i + 1 + } + } + if len(covering) > coverer.MaxCells() { + // If the covering has more than the requested number of cells, then check + // that the cell count cannot be reduced by using the parent of some cell. + for _, i := range minLevelCells { + if i != 1 { + panic("") + } + } + } + + if interior { + for i := 0; i < len(covering); i++ { + if !(region.ContainsCell(CellFromCellID(covering[i]))) { + panic("") + } + } + } else { + cellUnion := CellUnionFromCellIDs(covering) + checkCoveringRegion(region, *cellUnion, true, CellIDNone()) + } +} + +func TestSimpleCoverings(t *testing.T) { + coverer := NewRegionCoverer() + coverer.SetMaxCells(math.MaxInt32) + for i := 0; i < 1000; i++ { + level := int(random(MAX_LEVEL + 1)) + coverer.SetMinLevel(level) + coverer.SetMaxLevel(level) + maxArea := math.Min(4*math.Pi, 1000*AverageArea(level)) + cap := randomCap(0.1*AverageArea(MAX_LEVEL), maxArea) + covering := []CellID{} + GetSimpleCovering(cap, cap.Center(), level, &covering) + checkCovering(coverer, cap, covering, false) + } +} diff --git a/s2/s2.go b/s2/s2.go index 85069629..b023a571 100644 --- a/s2/s2.go +++ b/s2/s2.go @@ -518,7 +518,7 @@ 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 math.Pow(m.deriv, float64(int(m.dim)*(1-level))) + return m.deriv * math.Pow(2, float64(int(m.dim)*(1-level))) } /** diff --git a/s2/s2_test.go b/s2/s2_test.go index e92005bc..3e2e2719 100644 --- a/s2/s2_test.go +++ b/s2/s2_test.go @@ -76,6 +76,14 @@ func randomPoint() Point { randomUniformDouble(-1, 1), randomUniformDouble(-1, 1)).Normalize()} } +func randomCap(minArea, maxArea float64) Cap { + capArea := maxArea * math.Pow(minArea/maxArea, rand.Float64()) + if capArea < minArea || capArea > maxArea { + panic(fmt.Sprintf("should not happen %f %f %f %f\n", capArea, minArea, capArea, maxArea)) + } + return CapFromCenterArea(randomPoint(), capArea) +} + // parsePoint returns an Point from the latitude-longitude coordinate in degrees // in the given string, or the origin if the string was invalid. // e.g., "-20:150" From 40ad1ac887dc481d6219fded15de842e9d10837a Mon Sep 17 00:00:00 2001 From: Jeffrey Exterkate Date: Tue, 19 May 2015 10:08:30 -0700 Subject: [PATCH 23/23] Fixed s2.RectFromLatLngLoHi --- r1/interval.go | 6 ++++++ s2/rect.go | 2 +- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/r1/interval.go b/r1/interval.go index 256effbe..b00bd2e9 100644 --- a/r1/interval.go +++ b/r1/interval.go @@ -34,6 +34,12 @@ func EmptyInterval() Interval { return Interval{1, 0} } // IntervalFromPoint returns an interval representing a single point. func IntervalFromPoint(p float64) Interval { return Interval{p, p} } +// Convenience method to construct an interval from two points that are +// already ordered from lo to hi +func IntervalFromEndpoints(p1, p2 float64) Interval { + return Interval{p1, p2} +} + // Convenience method to construct the minimal interval containing the two // given points. This is equivalent to starting with an empty interval and // calling AddPoint() twice, but it is more efficient. diff --git a/s2/rect.go b/s2/rect.go index bad52e2d..5d04744b 100644 --- a/s2/rect.go +++ b/s2/rect.go @@ -52,7 +52,7 @@ func RectFromLatLng(p LatLng) Rect { func RectFromLatLngLoHi(lo, hi LatLng) Rect { // assert (p1.isValid() && p2.isValid()); return Rect{ - Lat: r1.IntervalFromPointPair(lo.Lat.Radians(), hi.Lat.Radians()), + Lat: r1.IntervalFromEndpoints(lo.Lat.Radians(), hi.Lat.Radians()), Lng: s1.IntervalFromEndpoints(lo.Lng.Radians(), hi.Lng.Radians()), } }