...

Source file src/github.com/golang/geo/s2/shapeindex.go

Documentation: github.com/golang/geo/s2

     1  // Copyright 2016 Google Inc. All rights reserved.
     2  //
     3  // Licensed under the Apache License, Version 2.0 (the "License");
     4  // you may not use this file except in compliance with the License.
     5  // You may obtain a copy of the License at
     6  //
     7  //     http://www.apache.org/licenses/LICENSE-2.0
     8  //
     9  // Unless required by applicable law or agreed to in writing, software
    10  // distributed under the License is distributed on an "AS IS" BASIS,
    11  // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    12  // See the License for the specific language governing permissions and
    13  // limitations under the License.
    14  
    15  package s2
    16  
    17  import (
    18  	"math"
    19  	"sort"
    20  	"sync"
    21  	"sync/atomic"
    22  
    23  	"github.com/golang/geo/r1"
    24  	"github.com/golang/geo/r2"
    25  )
    26  
    27  // CellRelation describes the possible relationships between a target cell
    28  // and the cells of the ShapeIndex. If the target is an index cell or is
    29  // contained by an index cell, it is Indexed. If the target is subdivided
    30  // into one or more index cells, it is Subdivided. Otherwise it is Disjoint.
    31  type CellRelation int
    32  
    33  // The possible CellRelations for a ShapeIndex.
    34  const (
    35  	Indexed CellRelation = iota
    36  	Subdivided
    37  	Disjoint
    38  )
    39  
    40  const (
    41  	// cellPadding defines the total error when clipping an edge which comes
    42  	// from two sources:
    43  	// (1) Clipping the original spherical edge to a cube face (the face edge).
    44  	//     The maximum error in this step is faceClipErrorUVCoord.
    45  	// (2) Clipping the face edge to the u- or v-coordinate of a cell boundary.
    46  	//     The maximum error in this step is edgeClipErrorUVCoord.
    47  	// Finally, since we encounter the same errors when clipping query edges, we
    48  	// double the total error so that we only need to pad edges during indexing
    49  	// and not at query time.
    50  	cellPadding = 2.0 * (faceClipErrorUVCoord + edgeClipErrorUVCoord)
    51  
    52  	// cellSizeToLongEdgeRatio defines the cell size relative to the length of an
    53  	// edge at which it is first considered to be long. Long edges do not
    54  	// contribute toward the decision to subdivide a cell further. For example,
    55  	// a value of 2.0 means that the cell must be at least twice the size of the
    56  	// edge in order for that edge to be counted. There are two reasons for not
    57  	// counting long edges: (1) such edges typically need to be propagated to
    58  	// several children, which increases time and memory costs without much benefit,
    59  	// and (2) in pathological cases, many long edges close together could force
    60  	// subdivision to continue all the way to the leaf cell level.
    61  	cellSizeToLongEdgeRatio = 1.0
    62  )
    63  
    64  // clippedShape represents the part of a shape that intersects a Cell.
    65  // It consists of the set of edge IDs that intersect that cell and a boolean
    66  // indicating whether the center of the cell is inside the shape (for shapes
    67  // that have an interior).
    68  //
    69  // Note that the edges themselves are not clipped; we always use the original
    70  // edges for intersection tests so that the results will be the same as the
    71  // original shape.
    72  type clippedShape struct {
    73  	// shapeID is the index of the shape this clipped shape is a part of.
    74  	shapeID int32
    75  
    76  	// containsCenter indicates if the center of the CellID this shape has been
    77  	// clipped to falls inside this shape. This is false for shapes that do not
    78  	// have an interior.
    79  	containsCenter bool
    80  
    81  	// edges is the ordered set of ShapeIndex original edge IDs. Edges
    82  	// are stored in increasing order of edge ID.
    83  	edges []int
    84  }
    85  
    86  // newClippedShape returns a new clipped shape for the given shapeID and number of expected edges.
    87  func newClippedShape(id int32, numEdges int) *clippedShape {
    88  	return &clippedShape{
    89  		shapeID: id,
    90  		edges:   make([]int, numEdges),
    91  	}
    92  }
    93  
    94  // numEdges returns the number of edges that intersect the CellID of the Cell this was clipped to.
    95  func (c *clippedShape) numEdges() int {
    96  	return len(c.edges)
    97  }
    98  
    99  // containsEdge reports if this clipped shape contains the given edge ID.
   100  func (c *clippedShape) containsEdge(id int) bool {
   101  	// Linear search is fast because the number of edges per shape is typically
   102  	// very small (less than 10).
   103  	for _, e := range c.edges {
   104  		if e == id {
   105  			return true
   106  		}
   107  	}
   108  	return false
   109  }
   110  
   111  // ShapeIndexCell stores the index contents for a particular CellID.
   112  type ShapeIndexCell struct {
   113  	shapes []*clippedShape
   114  }
   115  
   116  // NewShapeIndexCell creates a new cell that is sized to hold the given number of shapes.
   117  func NewShapeIndexCell(numShapes int) *ShapeIndexCell {
   118  	return &ShapeIndexCell{
   119  		shapes: make([]*clippedShape, numShapes),
   120  	}
   121  }
   122  
   123  // numEdges reports the total number of edges in all clipped shapes in this cell.
   124  func (s *ShapeIndexCell) numEdges() int {
   125  	var e int
   126  	for _, cs := range s.shapes {
   127  		e += cs.numEdges()
   128  	}
   129  	return e
   130  }
   131  
   132  // add adds the given clipped shape to this index cell.
   133  func (s *ShapeIndexCell) add(c *clippedShape) {
   134  	// C++ uses a set, so it's ordered and unique. We don't currently catch
   135  	// the case when a duplicate value is added.
   136  	s.shapes = append(s.shapes, c)
   137  }
   138  
   139  // findByShapeID returns the clipped shape that contains the given shapeID,
   140  // or nil if none of the clipped shapes contain it.
   141  func (s *ShapeIndexCell) findByShapeID(shapeID int32) *clippedShape {
   142  	// Linear search is fine because the number of shapes per cell is typically
   143  	// very small (most often 1), and is large only for pathological inputs
   144  	// (e.g. very deeply nested loops).
   145  	for _, clipped := range s.shapes {
   146  		if clipped.shapeID == shapeID {
   147  			return clipped
   148  		}
   149  	}
   150  	return nil
   151  }
   152  
   153  // faceEdge and clippedEdge store temporary edge data while the index is being
   154  // updated.
   155  //
   156  // While it would be possible to combine all the edge information into one
   157  // structure, there are two good reasons for separating it:
   158  //
   159  //  - Memory usage. Separating the two means that we only need to
   160  //    store one copy of the per-face data no matter how many times an edge is
   161  //    subdivided, and it also lets us delay computing bounding boxes until
   162  //    they are needed for processing each face (when the dataset spans
   163  //    multiple faces).
   164  //
   165  //  - Performance. UpdateEdges is significantly faster on large polygons when
   166  //    the data is separated, because it often only needs to access the data in
   167  //    clippedEdge and this data is cached more successfully.
   168  
   169  // faceEdge represents an edge that has been projected onto a given face,
   170  type faceEdge struct {
   171  	shapeID     int32    // The ID of shape that this edge belongs to
   172  	edgeID      int      // Edge ID within that shape
   173  	MaxLevel    int      // Not desirable to subdivide this edge beyond this level
   174  	hasInterior bool     // Belongs to a shape that has a dimension of 2
   175  	a, b        r2.Point // The edge endpoints, clipped to a given face
   176  	edge        Edge     // The original edge.
   177  }
   178  
   179  // clippedEdge represents the portion of that edge that has been clipped to a given Cell.
   180  type clippedEdge struct {
   181  	faceEdge *faceEdge // The original unclipped edge
   182  	bound    r2.Rect   // Bounding box for the clipped portion
   183  }
   184  
   185  // ShapeIndexIteratorPos defines the set of possible iterator starting positions. By
   186  // default iterators are unpositioned, since this avoids an extra seek in this
   187  // situation where one of the seek methods (such as Locate) is immediately called.
   188  type ShapeIndexIteratorPos int
   189  
   190  const (
   191  	// IteratorBegin specifies the iterator should be positioned at the beginning of the index.
   192  	IteratorBegin ShapeIndexIteratorPos = iota
   193  	// IteratorEnd specifies the iterator should be positioned at the end of the index.
   194  	IteratorEnd
   195  )
   196  
   197  // ShapeIndexIterator is an iterator that provides low-level access to
   198  // the cells of the index. Cells are returned in increasing order of CellID.
   199  //
   200  //	for it := index.Iterator(); !it.Done(); it.Next() {
   201  //	  fmt.Print(it.CellID())
   202  //	}
   203  type ShapeIndexIterator struct {
   204  	index    *ShapeIndex
   205  	position int
   206  	id       CellID
   207  	cell     *ShapeIndexCell
   208  }
   209  
   210  // NewShapeIndexIterator creates a new iterator for the given index. If a starting
   211  // position is specified, the iterator is positioned at the given spot.
   212  func NewShapeIndexIterator(index *ShapeIndex, pos ...ShapeIndexIteratorPos) *ShapeIndexIterator {
   213  	s := &ShapeIndexIterator{
   214  		index: index,
   215  	}
   216  
   217  	if len(pos) > 0 {
   218  		if len(pos) > 1 {
   219  			panic("too many ShapeIndexIteratorPos arguments")
   220  		}
   221  		switch pos[0] {
   222  		case IteratorBegin:
   223  			s.Begin()
   224  		case IteratorEnd:
   225  			s.End()
   226  		default:
   227  			panic("unknown ShapeIndexIteratorPos value")
   228  		}
   229  	}
   230  
   231  	return s
   232  }
   233  
   234  func (s *ShapeIndexIterator) clone() *ShapeIndexIterator {
   235  	return &ShapeIndexIterator{
   236  		index:    s.index,
   237  		position: s.position,
   238  		id:       s.id,
   239  		cell:     s.cell,
   240  	}
   241  }
   242  
   243  // CellID returns the CellID of the current index cell.
   244  // If s.Done() is true, a value larger than any valid CellID is returned.
   245  func (s *ShapeIndexIterator) CellID() CellID {
   246  	return s.id
   247  }
   248  
   249  // IndexCell returns the current index cell.
   250  func (s *ShapeIndexIterator) IndexCell() *ShapeIndexCell {
   251  	// TODO(roberts): C++ has this call a virtual method to allow subclasses
   252  	// of ShapeIndexIterator to do other work before returning the cell. Do
   253  	// we need such a thing?
   254  	return s.cell
   255  }
   256  
   257  // Center returns the Point at the center of the current position of the iterator.
   258  func (s *ShapeIndexIterator) Center() Point {
   259  	return s.CellID().Point()
   260  }
   261  
   262  // Begin positions the iterator at the beginning of the index.
   263  func (s *ShapeIndexIterator) Begin() {
   264  	if !s.index.IsFresh() {
   265  		s.index.maybeApplyUpdates()
   266  	}
   267  	s.position = 0
   268  	s.refresh()
   269  }
   270  
   271  // Next positions the iterator at the next index cell.
   272  func (s *ShapeIndexIterator) Next() {
   273  	s.position++
   274  	s.refresh()
   275  }
   276  
   277  // Prev advances the iterator to the previous cell in the index and returns true to
   278  // indicate it was not yet at the beginning of the index. If the iterator is at the
   279  // first cell the call does nothing and returns false.
   280  func (s *ShapeIndexIterator) Prev() bool {
   281  	if s.position <= 0 {
   282  		return false
   283  	}
   284  
   285  	s.position--
   286  	s.refresh()
   287  	return true
   288  }
   289  
   290  // End positions the iterator at the end of the index.
   291  func (s *ShapeIndexIterator) End() {
   292  	s.position = len(s.index.cells)
   293  	s.refresh()
   294  }
   295  
   296  // Done reports if the iterator is positioned at or after the last index cell.
   297  func (s *ShapeIndexIterator) Done() bool {
   298  	return s.id == SentinelCellID
   299  }
   300  
   301  // refresh updates the stored internal iterator values.
   302  func (s *ShapeIndexIterator) refresh() {
   303  	if s.position < len(s.index.cells) {
   304  		s.id = s.index.cells[s.position]
   305  		s.cell = s.index.cellMap[s.CellID()]
   306  	} else {
   307  		s.id = SentinelCellID
   308  		s.cell = nil
   309  	}
   310  }
   311  
   312  // seek positions the iterator at the first cell whose ID >= target, or at the
   313  // end of the index if no such cell exists.
   314  func (s *ShapeIndexIterator) seek(target CellID) {
   315  	s.position = sort.Search(len(s.index.cells), func(i int) bool {
   316  		return s.index.cells[i] >= target
   317  	})
   318  	s.refresh()
   319  }
   320  
   321  // LocatePoint positions the iterator at the cell that contains the given Point.
   322  // If no such cell exists, the iterator position is unspecified, and false is returned.
   323  // The cell at the matched position is guaranteed to contain all edges that might
   324  // intersect the line segment between target and the cell's center.
   325  func (s *ShapeIndexIterator) LocatePoint(p Point) bool {
   326  	// Let I = cellMap.LowerBound(T), where T is the leaf cell containing
   327  	// point P. Then if T is contained by an index cell, then the
   328  	// containing cell is either I or I'. We test for containment by comparing
   329  	// the ranges of leaf cells spanned by T, I, and I'.
   330  	target := cellIDFromPoint(p)
   331  	s.seek(target)
   332  	if !s.Done() && s.CellID().RangeMin() <= target {
   333  		return true
   334  	}
   335  
   336  	if s.Prev() && s.CellID().RangeMax() >= target {
   337  		return true
   338  	}
   339  	return false
   340  }
   341  
   342  // LocateCellID attempts to position the iterator at the first matching index cell
   343  // in the index that has some relation to the given CellID. Let T be the target CellID.
   344  // If T is contained by (or equal to) some index cell I, then the iterator is positioned
   345  // at I and returns Indexed. Otherwise if T contains one or more (smaller) index cells,
   346  // then the iterator is positioned at the first such cell I and return Subdivided.
   347  // Otherwise Disjoint is returned and the iterator position is undefined.
   348  func (s *ShapeIndexIterator) LocateCellID(target CellID) CellRelation {
   349  	// Let T be the target, let I = cellMap.LowerBound(T.RangeMin()), and
   350  	// let I' be the predecessor of I. If T contains any index cells, then T
   351  	// contains I. Similarly, if T is contained by an index cell, then the
   352  	// containing cell is either I or I'. We test for containment by comparing
   353  	// the ranges of leaf cells spanned by T, I, and I'.
   354  	s.seek(target.RangeMin())
   355  	if !s.Done() {
   356  		if s.CellID() >= target && s.CellID().RangeMin() <= target {
   357  			return Indexed
   358  		}
   359  		if s.CellID() <= target.RangeMax() {
   360  			return Subdivided
   361  		}
   362  	}
   363  	if s.Prev() && s.CellID().RangeMax() >= target {
   364  		return Indexed
   365  	}
   366  	return Disjoint
   367  }
   368  
   369  // tracker keeps track of which shapes in a given set contain a particular point
   370  // (the focus). It provides an efficient way to move the focus from one point
   371  // to another and incrementally update the set of shapes which contain it. We use
   372  // this to compute which shapes contain the center of every CellID in the index,
   373  // by advancing the focus from one cell center to the next.
   374  //
   375  // Initially the focus is at the start of the CellID space-filling curve. We then
   376  // visit all the cells that are being added to the ShapeIndex in increasing order
   377  // of CellID. For each cell, we draw two edges: one from the entry vertex to the
   378  // center, and another from the center to the exit vertex (where entry and exit
   379  // refer to the points where the space-filling curve enters and exits the cell).
   380  // By counting edge crossings we can incrementally compute which shapes contain
   381  // the cell center. Note that the same set of shapes will always contain the exit
   382  // point of one cell and the entry point of the next cell in the index, because
   383  // either (a) these two points are actually the same, or (b) the intervening
   384  // cells in CellID order are all empty, and therefore there are no edge crossings
   385  // if we follow this path from one cell to the other.
   386  //
   387  // In C++, this is S2ShapeIndex::InteriorTracker.
   388  type tracker struct {
   389  	isActive   bool
   390  	a          Point
   391  	b          Point
   392  	nextCellID CellID
   393  	crosser    *EdgeCrosser
   394  	shapeIDs   []int32
   395  
   396  	// Shape ids saved by saveAndClearStateBefore. The state is never saved
   397  	// recursively so we don't need to worry about maintaining a stack.
   398  	savedIDs []int32
   399  }
   400  
   401  // newTracker returns a new tracker with the appropriate defaults.
   402  func newTracker() *tracker {
   403  	// As shapes are added, we compute which ones contain the start of the
   404  	// CellID space-filling curve by drawing an edge from OriginPoint to this
   405  	// point and counting how many shape edges cross this edge.
   406  	t := &tracker{
   407  		isActive:   false,
   408  		b:          trackerOrigin(),
   409  		nextCellID: CellIDFromFace(0).ChildBeginAtLevel(MaxLevel),
   410  	}
   411  	t.drawTo(Point{faceUVToXYZ(0, -1, -1).Normalize()}) // CellID curve start
   412  
   413  	return t
   414  }
   415  
   416  // trackerOrigin returns the initial focus point when the tracker is created
   417  // (corresponding to the start of the CellID space-filling curve).
   418  func trackerOrigin() Point {
   419  	// The start of the S2CellId space-filling curve.
   420  	return Point{faceUVToXYZ(0, -1, -1).Normalize()}
   421  }
   422  
   423  // focus returns the current focus point of the tracker.
   424  func (t *tracker) focus() Point { return t.b }
   425  
   426  // addShape adds a shape whose interior should be tracked. containsOrigin indicates
   427  // whether the current focus point is inside the shape. Alternatively, if
   428  // the focus point is in the process of being moved (via moveTo/drawTo), you
   429  // can also specify containsOrigin at the old focus point and call testEdge
   430  // for every edge of the shape that might cross the current drawTo line.
   431  // This updates the state to correspond to the new focus point.
   432  //
   433  // This requires shape.HasInterior
   434  func (t *tracker) addShape(shapeID int32, containsFocus bool) {
   435  	t.isActive = true
   436  	if containsFocus {
   437  		t.toggleShape(shapeID)
   438  	}
   439  }
   440  
   441  // moveTo moves the focus of the tracker to the given point. This method should
   442  // only be used when it is known that there are no edge crossings between the old
   443  // and new focus locations; otherwise use drawTo.
   444  func (t *tracker) moveTo(b Point) { t.b = b }
   445  
   446  // drawTo moves the focus of the tracker to the given point. After this method is
   447  // called, testEdge should be called with all edges that may cross the line
   448  // segment between the old and new focus locations.
   449  func (t *tracker) drawTo(b Point) {
   450  	t.a = t.b
   451  	t.b = b
   452  	// TODO: the edge crosser may need an in-place Init method if this gets expensive
   453  	t.crosser = NewEdgeCrosser(t.a, t.b)
   454  }
   455  
   456  // testEdge checks if the given edge crosses the current edge, and if so, then
   457  // toggle the state of the given shapeID.
   458  // This requires shape to have an interior.
   459  func (t *tracker) testEdge(shapeID int32, edge Edge) {
   460  	if t.crosser.EdgeOrVertexCrossing(edge.V0, edge.V1) {
   461  		t.toggleShape(shapeID)
   462  	}
   463  }
   464  
   465  // setNextCellID is used to indicate that the last argument to moveTo or drawTo
   466  // was the entry vertex of the given CellID, i.e. the tracker is positioned at the
   467  // start of this cell. By using this method together with atCellID, the caller
   468  // can avoid calling moveTo in cases where the exit vertex of the previous cell
   469  // is the same as the entry vertex of the current cell.
   470  func (t *tracker) setNextCellID(nextCellID CellID) {
   471  	t.nextCellID = nextCellID.RangeMin()
   472  }
   473  
   474  // atCellID reports if the focus is already at the entry vertex of the given
   475  // CellID (provided that the caller calls setNextCellID as each cell is processed).
   476  func (t *tracker) atCellID(cellid CellID) bool {
   477  	return cellid.RangeMin() == t.nextCellID
   478  }
   479  
   480  // toggleShape adds or removes the given shapeID from the set of IDs it is tracking.
   481  func (t *tracker) toggleShape(shapeID int32) {
   482  	// Most shapeIDs slices are small, so special case the common steps.
   483  
   484  	// If there is nothing here, add it.
   485  	if len(t.shapeIDs) == 0 {
   486  		t.shapeIDs = append(t.shapeIDs, shapeID)
   487  		return
   488  	}
   489  
   490  	// If it's the first element, drop it from the slice.
   491  	if t.shapeIDs[0] == shapeID {
   492  		t.shapeIDs = t.shapeIDs[1:]
   493  		return
   494  	}
   495  
   496  	for i, s := range t.shapeIDs {
   497  		if s < shapeID {
   498  			continue
   499  		}
   500  
   501  		// If it's in the set, cut it out.
   502  		if s == shapeID {
   503  			copy(t.shapeIDs[i:], t.shapeIDs[i+1:]) // overwrite the ith element
   504  			t.shapeIDs = t.shapeIDs[:len(t.shapeIDs)-1]
   505  			return
   506  		}
   507  
   508  		// We've got to a point in the slice where we should be inserted.
   509  		// (the given shapeID is now less than the current positions id.)
   510  		t.shapeIDs = append(t.shapeIDs[0:i],
   511  			append([]int32{shapeID}, t.shapeIDs[i:len(t.shapeIDs)]...)...)
   512  		return
   513  	}
   514  
   515  	// We got to the end and didn't find it, so add it to the list.
   516  	t.shapeIDs = append(t.shapeIDs, shapeID)
   517  }
   518  
   519  // saveAndClearStateBefore makes an internal copy of the state for shape ids below
   520  // the given limit, and then clear the state for those shapes. This is used during
   521  // incremental updates to track the state of added and removed shapes separately.
   522  func (t *tracker) saveAndClearStateBefore(limitShapeID int32) {
   523  	limit := t.lowerBound(limitShapeID)
   524  	t.savedIDs = append([]int32(nil), t.shapeIDs[:limit]...)
   525  	t.shapeIDs = t.shapeIDs[limit:]
   526  }
   527  
   528  // restoreStateBefore restores the state previously saved by saveAndClearStateBefore.
   529  // This only affects the state for shapeIDs below "limitShapeID".
   530  func (t *tracker) restoreStateBefore(limitShapeID int32) {
   531  	limit := t.lowerBound(limitShapeID)
   532  	t.shapeIDs = append(append([]int32(nil), t.savedIDs...), t.shapeIDs[limit:]...)
   533  	t.savedIDs = nil
   534  }
   535  
   536  // lowerBound returns the shapeID of the first entry x where x >= shapeID.
   537  func (t *tracker) lowerBound(shapeID int32) int32 {
   538  	panic("not implemented")
   539  }
   540  
   541  // removedShape represents a set of edges from the given shape that is queued for removal.
   542  type removedShape struct {
   543  	shapeID               int32
   544  	hasInterior           bool
   545  	containsTrackerOrigin bool
   546  	edges                 []Edge
   547  }
   548  
   549  // There are three basic states the index can be in.
   550  const (
   551  	stale    int32 = iota // There are pending updates.
   552  	updating              // Updates are currently being applied.
   553  	fresh                 // There are no pending updates.
   554  )
   555  
   556  // ShapeIndex indexes a set of Shapes, where a Shape is some collection of edges
   557  // that optionally defines an interior. It can be used to represent a set of
   558  // points, a set of polylines, or a set of polygons. For Shapes that have
   559  // interiors, the index makes it very fast to determine which Shape(s) contain
   560  // a given point or region.
   561  //
   562  // The index can be updated incrementally by adding or removing shapes. It is
   563  // designed to handle up to hundreds of millions of edges. All data structures
   564  // are designed to be small, so the index is compact; generally it is smaller
   565  // than the underlying data being indexed. The index is also fast to construct.
   566  //
   567  // Polygon, Loop, and Polyline implement Shape which allows these objects to
   568  // be indexed easily. You can find useful query methods in CrossingEdgeQuery
   569  // and ClosestEdgeQuery (Not yet implemented in Go).
   570  //
   571  // Example showing how to build an index of Polylines:
   572  //
   573  //	index := NewShapeIndex()
   574  //	for _, polyline := range polylines {
   575  //	    index.Add(polyline);
   576  //	}
   577  //	// Now you can use a CrossingEdgeQuery or ClosestEdgeQuery here.
   578  type ShapeIndex struct {
   579  	// shapes is a map of shape ID to shape.
   580  	shapes map[int32]Shape
   581  
   582  	// The maximum number of edges per cell.
   583  	// TODO(roberts): Update the comments when the usage of this is implemented.
   584  	maxEdgesPerCell int
   585  
   586  	// nextID tracks the next ID to hand out. IDs are not reused when shapes
   587  	// are removed from the index.
   588  	nextID int32
   589  
   590  	// cellMap is a map from CellID to the set of clipped shapes that intersect that
   591  	// cell. The cell IDs cover a set of non-overlapping regions on the sphere.
   592  	// In C++, this is a BTree, so the cells are ordered naturally by the data structure.
   593  	cellMap map[CellID]*ShapeIndexCell
   594  	// Track the ordered list of cell IDs.
   595  	cells []CellID
   596  
   597  	// The current status of the index; accessed atomically.
   598  	status int32
   599  
   600  	// Additions and removals are queued and processed on the first subsequent
   601  	// query. There are several reasons to do this:
   602  	//
   603  	//  - It is significantly more efficient to process updates in batches if
   604  	//    the amount of entities added grows.
   605  	//  - Often the index will never be queried, in which case we can save both
   606  	//    the time and memory required to build it. Examples:
   607  	//     + Loops that are created simply to pass to an Polygon. (We don't
   608  	//       need the Loop index, because Polygon builds its own index.)
   609  	//     + Applications that load a database of geometry and then query only
   610  	//       a small fraction of it.
   611  	//
   612  	// The main drawback is that we need to go to some extra work to ensure that
   613  	// some methods are still thread-safe. Note that the goal is *not* to
   614  	// make this thread-safe in general, but simply to hide the fact that
   615  	// we defer some of the indexing work until query time.
   616  	//
   617  	// This mutex protects all of following fields in the index.
   618  	mu sync.RWMutex
   619  
   620  	// pendingAdditionsPos is the index of the first entry that has not been processed
   621  	// via applyUpdatesInternal.
   622  	pendingAdditionsPos int32
   623  
   624  	// The set of shapes that have been queued for removal but not processed yet by
   625  	// applyUpdatesInternal.
   626  	pendingRemovals []*removedShape
   627  }
   628  
   629  // NewShapeIndex creates a new ShapeIndex.
   630  func NewShapeIndex() *ShapeIndex {
   631  	return &ShapeIndex{
   632  		maxEdgesPerCell: 10,
   633  		shapes:          make(map[int32]Shape),
   634  		cellMap:         make(map[CellID]*ShapeIndexCell),
   635  		cells:           nil,
   636  		status:          fresh,
   637  	}
   638  }
   639  
   640  // Iterator returns an iterator for this index.
   641  func (s *ShapeIndex) Iterator() *ShapeIndexIterator {
   642  	s.maybeApplyUpdates()
   643  	return NewShapeIndexIterator(s, IteratorBegin)
   644  }
   645  
   646  // Begin positions the iterator at the first cell in the index.
   647  func (s *ShapeIndex) Begin() *ShapeIndexIterator {
   648  	s.maybeApplyUpdates()
   649  	return NewShapeIndexIterator(s, IteratorBegin)
   650  }
   651  
   652  // End positions the iterator at the last cell in the index.
   653  func (s *ShapeIndex) End() *ShapeIndexIterator {
   654  	// TODO(roberts): It's possible that updates could happen to the index between
   655  	// the time this is called and the time the iterators position is used and this
   656  	// will be invalid or not the end. For now, things will be undefined if this
   657  	// happens. See about referencing the IsFresh to guard for this in the future.
   658  	s.maybeApplyUpdates()
   659  	return NewShapeIndexIterator(s, IteratorEnd)
   660  }
   661  
   662  // Region returns a new ShapeIndexRegion for this ShapeIndex.
   663  func (s *ShapeIndex) Region() *ShapeIndexRegion {
   664  	return &ShapeIndexRegion{
   665  		index:         s,
   666  		containsQuery: NewContainsPointQuery(s, VertexModelSemiOpen),
   667  		iter:          s.Iterator(),
   668  	}
   669  }
   670  
   671  // Len reports the number of Shapes in this index.
   672  func (s *ShapeIndex) Len() int {
   673  	return len(s.shapes)
   674  }
   675  
   676  // Reset resets the index to its original state.
   677  func (s *ShapeIndex) Reset() {
   678  	s.shapes = make(map[int32]Shape)
   679  	s.nextID = 0
   680  	s.cellMap = make(map[CellID]*ShapeIndexCell)
   681  	s.cells = nil
   682  	atomic.StoreInt32(&s.status, fresh)
   683  }
   684  
   685  // NumEdges returns the number of edges in this index.
   686  func (s *ShapeIndex) NumEdges() int {
   687  	numEdges := 0
   688  	for _, shape := range s.shapes {
   689  		numEdges += shape.NumEdges()
   690  	}
   691  	return numEdges
   692  }
   693  
   694  // NumEdgesUpTo returns the number of edges in the given index, up to the given
   695  // limit. If the limit is encountered, the current running total is returned,
   696  // which may be more than the limit.
   697  func (s *ShapeIndex) NumEdgesUpTo(limit int) int {
   698  	var numEdges int
   699  	// We choose to iterate over the shapes in order to match the counting
   700  	// up behavior in C++ and for test compatibility instead of using a
   701  	// more idiomatic range over the shape map.
   702  	for i := int32(0); i <= s.nextID; i++ {
   703  		s := s.Shape(i)
   704  		if s == nil {
   705  			continue
   706  		}
   707  		numEdges += s.NumEdges()
   708  		if numEdges >= limit {
   709  			break
   710  		}
   711  	}
   712  
   713  	return numEdges
   714  }
   715  
   716  // Shape returns the shape with the given ID, or nil if the shape has been removed from the index.
   717  func (s *ShapeIndex) Shape(id int32) Shape { return s.shapes[id] }
   718  
   719  // idForShape returns the id of the given shape in this index, or -1 if it is
   720  // not in the index.
   721  //
   722  // TODO(roberts): Need to figure out an appropriate way to expose this on a Shape.
   723  // C++ allows a given S2 type (Loop, Polygon, etc) to be part of multiple indexes.
   724  // By having each type extend S2Shape which has an id element, they all inherit their
   725  // own id field rather than having to track it themselves.
   726  func (s *ShapeIndex) idForShape(shape Shape) int32 {
   727  	for k, v := range s.shapes {
   728  		if v == shape {
   729  			return k
   730  		}
   731  	}
   732  	return -1
   733  }
   734  
   735  // Add adds the given shape to the index and returns the assigned ID..
   736  func (s *ShapeIndex) Add(shape Shape) int32 {
   737  	s.shapes[s.nextID] = shape
   738  	s.nextID++
   739  	atomic.StoreInt32(&s.status, stale)
   740  	return s.nextID - 1
   741  }
   742  
   743  // Remove removes the given shape from the index.
   744  func (s *ShapeIndex) Remove(shape Shape) {
   745  	// The index updates itself lazily because it is much more efficient to
   746  	// process additions and removals in batches.
   747  	id := s.idForShape(shape)
   748  
   749  	// If the shape wasn't found, it's already been removed or was not in the index.
   750  	if s.shapes[id] == nil {
   751  		return
   752  	}
   753  
   754  	// Remove the shape from the shapes map.
   755  	delete(s.shapes, id)
   756  
   757  	// We are removing a shape that has not yet been added to the index,
   758  	// so there is nothing else to do.
   759  	if id >= s.pendingAdditionsPos {
   760  		return
   761  	}
   762  
   763  	numEdges := shape.NumEdges()
   764  	removed := &removedShape{
   765  		shapeID:               id,
   766  		hasInterior:           shape.Dimension() == 2,
   767  		containsTrackerOrigin: shape.ReferencePoint().Contained,
   768  		edges:                 make([]Edge, numEdges),
   769  	}
   770  
   771  	for e := 0; e < numEdges; e++ {
   772  		removed.edges[e] = shape.Edge(e)
   773  	}
   774  
   775  	s.pendingRemovals = append(s.pendingRemovals, removed)
   776  	atomic.StoreInt32(&s.status, stale)
   777  }
   778  
   779  // Build triggers the update of the index. Calls to Add and Release are normally
   780  // queued and processed on the first subsequent query. This has many advantages,
   781  // the most important of which is that sometimes there *is* no subsequent
   782  // query, which lets us avoid building the index completely.
   783  //
   784  // This method forces any pending updates to be applied immediately.
   785  func (s *ShapeIndex) Build() {
   786  	s.maybeApplyUpdates()
   787  }
   788  
   789  // IsFresh reports if there are no pending updates that need to be applied.
   790  // This can be useful to avoid building the index unnecessarily, or for
   791  // choosing between two different algorithms depending on whether the index
   792  // is available.
   793  //
   794  // The returned index status may be slightly out of date if the index was
   795  // built in a different thread. This is fine for the intended use (as an
   796  // efficiency hint), but it should not be used by internal methods.
   797  func (s *ShapeIndex) IsFresh() bool {
   798  	return atomic.LoadInt32(&s.status) == fresh
   799  }
   800  
   801  // isFirstUpdate reports if this is the first update to the index.
   802  func (s *ShapeIndex) isFirstUpdate() bool {
   803  	// Note that it is not sufficient to check whether cellMap is empty, since
   804  	// entries are added to it during the update process.
   805  	return s.pendingAdditionsPos == 0
   806  }
   807  
   808  // isShapeBeingRemoved reports if the shape with the given ID is currently slated for removal.
   809  func (s *ShapeIndex) isShapeBeingRemoved(shapeID int32) bool {
   810  	// All shape ids being removed fall below the index position of shapes being added.
   811  	return shapeID < s.pendingAdditionsPos
   812  }
   813  
   814  // maybeApplyUpdates checks if the index pieces have changed, and if so, applies pending updates.
   815  func (s *ShapeIndex) maybeApplyUpdates() {
   816  	// TODO(roberts): To avoid acquiring and releasing the mutex on every
   817  	// query, we should use atomic operations when testing whether the status
   818  	// is fresh and when updating the status to be fresh. This guarantees
   819  	// that any thread that sees a status of fresh will also see the
   820  	// corresponding index updates.
   821  	if atomic.LoadInt32(&s.status) != fresh {
   822  		s.mu.Lock()
   823  		s.applyUpdatesInternal()
   824  		atomic.StoreInt32(&s.status, fresh)
   825  		s.mu.Unlock()
   826  	}
   827  }
   828  
   829  // applyUpdatesInternal does the actual work of updating the index by applying all
   830  // pending additions and removals. It does *not* update the indexes status.
   831  func (s *ShapeIndex) applyUpdatesInternal() {
   832  	// TODO(roberts): Building the index can use up to 20x as much memory per
   833  	// edge as the final index memory size. If this causes issues, add in
   834  	// batched updating to limit the amount of items per batch to a
   835  	// configurable memory footprint overhead.
   836  	t := newTracker()
   837  
   838  	// allEdges maps a Face to a collection of faceEdges.
   839  	allEdges := make([][]faceEdge, 6)
   840  
   841  	for _, p := range s.pendingRemovals {
   842  		s.removeShapeInternal(p, allEdges, t)
   843  	}
   844  
   845  	for id := s.pendingAdditionsPos; id < int32(len(s.shapes)); id++ {
   846  		s.addShapeInternal(id, allEdges, t)
   847  	}
   848  
   849  	for face := 0; face < 6; face++ {
   850  		s.updateFaceEdges(face, allEdges[face], t)
   851  	}
   852  
   853  	s.pendingRemovals = s.pendingRemovals[:0]
   854  	s.pendingAdditionsPos = int32(len(s.shapes))
   855  	// It is the caller's responsibility to update the index status.
   856  }
   857  
   858  // addShapeInternal clips all edges of the given shape to the six cube faces,
   859  // adds the clipped edges to the set of allEdges, and starts tracking its
   860  // interior if necessary.
   861  func (s *ShapeIndex) addShapeInternal(shapeID int32, allEdges [][]faceEdge, t *tracker) {
   862  	shape, ok := s.shapes[shapeID]
   863  	if !ok {
   864  		// This shape has already been removed.
   865  		return
   866  	}
   867  
   868  	faceEdge := faceEdge{
   869  		shapeID:     shapeID,
   870  		hasInterior: shape.Dimension() == 2,
   871  	}
   872  
   873  	if faceEdge.hasInterior {
   874  		t.addShape(shapeID, containsBruteForce(shape, t.focus()))
   875  	}
   876  
   877  	numEdges := shape.NumEdges()
   878  	for e := 0; e < numEdges; e++ {
   879  		edge := shape.Edge(e)
   880  
   881  		faceEdge.edgeID = e
   882  		faceEdge.edge = edge
   883  		faceEdge.MaxLevel = maxLevelForEdge(edge)
   884  		s.addFaceEdge(faceEdge, allEdges)
   885  	}
   886  }
   887  
   888  // addFaceEdge adds the given faceEdge into the collection of all edges.
   889  func (s *ShapeIndex) addFaceEdge(fe faceEdge, allEdges [][]faceEdge) {
   890  	aFace := face(fe.edge.V0.Vector)
   891  	// See if both endpoints are on the same face, and are far enough from
   892  	// the edge of the face that they don't intersect any (padded) adjacent face.
   893  	if aFace == face(fe.edge.V1.Vector) {
   894  		x, y := validFaceXYZToUV(aFace, fe.edge.V0.Vector)
   895  		fe.a = r2.Point{x, y}
   896  		x, y = validFaceXYZToUV(aFace, fe.edge.V1.Vector)
   897  		fe.b = r2.Point{x, y}
   898  
   899  		maxUV := 1 - cellPadding
   900  		if math.Abs(fe.a.X) <= maxUV && math.Abs(fe.a.Y) <= maxUV &&
   901  			math.Abs(fe.b.X) <= maxUV && math.Abs(fe.b.Y) <= maxUV {
   902  			allEdges[aFace] = append(allEdges[aFace], fe)
   903  			return
   904  		}
   905  	}
   906  
   907  	// Otherwise, we simply clip the edge to all six faces.
   908  	for face := 0; face < 6; face++ {
   909  		if aClip, bClip, intersects := ClipToPaddedFace(fe.edge.V0, fe.edge.V1, face, cellPadding); intersects {
   910  			fe.a = aClip
   911  			fe.b = bClip
   912  			allEdges[face] = append(allEdges[face], fe)
   913  		}
   914  	}
   915  }
   916  
   917  // updateFaceEdges adds or removes the various edges from the index.
   918  // An edge is added if shapes[id] is not nil, and removed otherwise.
   919  func (s *ShapeIndex) updateFaceEdges(face int, faceEdges []faceEdge, t *tracker) {
   920  	numEdges := len(faceEdges)
   921  	if numEdges == 0 && len(t.shapeIDs) == 0 {
   922  		return
   923  	}
   924  
   925  	// Create the initial clippedEdge for each faceEdge. Additional clipped
   926  	// edges are created when edges are split between child cells. We create
   927  	// two arrays, one containing the edge data and another containing pointers
   928  	// to those edges, so that during the recursion we only need to copy
   929  	// pointers in order to propagate an edge to the correct child.
   930  	clippedEdges := make([]*clippedEdge, numEdges)
   931  	bound := r2.EmptyRect()
   932  	for e := 0; e < numEdges; e++ {
   933  		clipped := &clippedEdge{
   934  			faceEdge: &faceEdges[e],
   935  		}
   936  		clipped.bound = r2.RectFromPoints(faceEdges[e].a, faceEdges[e].b)
   937  		clippedEdges[e] = clipped
   938  		bound = bound.AddRect(clipped.bound)
   939  	}
   940  
   941  	// Construct the initial face cell containing all the edges, and then update
   942  	// all the edges in the index recursively.
   943  	faceID := CellIDFromFace(face)
   944  	pcell := PaddedCellFromCellID(faceID, cellPadding)
   945  
   946  	disjointFromIndex := s.isFirstUpdate()
   947  	if numEdges > 0 {
   948  		shrunkID := s.shrinkToFit(pcell, bound)
   949  		if shrunkID != pcell.id {
   950  			// All the edges are contained by some descendant of the face cell. We
   951  			// can save a lot of work by starting directly with that cell, but if we
   952  			// are in the interior of at least one shape then we need to create
   953  			// index entries for the cells we are skipping over.
   954  			s.skipCellRange(faceID.RangeMin(), shrunkID.RangeMin(), t, disjointFromIndex)
   955  			pcell = PaddedCellFromCellID(shrunkID, cellPadding)
   956  			s.updateEdges(pcell, clippedEdges, t, disjointFromIndex)
   957  			s.skipCellRange(shrunkID.RangeMax().Next(), faceID.RangeMax().Next(), t, disjointFromIndex)
   958  			return
   959  		}
   960  	}
   961  
   962  	// Otherwise (no edges, or no shrinking is possible), subdivide normally.
   963  	s.updateEdges(pcell, clippedEdges, t, disjointFromIndex)
   964  }
   965  
   966  // shrinkToFit shrinks the PaddedCell to fit within the given bounds.
   967  func (s *ShapeIndex) shrinkToFit(pcell *PaddedCell, bound r2.Rect) CellID {
   968  	shrunkID := pcell.ShrinkToFit(bound)
   969  
   970  	if !s.isFirstUpdate() && shrunkID != pcell.CellID() {
   971  		// Don't shrink any smaller than the existing index cells, since we need
   972  		// to combine the new edges with those cells.
   973  		iter := s.Iterator()
   974  		if iter.LocateCellID(shrunkID) == Indexed {
   975  			shrunkID = iter.CellID()
   976  		}
   977  	}
   978  	return shrunkID
   979  }
   980  
   981  // skipCellRange skips over the cells in the given range, creating index cells if we are
   982  // currently in the interior of at least one shape.
   983  func (s *ShapeIndex) skipCellRange(begin, end CellID, t *tracker, disjointFromIndex bool) {
   984  	// If we aren't in the interior of a shape, then skipping over cells is easy.
   985  	if len(t.shapeIDs) == 0 {
   986  		return
   987  	}
   988  
   989  	// Otherwise generate the list of cell ids that we need to visit, and create
   990  	// an index entry for each one.
   991  	skipped := CellUnionFromRange(begin, end)
   992  	for _, cell := range skipped {
   993  		var clippedEdges []*clippedEdge
   994  		s.updateEdges(PaddedCellFromCellID(cell, cellPadding), clippedEdges, t, disjointFromIndex)
   995  	}
   996  }
   997  
   998  // updateEdges adds or removes the given edges whose bounding boxes intersect a
   999  // given cell. disjointFromIndex is an optimization hint indicating that cellMap
  1000  // does not contain any entries that overlap the given cell.
  1001  func (s *ShapeIndex) updateEdges(pcell *PaddedCell, edges []*clippedEdge, t *tracker, disjointFromIndex bool) {
  1002  	// This function is recursive with a maximum recursion depth of 30 (MaxLevel).
  1003  
  1004  	// Incremental updates are handled as follows. All edges being added or
  1005  	// removed are combined together in edges, and all shapes with interiors
  1006  	// are tracked using tracker. We subdivide recursively as usual until we
  1007  	// encounter an existing index cell. At this point we absorb the index
  1008  	// cell as follows:
  1009  	//
  1010  	//   - Edges and shapes that are being removed are deleted from edges and
  1011  	//     tracker.
  1012  	//   - All remaining edges and shapes from the index cell are added to
  1013  	//     edges and tracker.
  1014  	//   - Continue subdividing recursively, creating new index cells as needed.
  1015  	//   - When the recursion gets back to the cell that was absorbed, we
  1016  	//     restore edges and tracker to their previous state.
  1017  	//
  1018  	// Note that the only reason that we include removed shapes in the recursive
  1019  	// subdivision process is so that we can find all of the index cells that
  1020  	// contain those shapes efficiently, without maintaining an explicit list of
  1021  	// index cells for each shape (which would be expensive in terms of memory).
  1022  	indexCellAbsorbed := false
  1023  	if !disjointFromIndex {
  1024  		// There may be existing index cells contained inside pcell. If we
  1025  		// encounter such a cell, we need to combine the edges being updated with
  1026  		// the existing cell contents by absorbing the cell.
  1027  		iter := s.Iterator()
  1028  		r := iter.LocateCellID(pcell.id)
  1029  		if r == Disjoint {
  1030  			disjointFromIndex = true
  1031  		} else if r == Indexed {
  1032  			// Absorb the index cell by transferring its contents to edges and
  1033  			// deleting it. We also start tracking the interior of any new shapes.
  1034  			s.absorbIndexCell(pcell, iter, edges, t)
  1035  			indexCellAbsorbed = true
  1036  			disjointFromIndex = true
  1037  		} else {
  1038  			// DCHECK_EQ(SUBDIVIDED, r)
  1039  		}
  1040  	}
  1041  
  1042  	// If there are existing index cells below us, then we need to keep
  1043  	// subdividing so that we can merge with those cells. Otherwise,
  1044  	// makeIndexCell checks if the number of edges is small enough, and creates
  1045  	// an index cell if possible (returning true when it does so).
  1046  	if !disjointFromIndex || !s.makeIndexCell(pcell, edges, t) {
  1047  		// TODO(roberts): If it turns out to have memory problems when there
  1048  		// are 10M+ edges in the index, look into pre-allocating space so we
  1049  		// are not always appending.
  1050  		childEdges := [2][2][]*clippedEdge{} // [i][j]
  1051  
  1052  		// Compute the middle of the padded cell, defined as the rectangle in
  1053  		// (u,v)-space that belongs to all four (padded) children. By comparing
  1054  		// against the four boundaries of middle we can determine which children
  1055  		// each edge needs to be propagated to.
  1056  		middle := pcell.Middle()
  1057  
  1058  		// Build up a vector edges to be passed to each child cell. The (i,j)
  1059  		// directions are left (i=0), right (i=1), lower (j=0), and upper (j=1).
  1060  		// Note that the vast majority of edges are propagated to a single child.
  1061  		for _, edge := range edges {
  1062  			if edge.bound.X.Hi <= middle.X.Lo {
  1063  				// Edge is entirely contained in the two left children.
  1064  				a, b := s.clipVAxis(edge, middle.Y)
  1065  				if a != nil {
  1066  					childEdges[0][0] = append(childEdges[0][0], a)
  1067  				}
  1068  				if b != nil {
  1069  					childEdges[0][1] = append(childEdges[0][1], b)
  1070  				}
  1071  			} else if edge.bound.X.Lo >= middle.X.Hi {
  1072  				// Edge is entirely contained in the two right children.
  1073  				a, b := s.clipVAxis(edge, middle.Y)
  1074  				if a != nil {
  1075  					childEdges[1][0] = append(childEdges[1][0], a)
  1076  				}
  1077  				if b != nil {
  1078  					childEdges[1][1] = append(childEdges[1][1], b)
  1079  				}
  1080  			} else if edge.bound.Y.Hi <= middle.Y.Lo {
  1081  				// Edge is entirely contained in the two lower children.
  1082  				if a := s.clipUBound(edge, 1, middle.X.Hi); a != nil {
  1083  					childEdges[0][0] = append(childEdges[0][0], a)
  1084  				}
  1085  				if b := s.clipUBound(edge, 0, middle.X.Lo); b != nil {
  1086  					childEdges[1][0] = append(childEdges[1][0], b)
  1087  				}
  1088  			} else if edge.bound.Y.Lo >= middle.Y.Hi {
  1089  				// Edge is entirely contained in the two upper children.
  1090  				if a := s.clipUBound(edge, 1, middle.X.Hi); a != nil {
  1091  					childEdges[0][1] = append(childEdges[0][1], a)
  1092  				}
  1093  				if b := s.clipUBound(edge, 0, middle.X.Lo); b != nil {
  1094  					childEdges[1][1] = append(childEdges[1][1], b)
  1095  				}
  1096  			} else {
  1097  				// The edge bound spans all four children. The edge
  1098  				// itself intersects either three or four padded children.
  1099  				left := s.clipUBound(edge, 1, middle.X.Hi)
  1100  				a, b := s.clipVAxis(left, middle.Y)
  1101  				if a != nil {
  1102  					childEdges[0][0] = append(childEdges[0][0], a)
  1103  				}
  1104  				if b != nil {
  1105  					childEdges[0][1] = append(childEdges[0][1], b)
  1106  				}
  1107  				right := s.clipUBound(edge, 0, middle.X.Lo)
  1108  				a, b = s.clipVAxis(right, middle.Y)
  1109  				if a != nil {
  1110  					childEdges[1][0] = append(childEdges[1][0], a)
  1111  				}
  1112  				if b != nil {
  1113  					childEdges[1][1] = append(childEdges[1][1], b)
  1114  				}
  1115  			}
  1116  		}
  1117  
  1118  		// Now recursively update the edges in each child. We call the children in
  1119  		// increasing order of CellID so that when the index is first constructed,
  1120  		// all insertions into cellMap are at the end (which is much faster).
  1121  		for pos := 0; pos < 4; pos++ {
  1122  			i, j := pcell.ChildIJ(pos)
  1123  			if len(childEdges[i][j]) > 0 || len(t.shapeIDs) > 0 {
  1124  				s.updateEdges(PaddedCellFromParentIJ(pcell, i, j), childEdges[i][j],
  1125  					t, disjointFromIndex)
  1126  			}
  1127  		}
  1128  	}
  1129  
  1130  	if indexCellAbsorbed {
  1131  		// Restore the state for any edges being removed that we are tracking.
  1132  		t.restoreStateBefore(s.pendingAdditionsPos)
  1133  	}
  1134  }
  1135  
  1136  // makeIndexCell builds an indexCell from the given padded cell and set of edges and adds
  1137  // it to the index. If the cell or edges are empty, no cell is added.
  1138  func (s *ShapeIndex) makeIndexCell(p *PaddedCell, edges []*clippedEdge, t *tracker) bool {
  1139  	// If the cell is empty, no index cell is needed. (In most cases this
  1140  	// situation is detected before we get to this point, but this can happen
  1141  	// when all shapes in a cell are removed.)
  1142  	if len(edges) == 0 && len(t.shapeIDs) == 0 {
  1143  		return true
  1144  	}
  1145  
  1146  	// Count the number of edges that have not reached their maximum level yet.
  1147  	// Return false if there are too many such edges.
  1148  	count := 0
  1149  	for _, ce := range edges {
  1150  		if p.Level() < ce.faceEdge.MaxLevel {
  1151  			count++
  1152  		}
  1153  
  1154  		if count > s.maxEdgesPerCell {
  1155  			return false
  1156  		}
  1157  	}
  1158  
  1159  	// Possible optimization: Continue subdividing as long as exactly one child
  1160  	// of the padded cell intersects the given edges. This can be done by finding
  1161  	// the bounding box of all the edges and calling ShrinkToFit:
  1162  	//
  1163  	// cellID = p.ShrinkToFit(RectBound(edges));
  1164  	//
  1165  	// Currently this is not beneficial; it slows down construction by 4-25%
  1166  	// (mainly computing the union of the bounding rectangles) and also slows
  1167  	// down queries (since more recursive clipping is required to get down to
  1168  	// the level of a spatial index cell). But it may be worth trying again
  1169  	// once containsCenter is computed and all algorithms are modified to
  1170  	// take advantage of it.
  1171  
  1172  	// We update the InteriorTracker as follows. For every Cell in the index
  1173  	// we construct two edges: one edge from entry vertex of the cell to its
  1174  	// center, and one from the cell center to its exit vertex. Here entry
  1175  	// and exit refer the CellID ordering, i.e. the order in which points
  1176  	// are encountered along the 2 space-filling curve. The exit vertex then
  1177  	// becomes the entry vertex for the next cell in the index, unless there are
  1178  	// one or more empty intervening cells, in which case the InteriorTracker
  1179  	// state is unchanged because the intervening cells have no edges.
  1180  
  1181  	// Shift the InteriorTracker focus point to the center of the current cell.
  1182  	if t.isActive && len(edges) != 0 {
  1183  		if !t.atCellID(p.id) {
  1184  			t.moveTo(p.EntryVertex())
  1185  		}
  1186  		t.drawTo(p.Center())
  1187  		s.testAllEdges(edges, t)
  1188  	}
  1189  
  1190  	// Allocate and fill a new index cell. To get the total number of shapes we
  1191  	// need to merge the shapes associated with the intersecting edges together
  1192  	// with the shapes that happen to contain the cell center.
  1193  	cshapeIDs := t.shapeIDs
  1194  	numShapes := s.countShapes(edges, cshapeIDs)
  1195  	cell := NewShapeIndexCell(numShapes)
  1196  
  1197  	// To fill the index cell we merge the two sources of shapes: edge shapes
  1198  	// (those that have at least one edge that intersects this cell), and
  1199  	// containing shapes (those that contain the cell center). We keep track
  1200  	// of the index of the next intersecting edge and the next containing shape
  1201  	// as we go along. Both sets of shape ids are already sorted.
  1202  	eNext := 0
  1203  	cNextIdx := 0
  1204  	for i := 0; i < numShapes; i++ {
  1205  		var clipped *clippedShape
  1206  		// advance to next value base + i
  1207  		eshapeID := int32(s.Len())
  1208  		cshapeID := eshapeID // Sentinels
  1209  
  1210  		if eNext != len(edges) {
  1211  			eshapeID = edges[eNext].faceEdge.shapeID
  1212  		}
  1213  		if cNextIdx < len(cshapeIDs) {
  1214  			cshapeID = cshapeIDs[cNextIdx]
  1215  		}
  1216  		eBegin := eNext
  1217  		if cshapeID < eshapeID {
  1218  			// The entire cell is in the shape interior.
  1219  			clipped = newClippedShape(cshapeID, 0)
  1220  			clipped.containsCenter = true
  1221  			cNextIdx++
  1222  		} else {
  1223  			// Count the number of edges for this shape and allocate space for them.
  1224  			for eNext < len(edges) && edges[eNext].faceEdge.shapeID == eshapeID {
  1225  				eNext++
  1226  			}
  1227  			clipped = newClippedShape(eshapeID, eNext-eBegin)
  1228  			for e := eBegin; e < eNext; e++ {
  1229  				clipped.edges[e-eBegin] = edges[e].faceEdge.edgeID
  1230  			}
  1231  			if cshapeID == eshapeID {
  1232  				clipped.containsCenter = true
  1233  				cNextIdx++
  1234  			}
  1235  		}
  1236  		cell.shapes[i] = clipped
  1237  	}
  1238  
  1239  	// Add this cell to the map.
  1240  	s.cellMap[p.id] = cell
  1241  	s.cells = append(s.cells, p.id)
  1242  
  1243  	// Shift the tracker focus point to the exit vertex of this cell.
  1244  	if t.isActive && len(edges) != 0 {
  1245  		t.drawTo(p.ExitVertex())
  1246  		s.testAllEdges(edges, t)
  1247  		t.setNextCellID(p.id.Next())
  1248  	}
  1249  	return true
  1250  }
  1251  
  1252  // updateBound updates the specified endpoint of the given clipped edge and returns the
  1253  // resulting clipped edge.
  1254  func (s *ShapeIndex) updateBound(edge *clippedEdge, uEnd int, u float64, vEnd int, v float64) *clippedEdge {
  1255  	c := &clippedEdge{faceEdge: edge.faceEdge}
  1256  	if uEnd == 0 {
  1257  		c.bound.X.Lo = u
  1258  		c.bound.X.Hi = edge.bound.X.Hi
  1259  	} else {
  1260  		c.bound.X.Lo = edge.bound.X.Lo
  1261  		c.bound.X.Hi = u
  1262  	}
  1263  
  1264  	if vEnd == 0 {
  1265  		c.bound.Y.Lo = v
  1266  		c.bound.Y.Hi = edge.bound.Y.Hi
  1267  	} else {
  1268  		c.bound.Y.Lo = edge.bound.Y.Lo
  1269  		c.bound.Y.Hi = v
  1270  	}
  1271  
  1272  	return c
  1273  }
  1274  
  1275  // clipUBound clips the given endpoint (lo=0, hi=1) of the u-axis so that
  1276  // it does not extend past the given value of the given edge.
  1277  func (s *ShapeIndex) clipUBound(edge *clippedEdge, uEnd int, u float64) *clippedEdge {
  1278  	// First check whether the edge actually requires any clipping. (Sometimes
  1279  	// this method is called when clipping is not necessary, e.g. when one edge
  1280  	// endpoint is in the overlap area between two padded child cells.)
  1281  	if uEnd == 0 {
  1282  		if edge.bound.X.Lo >= u {
  1283  			return edge
  1284  		}
  1285  	} else {
  1286  		if edge.bound.X.Hi <= u {
  1287  			return edge
  1288  		}
  1289  	}
  1290  	// We interpolate the new v-value from the endpoints of the original edge.
  1291  	// This has two advantages: (1) we don't need to store the clipped endpoints
  1292  	// at all, just their bounding box; and (2) it avoids the accumulation of
  1293  	// roundoff errors due to repeated interpolations. The result needs to be
  1294  	// clamped to ensure that it is in the appropriate range.
  1295  	e := edge.faceEdge
  1296  	v := edge.bound.Y.ClampPoint(interpolateFloat64(u, e.a.X, e.b.X, e.a.Y, e.b.Y))
  1297  
  1298  	// Determine which endpoint of the v-axis bound to update. If the edge
  1299  	// slope is positive we update the same endpoint, otherwise we update the
  1300  	// opposite endpoint.
  1301  	var vEnd int
  1302  	positiveSlope := (e.a.X > e.b.X) == (e.a.Y > e.b.Y)
  1303  	if (uEnd == 1) == positiveSlope {
  1304  		vEnd = 1
  1305  	}
  1306  	return s.updateBound(edge, uEnd, u, vEnd, v)
  1307  }
  1308  
  1309  // clipVBound clips the given endpoint (lo=0, hi=1) of the v-axis so that
  1310  // it does not extend past the given value of the given edge.
  1311  func (s *ShapeIndex) clipVBound(edge *clippedEdge, vEnd int, v float64) *clippedEdge {
  1312  	if vEnd == 0 {
  1313  		if edge.bound.Y.Lo >= v {
  1314  			return edge
  1315  		}
  1316  	} else {
  1317  		if edge.bound.Y.Hi <= v {
  1318  			return edge
  1319  		}
  1320  	}
  1321  
  1322  	// We interpolate the new v-value from the endpoints of the original edge.
  1323  	// This has two advantages: (1) we don't need to store the clipped endpoints
  1324  	// at all, just their bounding box; and (2) it avoids the accumulation of
  1325  	// roundoff errors due to repeated interpolations. The result needs to be
  1326  	// clamped to ensure that it is in the appropriate range.
  1327  	e := edge.faceEdge
  1328  	u := edge.bound.X.ClampPoint(interpolateFloat64(v, e.a.Y, e.b.Y, e.a.X, e.b.X))
  1329  
  1330  	// Determine which endpoint of the v-axis bound to update. If the edge
  1331  	// slope is positive we update the same endpoint, otherwise we update the
  1332  	// opposite endpoint.
  1333  	var uEnd int
  1334  	positiveSlope := (e.a.X > e.b.X) == (e.a.Y > e.b.Y)
  1335  	if (vEnd == 1) == positiveSlope {
  1336  		uEnd = 1
  1337  	}
  1338  	return s.updateBound(edge, uEnd, u, vEnd, v)
  1339  }
  1340  
  1341  // cliupVAxis returns the given edge clipped to within the boundaries of the middle
  1342  // interval along the v-axis, and adds the result to its children.
  1343  func (s *ShapeIndex) clipVAxis(edge *clippedEdge, middle r1.Interval) (a, b *clippedEdge) {
  1344  	if edge.bound.Y.Hi <= middle.Lo {
  1345  		// Edge is entirely contained in the lower child.
  1346  		return edge, nil
  1347  	} else if edge.bound.Y.Lo >= middle.Hi {
  1348  		// Edge is entirely contained in the upper child.
  1349  		return nil, edge
  1350  	}
  1351  	// The edge bound spans both children.
  1352  	return s.clipVBound(edge, 1, middle.Hi), s.clipVBound(edge, 0, middle.Lo)
  1353  }
  1354  
  1355  // absorbIndexCell absorbs an index cell by transferring its contents to edges
  1356  // and/or "tracker", and then delete this cell from the index. If edges includes
  1357  // any edges that are being removed, this method also updates their
  1358  // InteriorTracker state to correspond to the exit vertex of this cell.
  1359  func (s *ShapeIndex) absorbIndexCell(p *PaddedCell, iter *ShapeIndexIterator, edges []*clippedEdge, t *tracker) {
  1360  	// When we absorb a cell, we erase all the edges that are being removed.
  1361  	// However when we are finished with this cell, we want to restore the state
  1362  	// of those edges (since that is how we find all the index cells that need
  1363  	// to be updated).  The edges themselves are restored automatically when
  1364  	// UpdateEdges returns from its recursive call, but the InteriorTracker
  1365  	// state needs to be restored explicitly.
  1366  	//
  1367  	// Here we first update the InteriorTracker state for removed edges to
  1368  	// correspond to the exit vertex of this cell, and then save the
  1369  	// InteriorTracker state.  This state will be restored by UpdateEdges when
  1370  	// it is finished processing the contents of this cell.
  1371  	if t.isActive && len(edges) != 0 && s.isShapeBeingRemoved(edges[0].faceEdge.shapeID) {
  1372  		// We probably need to update the tracker. ("Probably" because
  1373  		// it's possible that all shapes being removed do not have interiors.)
  1374  		if !t.atCellID(p.id) {
  1375  			t.moveTo(p.EntryVertex())
  1376  		}
  1377  		t.drawTo(p.ExitVertex())
  1378  		t.setNextCellID(p.id.Next())
  1379  		for _, edge := range edges {
  1380  			fe := edge.faceEdge
  1381  			if !s.isShapeBeingRemoved(fe.shapeID) {
  1382  				break // All shapes being removed come first.
  1383  			}
  1384  			if fe.hasInterior {
  1385  				t.testEdge(fe.shapeID, fe.edge)
  1386  			}
  1387  		}
  1388  	}
  1389  
  1390  	// Save the state of the edges being removed, so that it can be restored
  1391  	// when we are finished processing this cell and its children.  We don't
  1392  	// need to save the state of the edges being added because they aren't being
  1393  	// removed from "edges" and will therefore be updated normally as we visit
  1394  	// this cell and its children.
  1395  	t.saveAndClearStateBefore(s.pendingAdditionsPos)
  1396  
  1397  	// Create a faceEdge for each edge in this cell that isn't being removed.
  1398  	var faceEdges []*faceEdge
  1399  	trackerMoved := false
  1400  
  1401  	cell := iter.IndexCell()
  1402  	for _, clipped := range cell.shapes {
  1403  		shapeID := clipped.shapeID
  1404  		shape := s.Shape(shapeID)
  1405  		if shape == nil {
  1406  			continue // This shape is being removed.
  1407  		}
  1408  
  1409  		numClipped := clipped.numEdges()
  1410  
  1411  		// If this shape has an interior, start tracking whether we are inside the
  1412  		// shape. updateEdges wants to know whether the entry vertex of this
  1413  		// cell is inside the shape, but we only know whether the center of the
  1414  		// cell is inside the shape, so we need to test all the edges against the
  1415  		// line segment from the cell center to the entry vertex.
  1416  		edge := &faceEdge{
  1417  			shapeID:     shapeID,
  1418  			hasInterior: shape.Dimension() == 2,
  1419  		}
  1420  
  1421  		if edge.hasInterior {
  1422  			t.addShape(shapeID, clipped.containsCenter)
  1423  			// There might not be any edges in this entire cell (i.e., it might be
  1424  			// in the interior of all shapes), so we delay updating the tracker
  1425  			// until we see the first edge.
  1426  			if !trackerMoved && numClipped > 0 {
  1427  				t.moveTo(p.Center())
  1428  				t.drawTo(p.EntryVertex())
  1429  				t.setNextCellID(p.id)
  1430  				trackerMoved = true
  1431  			}
  1432  		}
  1433  		for i := 0; i < numClipped; i++ {
  1434  			edgeID := clipped.edges[i]
  1435  			edge.edgeID = edgeID
  1436  			edge.edge = shape.Edge(edgeID)
  1437  			edge.MaxLevel = maxLevelForEdge(edge.edge)
  1438  			if edge.hasInterior {
  1439  				t.testEdge(shapeID, edge.edge)
  1440  			}
  1441  			var ok bool
  1442  			edge.a, edge.b, ok = ClipToPaddedFace(edge.edge.V0, edge.edge.V1, p.id.Face(), cellPadding)
  1443  			if !ok {
  1444  				panic("invariant failure in ShapeIndex")
  1445  			}
  1446  			faceEdges = append(faceEdges, edge)
  1447  		}
  1448  	}
  1449  	// Now create a clippedEdge for each faceEdge, and put them in "new_edges".
  1450  	var newEdges []*clippedEdge
  1451  	for _, faceEdge := range faceEdges {
  1452  		clipped := &clippedEdge{
  1453  			faceEdge: faceEdge,
  1454  			bound:    clippedEdgeBound(faceEdge.a, faceEdge.b, p.bound),
  1455  		}
  1456  		newEdges = append(newEdges, clipped)
  1457  	}
  1458  
  1459  	// Discard any edges from "edges" that are being removed, and append the
  1460  	// remainder to "newEdges"  (This keeps the edges sorted by shape id.)
  1461  	for i, clipped := range edges {
  1462  		if !s.isShapeBeingRemoved(clipped.faceEdge.shapeID) {
  1463  			newEdges = append(newEdges, edges[i:]...)
  1464  			break
  1465  		}
  1466  	}
  1467  
  1468  	// Update the edge list and delete this cell from the index.
  1469  	edges, newEdges = newEdges, edges
  1470  	delete(s.cellMap, p.id)
  1471  	// TODO(roberts): delete from s.Cells
  1472  }
  1473  
  1474  // testAllEdges calls the trackers testEdge on all edges from shapes that have interiors.
  1475  func (s *ShapeIndex) testAllEdges(edges []*clippedEdge, t *tracker) {
  1476  	for _, edge := range edges {
  1477  		if edge.faceEdge.hasInterior {
  1478  			t.testEdge(edge.faceEdge.shapeID, edge.faceEdge.edge)
  1479  		}
  1480  	}
  1481  }
  1482  
  1483  // countShapes reports the number of distinct shapes that are either associated with the
  1484  // given edges, or that are currently stored in the InteriorTracker.
  1485  func (s *ShapeIndex) countShapes(edges []*clippedEdge, shapeIDs []int32) int {
  1486  	count := 0
  1487  	lastShapeID := int32(-1)
  1488  
  1489  	// next clipped shape id in the shapeIDs list.
  1490  	clippedNext := int32(0)
  1491  	// index of the current element in the shapeIDs list.
  1492  	shapeIDidx := 0
  1493  	for _, edge := range edges {
  1494  		if edge.faceEdge.shapeID == lastShapeID {
  1495  			continue
  1496  		}
  1497  
  1498  		count++
  1499  		lastShapeID = edge.faceEdge.shapeID
  1500  
  1501  		// Skip over any containing shapes up to and including this one,
  1502  		// updating count as appropriate.
  1503  		for ; shapeIDidx < len(shapeIDs); shapeIDidx++ {
  1504  			clippedNext = shapeIDs[shapeIDidx]
  1505  			if clippedNext > lastShapeID {
  1506  				break
  1507  			}
  1508  			if clippedNext < lastShapeID {
  1509  				count++
  1510  			}
  1511  		}
  1512  	}
  1513  
  1514  	// Count any remaining containing shapes.
  1515  	count += len(shapeIDs) - shapeIDidx
  1516  	return count
  1517  }
  1518  
  1519  // maxLevelForEdge reports the maximum level for a given edge.
  1520  func maxLevelForEdge(edge Edge) int {
  1521  	// Compute the maximum cell size for which this edge is considered long.
  1522  	// The calculation does not need to be perfectly accurate, so we use Norm
  1523  	// rather than Angle for speed.
  1524  	cellSize := edge.V0.Sub(edge.V1.Vector).Norm() * cellSizeToLongEdgeRatio
  1525  	// Now return the first level encountered during subdivision where the
  1526  	// average cell size is at most cellSize.
  1527  	return AvgEdgeMetric.MinLevel(cellSize)
  1528  }
  1529  
  1530  // removeShapeInternal does the actual work for removing a given shape from the index.
  1531  func (s *ShapeIndex) removeShapeInternal(removed *removedShape, allEdges [][]faceEdge, t *tracker) {
  1532  	// TODO(roberts): finish the implementation of this.
  1533  }
  1534  

View as plain text