edge_query.go (30347B)
1 // Copyright 2019 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 "sort" 19 20 "github.com/golang/geo/s1" 21 ) 22 23 // EdgeQueryOptions holds the options for controlling how EdgeQuery operates. 24 // 25 // Options can be chained together builder-style: 26 // 27 // opts = NewClosestEdgeQueryOptions(). 28 // MaxResults(1). 29 // DistanceLimit(s1.ChordAngleFromAngle(3 * s1.Degree)). 30 // MaxError(s1.ChordAngleFromAngle(0.001 * s1.Degree)) 31 // query = NewClosestEdgeQuery(index, opts) 32 // 33 // or set individually: 34 // 35 // opts = NewClosestEdgeQueryOptions() 36 // opts.IncludeInteriors(true) 37 // 38 // or just inline: 39 // 40 // query = NewClosestEdgeQuery(index, NewClosestEdgeQueryOptions().MaxResults(3)) 41 // 42 // If you pass a nil as the options you get the default values for the options. 43 type EdgeQueryOptions struct { 44 common *queryOptions 45 } 46 47 // DistanceLimit specifies that only edges whose distance to the target is 48 // within, this distance should be returned. Edges whose distance is equal 49 // are not returned. To include values that are equal, specify the limit with 50 // the next largest representable distance. i.e. limit.Successor(). 51 func (e *EdgeQueryOptions) DistanceLimit(limit s1.ChordAngle) *EdgeQueryOptions { 52 e.common = e.common.DistanceLimit(limit) 53 return e 54 } 55 56 // IncludeInteriors specifies whether polygon interiors should be 57 // included when measuring distances. 58 func (e *EdgeQueryOptions) IncludeInteriors(x bool) *EdgeQueryOptions { 59 e.common = e.common.IncludeInteriors(x) 60 return e 61 } 62 63 // UseBruteForce sets or disables the use of brute force in a query. 64 func (e *EdgeQueryOptions) UseBruteForce(x bool) *EdgeQueryOptions { 65 e.common = e.common.UseBruteForce(x) 66 return e 67 } 68 69 // MaxError specifies that edges up to dist away than the true 70 // matching edges may be substituted in the result set, as long as such 71 // edges satisfy all the remaining search criteria (such as DistanceLimit). 72 // This option only has an effect if MaxResults is also specified; 73 // otherwise all edges closer than MaxDistance will always be returned. 74 func (e *EdgeQueryOptions) MaxError(dist s1.ChordAngle) *EdgeQueryOptions { 75 e.common = e.common.MaxError(dist) 76 return e 77 } 78 79 // MaxResults specifies that at most MaxResults edges should be returned. 80 // This must be at least 1. 81 func (e *EdgeQueryOptions) MaxResults(n int) *EdgeQueryOptions { 82 e.common = e.common.MaxResults(n) 83 return e 84 } 85 86 // NewClosestEdgeQueryOptions returns a set of edge query options suitable 87 // for performing closest edge queries. 88 func NewClosestEdgeQueryOptions() *EdgeQueryOptions { 89 return &EdgeQueryOptions{ 90 common: newQueryOptions(minDistance(0)), 91 } 92 } 93 94 // NewFurthestEdgeQueryOptions returns a set of edge query options suitable 95 // for performing furthest edge queries. 96 func NewFurthestEdgeQueryOptions() *EdgeQueryOptions { 97 return &EdgeQueryOptions{ 98 common: newQueryOptions(maxDistance(0)), 99 } 100 } 101 102 // EdgeQueryResult represents an edge that meets the target criteria for the 103 // query. Note the following special cases: 104 // 105 // - ShapeID >= 0 && EdgeID < 0 represents the interior of a shape. 106 // Such results may be returned when the option IncludeInteriors is true. 107 // 108 // - ShapeID < 0 && EdgeID < 0 is returned to indicate that no edge 109 // satisfies the requested query options. 110 type EdgeQueryResult struct { 111 distance distance 112 shapeID int32 113 edgeID int32 114 } 115 116 // Distance reports the distance between the edge in this shape that satisfied 117 // the query's parameters. 118 func (e EdgeQueryResult) Distance() s1.ChordAngle { return e.distance.chordAngle() } 119 120 // ShapeID reports the ID of the Shape this result is for. 121 func (e EdgeQueryResult) ShapeID() int32 { return e.shapeID } 122 123 // EdgeID reports the ID of the edge in the results Shape. 124 func (e EdgeQueryResult) EdgeID() int32 { return e.edgeID } 125 126 // newEdgeQueryResult returns a result instance with default values. 127 func newEdgeQueryResult(target distanceTarget) EdgeQueryResult { 128 return EdgeQueryResult{ 129 distance: target.distance().infinity(), 130 shapeID: -1, 131 edgeID: -1, 132 } 133 } 134 135 // IsInterior reports if this result represents the interior of a Shape. 136 func (e EdgeQueryResult) IsInterior() bool { 137 return e.shapeID >= 0 && e.edgeID < 0 138 } 139 140 // IsEmpty reports if this has no edge that satisfies the given edge query options. 141 // This result is only returned in one special case, namely when FindEdge() does 142 // not find any suitable edges. 143 func (e EdgeQueryResult) IsEmpty() bool { 144 return e.shapeID < 0 145 } 146 147 // Less reports if this results is less that the other first by distance, 148 // then by (shapeID, edgeID). This is used for sorting. 149 func (e EdgeQueryResult) Less(other EdgeQueryResult) bool { 150 if e.distance.chordAngle() != other.distance.chordAngle() { 151 return e.distance.less(other.distance) 152 } 153 if e.shapeID != other.shapeID { 154 return e.shapeID < other.shapeID 155 } 156 return e.edgeID < other.edgeID 157 } 158 159 // EdgeQuery is used to find the edge(s) between two geometries that match a 160 // given set of options. It is flexible enough so that it can be adapted to 161 // compute maximum distances and even potentially Hausdorff distances. 162 // 163 // By using the appropriate options, this type can answer questions such as: 164 // 165 // - Find the minimum distance between two geometries A and B. 166 // - Find all edges of geometry A that are within a distance D of geometry B. 167 // - Find the k edges of geometry A that are closest to a given point P. 168 // 169 // You can also specify whether polygons should include their interiors (i.e., 170 // if a point is contained by a polygon, should the distance be zero or should 171 // it be measured to the polygon boundary?) 172 // 173 // The input geometries may consist of any number of points, polylines, and 174 // polygons (collectively referred to as "shapes"). Shapes do not need to be 175 // disjoint; they may overlap or intersect arbitrarily. The implementation is 176 // designed to be fast for both simple and complex geometries. 177 type EdgeQuery struct { 178 index *ShapeIndex 179 opts *queryOptions 180 target distanceTarget 181 182 // True if opts.maxError must be subtracted from ShapeIndex cell distances 183 // in order to ensure that such distances are measured conservatively. This 184 // is true only if the target takes advantage of maxError in order to 185 // return faster results, and 0 < maxError < distanceLimit. 186 useConservativeCellDistance bool 187 188 // The decision about whether to use the brute force algorithm is based on 189 // counting the total number of edges in the index. However if the index 190 // contains a large number of shapes, this in itself might take too long. 191 // So instead we only count edges up to (maxBruteForceIndexSize() + 1) 192 // for the current target type (stored as indexNumEdgesLimit). 193 indexNumEdges int 194 indexNumEdgesLimit int 195 196 // The distance beyond which we can safely ignore further candidate edges. 197 // (Candidates that are exactly at the limit are ignored; this is more 198 // efficient for UpdateMinDistance and should not affect clients since 199 // distance measurements have a small amount of error anyway.) 200 // 201 // Initially this is the same as the maximum distance specified by the user, 202 // but it can also be updated by the algorithm (see maybeAddResult). 203 distanceLimit distance 204 205 // The current set of results of the query. 206 results []EdgeQueryResult 207 208 // This field is true when duplicates must be avoided explicitly. This 209 // is achieved by maintaining a separate set keyed by (shapeID, edgeID) 210 // only, and checking whether each edge is in that set before computing the 211 // distance to it. 212 avoidDuplicates bool 213 214 // testedEdges tracks the set of shape and edges that have already been tested. 215 testedEdges map[ShapeEdgeID]uint32 216 217 // For the optimized algorihm we precompute the top-level CellIDs that 218 // will be added to the priority queue. There can be at most 6 of these 219 // cells. Essentially this is just a covering of the indexed edges, except 220 // that we also store pointers to the corresponding ShapeIndexCells to 221 // reduce the number of index seeks required. 222 indexCovering []CellID 223 indexCells []*ShapeIndexCell 224 225 // The algorithm maintains a priority queue of unprocessed CellIDs, sorted 226 // in increasing order of distance from the target. 227 queue *queryQueue 228 229 iter *ShapeIndexIterator 230 maxDistanceCovering []CellID 231 initialCells []CellID 232 } 233 234 // NewClosestEdgeQuery returns an EdgeQuery that is used for finding the 235 // closest edge(s) to a given Point, Edge, Cell, or geometry collection. 236 // 237 // You can find either the k closest edges, or all edges within a given 238 // radius, or both (i.e., the k closest edges up to a given maximum radius). 239 // E.g. to find all the edges within 5 kilometers, set the DistanceLimit in 240 // the options. 241 // 242 // By default *all* edges are returned, so you should always specify either 243 // MaxResults or DistanceLimit options or both. 244 // 245 // Note that by default, distances are measured to the boundary and interior 246 // of polygons. For example, if a point is inside a polygon then its distance 247 // is zero. To change this behavior, set the IncludeInteriors option to false. 248 // 249 // If you only need to test whether the distance is above or below a given 250 // threshold (e.g., 10 km), you can use the IsDistanceLess() method. This is 251 // much faster than actually calculating the distance with FindEdge, 252 // since the implementation can stop as soon as it can prove that the minimum 253 // distance is either above or below the threshold. 254 func NewClosestEdgeQuery(index *ShapeIndex, opts *EdgeQueryOptions) *EdgeQuery { 255 if opts == nil { 256 opts = NewClosestEdgeQueryOptions() 257 } 258 e := &EdgeQuery{ 259 testedEdges: make(map[ShapeEdgeID]uint32), 260 index: index, 261 opts: opts.common, 262 queue: newQueryQueue(), 263 } 264 265 return e 266 } 267 268 // NewFurthestEdgeQuery returns an EdgeQuery that is used for finding the 269 // furthest edge(s) to a given Point, Edge, Cell, or geometry collection. 270 // 271 // The furthest edge is defined as the one which maximizes the 272 // distance from any point on that edge to any point on the target geometry. 273 // 274 // Similar to the example in NewClosestEdgeQuery, to find the 5 furthest edges 275 // from a given Point: 276 func NewFurthestEdgeQuery(index *ShapeIndex, opts *EdgeQueryOptions) *EdgeQuery { 277 if opts == nil { 278 opts = NewFurthestEdgeQueryOptions() 279 } 280 e := &EdgeQuery{ 281 testedEdges: make(map[ShapeEdgeID]uint32), 282 index: index, 283 opts: opts.common, 284 queue: newQueryQueue(), 285 } 286 287 return e 288 } 289 290 // Reset resets the state of this EdgeQuery. 291 func (e *EdgeQuery) Reset() { 292 e.indexNumEdges = 0 293 e.indexNumEdgesLimit = 0 294 e.indexCovering = nil 295 e.indexCells = nil 296 } 297 298 // FindEdges returns the edges for the given target that satisfy the current options. 299 // 300 // Note that if opts.IncludeInteriors is true, the results may include some 301 // entries with edge_id == -1. This indicates that the target intersects 302 // the indexed polygon with the given ShapeID. 303 func (e *EdgeQuery) FindEdges(target distanceTarget) []EdgeQueryResult { 304 return e.findEdges(target, e.opts) 305 } 306 307 // Distance reports the distance to the target. If the index or target is empty, 308 // returns the EdgeQuery's maximal sentinel. 309 // 310 // Use IsDistanceLess()/IsDistanceGreater() if you only want to compare the 311 // distance against a threshold value, since it is often much faster. 312 func (e *EdgeQuery) Distance(target distanceTarget) s1.ChordAngle { 313 return e.findEdge(target, e.opts).Distance() 314 } 315 316 // IsDistanceLess reports if the distance to target is less than the given limit. 317 // 318 // This method is usually much faster than Distance(), since it is much 319 // less work to determine whether the minimum distance is above or below a 320 // threshold than it is to calculate the actual minimum distance. 321 // 322 // If you wish to check if the distance is less than or equal to the limit, use: 323 // 324 // query.IsDistanceLess(target, limit.Successor()) 325 // 326 func (e *EdgeQuery) IsDistanceLess(target distanceTarget, limit s1.ChordAngle) bool { 327 opts := e.opts 328 opts = opts.MaxResults(1). 329 DistanceLimit(limit). 330 MaxError(s1.StraightChordAngle) 331 return !e.findEdge(target, opts).IsEmpty() 332 } 333 334 // IsDistanceGreater reports if the distance to target is greater than limit. 335 // 336 // This method is usually much faster than Distance, since it is much 337 // less work to determine whether the maximum distance is above or below a 338 // threshold than it is to calculate the actual maximum distance. 339 // If you wish to check if the distance is less than or equal to the limit, use: 340 // 341 // query.IsDistanceGreater(target, limit.Predecessor()) 342 // 343 func (e *EdgeQuery) IsDistanceGreater(target distanceTarget, limit s1.ChordAngle) bool { 344 return e.IsDistanceLess(target, limit) 345 } 346 347 // IsConservativeDistanceLessOrEqual reports if the distance to target is less 348 // or equal to the limit, where the limit has been expanded by the maximum error 349 // for the distance calculation. 350 // 351 // For example, suppose that we want to test whether two geometries might 352 // intersect each other after they are snapped together using Builder 353 // (using the IdentitySnapFunction with a given "snap radius"). Since 354 // Builder uses exact distance predicates (s2predicates), we need to 355 // measure the distance between the two geometries conservatively. If the 356 // distance is definitely greater than "snap radius", then the geometries 357 // are guaranteed to not intersect after snapping. 358 func (e *EdgeQuery) IsConservativeDistanceLessOrEqual(target distanceTarget, limit s1.ChordAngle) bool { 359 return e.IsDistanceLess(target, limit.Expanded(minUpdateDistanceMaxError(limit))) 360 } 361 362 // IsConservativeDistanceGreaterOrEqual reports if the distance to the target is greater 363 // than or equal to the given limit with some small tolerance. 364 func (e *EdgeQuery) IsConservativeDistanceGreaterOrEqual(target distanceTarget, limit s1.ChordAngle) bool { 365 return e.IsDistanceGreater(target, limit.Expanded(-minUpdateDistanceMaxError(limit))) 366 } 367 368 // findEdges returns the closest edges to the given target that satisfy the given options. 369 // 370 // Note that if opts.includeInteriors is true, the results may include some 371 // entries with edgeID == -1. This indicates that the target intersects the 372 // indexed polygon with the given shapeID. 373 func (e *EdgeQuery) findEdges(target distanceTarget, opts *queryOptions) []EdgeQueryResult { 374 e.findEdgesInternal(target, opts) 375 // TODO(roberts): Revisit this if there is a heap or other sorted and 376 // uniquing datastructure we can use instead of just a slice. 377 e.results = sortAndUniqueResults(e.results) 378 if len(e.results) > e.opts.maxResults { 379 e.results = e.results[:e.opts.maxResults] 380 } 381 return e.results 382 } 383 384 func sortAndUniqueResults(results []EdgeQueryResult) []EdgeQueryResult { 385 if len(results) <= 1 { 386 return results 387 } 388 sort.Slice(results, func(i, j int) bool { return results[i].Less(results[j]) }) 389 j := 0 390 for i := 1; i < len(results); i++ { 391 if results[j] == results[i] { 392 continue 393 } 394 j++ 395 results[j] = results[i] 396 } 397 return results[:j+1] 398 } 399 400 // findEdge is a convenience method that returns exactly one edge, and if no 401 // edges satisfy the given search criteria, then a default Result is returned. 402 // 403 // This is primarily to ease the usage of a number of the methods in the DistanceTargets 404 // and in EdgeQuery. 405 func (e *EdgeQuery) findEdge(target distanceTarget, opts *queryOptions) EdgeQueryResult { 406 opts.MaxResults(1) 407 e.findEdges(target, opts) 408 if len(e.results) > 0 { 409 return e.results[0] 410 } 411 412 return newEdgeQueryResult(target) 413 } 414 415 // findEdgesInternal does the actual work for find edges that match the given options. 416 func (e *EdgeQuery) findEdgesInternal(target distanceTarget, opts *queryOptions) { 417 e.target = target 418 e.opts = opts 419 420 e.testedEdges = make(map[ShapeEdgeID]uint32) 421 e.distanceLimit = target.distance().fromChordAngle(opts.distanceLimit) 422 e.results = make([]EdgeQueryResult, 0) 423 424 if e.distanceLimit == target.distance().zero() { 425 return 426 } 427 428 if opts.includeInteriors { 429 shapeIDs := map[int32]struct{}{} 430 e.target.visitContainingShapes(e.index, func(containingShape Shape, targetPoint Point) bool { 431 shapeIDs[e.index.idForShape(containingShape)] = struct{}{} 432 return len(shapeIDs) < opts.maxResults 433 }) 434 for shapeID := range shapeIDs { 435 e.addResult(EdgeQueryResult{target.distance().zero(), shapeID, -1}) 436 } 437 438 if e.distanceLimit == target.distance().zero() { 439 return 440 } 441 } 442 443 // If maxError > 0 and the target takes advantage of this, then we may 444 // need to adjust the distance estimates to ShapeIndex cells to ensure 445 // that they are always a lower bound on the true distance. For example, 446 // suppose max_distance == 100, maxError == 30, and we compute the distance 447 // to the target from some cell C0 as d(C0) == 80. Then because the target 448 // takes advantage of maxError, the true distance could be as low as 50. 449 // In order not to miss edges contained by such cells, we need to subtract 450 // maxError from the distance estimates. This behavior is controlled by 451 // the useConservativeCellDistance flag. 452 // 453 // However there is one important case where this adjustment is not 454 // necessary, namely when distanceLimit < maxError, This is because 455 // maxError only affects the algorithm once at least maxEdges edges 456 // have been found that satisfy the given distance limit. At that point, 457 // maxError is subtracted from distanceLimit in order to ensure that 458 // any further matches are closer by at least that amount. But when 459 // distanceLimit < maxError, this reduces the distance limit to 0, 460 // i.e. all remaining candidate cells and edges can safely be discarded. 461 // (This is how IsDistanceLess() and friends are implemented.) 462 targetUsesMaxError := opts.maxError != target.distance().zero().chordAngle() && 463 e.target.setMaxError(opts.maxError) 464 465 // Note that we can't compare maxError and distanceLimit directly 466 // because one is a Delta and one is a Distance. Instead we subtract them. 467 e.useConservativeCellDistance = targetUsesMaxError && 468 (e.distanceLimit == target.distance().infinity() || 469 target.distance().zero().less(e.distanceLimit.sub(target.distance().fromChordAngle(opts.maxError)))) 470 471 // Use the brute force algorithm if the index is small enough. To avoid 472 // spending too much time counting edges when there are many shapes, we stop 473 // counting once there are too many edges. We may need to recount the edges 474 // if we later see a target with a larger brute force edge threshold. 475 minOptimizedEdges := e.target.maxBruteForceIndexSize() + 1 476 if minOptimizedEdges > e.indexNumEdgesLimit && e.indexNumEdges >= e.indexNumEdgesLimit { 477 e.indexNumEdges = e.index.NumEdgesUpTo(minOptimizedEdges) 478 e.indexNumEdgesLimit = minOptimizedEdges 479 } 480 481 if opts.useBruteForce || e.indexNumEdges < minOptimizedEdges { 482 // The brute force algorithm already considers each edge exactly once. 483 e.avoidDuplicates = false 484 e.findEdgesBruteForce() 485 } else { 486 // If the target takes advantage of maxError then we need to avoid 487 // duplicate edges explicitly. (Otherwise it happens automatically.) 488 e.avoidDuplicates = targetUsesMaxError && opts.maxResults > 1 489 e.findEdgesOptimized() 490 } 491 } 492 493 func (e *EdgeQuery) addResult(r EdgeQueryResult) { 494 e.results = append(e.results, r) 495 if e.opts.maxResults == 1 { 496 // Optimization for the common case where only the closest edge is wanted. 497 e.distanceLimit = r.distance.sub(e.target.distance().fromChordAngle(e.opts.maxError)) 498 } 499 // TODO(roberts): Add the other if/else cases when a different data structure 500 // is used for the results. 501 } 502 503 func (e *EdgeQuery) maybeAddResult(shape Shape, edgeID int32) { 504 if _, ok := e.testedEdges[ShapeEdgeID{e.index.idForShape(shape), edgeID}]; e.avoidDuplicates && !ok { 505 return 506 } 507 edge := shape.Edge(int(edgeID)) 508 dist := e.distanceLimit 509 510 if dist, ok := e.target.updateDistanceToEdge(edge, dist); ok { 511 e.addResult(EdgeQueryResult{dist, e.index.idForShape(shape), edgeID}) 512 } 513 } 514 515 func (e *EdgeQuery) findEdgesBruteForce() { 516 // Range over all shapes in the index. Does order matter here? if so 517 // switch to for i = 0 .. n? 518 for _, shape := range e.index.shapes { 519 // TODO(roberts): can this happen if we are only ranging over current entries? 520 if shape == nil { 521 continue 522 } 523 for edgeID := int32(0); edgeID < int32(shape.NumEdges()); edgeID++ { 524 e.maybeAddResult(shape, edgeID) 525 } 526 } 527 } 528 529 func (e *EdgeQuery) findEdgesOptimized() { 530 e.initQueue() 531 // Repeatedly find the closest Cell to "target" and either split it into 532 // its four children or process all of its edges. 533 for e.queue.size() > 0 { 534 // We need to copy the top entry before removing it, and we need to 535 // remove it before adding any new entries to the queue. 536 entry := e.queue.pop() 537 538 if !entry.distance.less(e.distanceLimit) { 539 e.queue.reset() // Clear any remaining entries. 540 break 541 } 542 // If this is already known to be an index cell, just process it. 543 if entry.indexCell != nil { 544 e.processEdges(entry) 545 continue 546 } 547 // Otherwise split the cell into its four children. Before adding a 548 // child back to the queue, we first check whether it is empty. We do 549 // this in two seek operations rather than four by seeking to the key 550 // between children 0 and 1 and to the key between children 2 and 3. 551 id := entry.id 552 ch := id.Children() 553 e.iter.seek(ch[1].RangeMin()) 554 555 if !e.iter.Done() && e.iter.CellID() <= ch[1].RangeMax() { 556 e.processOrEnqueueCell(ch[1]) 557 } 558 if e.iter.Prev() && e.iter.CellID() >= id.RangeMin() { 559 e.processOrEnqueueCell(ch[0]) 560 } 561 562 e.iter.seek(ch[3].RangeMin()) 563 if !e.iter.Done() && e.iter.CellID() <= id.RangeMax() { 564 e.processOrEnqueueCell(ch[3]) 565 } 566 if e.iter.Prev() && e.iter.CellID() >= ch[2].RangeMin() { 567 e.processOrEnqueueCell(ch[2]) 568 } 569 } 570 } 571 572 func (e *EdgeQuery) processOrEnqueueCell(id CellID) { 573 if e.iter.CellID() == id { 574 e.processOrEnqueue(id, e.iter.IndexCell()) 575 } else { 576 e.processOrEnqueue(id, nil) 577 } 578 } 579 580 func (e *EdgeQuery) initQueue() { 581 if len(e.indexCovering) == 0 { 582 // We delay iterator initialization until now to make queries on very 583 // small indexes a bit faster (i.e., where brute force is used). 584 e.iter = NewShapeIndexIterator(e.index) 585 } 586 587 // Optimization: if the user is searching for just the closest edge, and the 588 // center of the target's bounding cap happens to intersect an index cell, 589 // then we try to limit the search region to a small disc by first 590 // processing the edges in that cell. This sets distance_limit_ based on 591 // the closest edge in that cell, which we can then use to limit the search 592 // area. This means that the cell containing "target" will be processed 593 // twice, but in general this is still faster. 594 // 595 // TODO(roberts): Even if the cap center is not contained, we could still 596 // process one or both of the adjacent index cells in CellID order, 597 // provided that those cells are closer than distanceLimit. 598 cb := e.target.capBound() 599 if cb.IsEmpty() { 600 return // Empty target. 601 } 602 603 if e.opts.maxResults == 1 && e.iter.LocatePoint(cb.Center()) { 604 e.processEdges(&queryQueueEntry{ 605 distance: e.target.distance().zero(), 606 id: e.iter.CellID(), 607 indexCell: e.iter.IndexCell(), 608 }) 609 // Skip the rest of the algorithm if we found an intersecting edge. 610 if e.distanceLimit == e.target.distance().zero() { 611 return 612 } 613 } 614 if len(e.indexCovering) == 0 { 615 e.initCovering() 616 } 617 if e.distanceLimit == e.target.distance().infinity() { 618 // Start with the precomputed index covering. 619 for i := range e.indexCovering { 620 e.processOrEnqueue(e.indexCovering[i], e.indexCells[i]) 621 } 622 } else { 623 // Compute a covering of the search disc and intersect it with the 624 // precomputed index covering. 625 coverer := &RegionCoverer{MaxCells: 4, LevelMod: 1, MaxLevel: maxLevel} 626 627 radius := cb.Radius() + e.distanceLimit.chordAngleBound().Angle() 628 searchCB := CapFromCenterAngle(cb.Center(), radius) 629 maxDistCover := coverer.FastCovering(searchCB) 630 e.initialCells = CellUnionFromIntersection(e.indexCovering, maxDistCover) 631 632 // Now we need to clean up the initial cells to ensure that they all 633 // contain at least one cell of the ShapeIndex. (Some may not intersect 634 // the index at all, while other may be descendants of an index cell.) 635 i, j := 0, 0 636 for i < len(e.initialCells) { 637 idI := e.initialCells[i] 638 // Find the top-level cell that contains this initial cell. 639 for e.indexCovering[j].RangeMax() < idI { 640 j++ 641 } 642 643 idJ := e.indexCovering[j] 644 if idI == idJ { 645 // This initial cell is one of the top-level cells. Use the 646 // precomputed ShapeIndexCell pointer to avoid an index seek. 647 e.processOrEnqueue(idJ, e.indexCells[j]) 648 i++ 649 j++ 650 } else { 651 // This initial cell is a proper descendant of a top-level cell. 652 // Check how it is related to the cells of the ShapeIndex. 653 r := e.iter.LocateCellID(idI) 654 if r == Indexed { 655 // This cell is a descendant of an index cell. 656 // Enqueue it and skip any other initial cells 657 // that are also descendants of this cell. 658 e.processOrEnqueue(e.iter.CellID(), e.iter.IndexCell()) 659 lastID := e.iter.CellID().RangeMax() 660 for i < len(e.initialCells) && e.initialCells[i] <= lastID { 661 i++ 662 } 663 } else { 664 // Enqueue the cell only if it contains at least one index cell. 665 if r == Subdivided { 666 e.processOrEnqueue(idI, nil) 667 } 668 i++ 669 } 670 } 671 } 672 } 673 } 674 675 func (e *EdgeQuery) initCovering() { 676 // Find the range of Cells spanned by the index and choose a level such 677 // that the entire index can be covered with just a few cells. These are 678 // the "top-level" cells. There are two cases: 679 // 680 // - If the index spans more than one face, then there is one top-level cell 681 // per spanned face, just big enough to cover the index cells on that face. 682 // 683 // - If the index spans only one face, then we find the smallest cell "C" 684 // that covers the index cells on that face (just like the case above). 685 // Then for each of the 4 children of "C", if the child contains any index 686 // cells then we create a top-level cell that is big enough to just fit 687 // those index cells (i.e., shrinking the child as much as possible to fit 688 // its contents). This essentially replicates what would happen if we 689 // started with "C" as the top-level cell, since "C" would immediately be 690 // split, except that we take the time to prune the children further since 691 // this will save work on every subsequent query. 692 e.indexCovering = make([]CellID, 0, 6) 693 694 // TODO(roberts): Use a single iterator below and save position 695 // information using pair {CellID, ShapeIndexCell}. 696 next := NewShapeIndexIterator(e.index, IteratorBegin) 697 last := NewShapeIndexIterator(e.index, IteratorEnd) 698 last.Prev() 699 if next.CellID() != last.CellID() { 700 // The index has at least two cells. Choose a level such that the entire 701 // index can be spanned with at most 6 cells (if the index spans multiple 702 // faces) or 4 cells (it the index spans a single face). 703 level, ok := next.CellID().CommonAncestorLevel(last.CellID()) 704 if !ok { 705 level = 0 706 } else { 707 level++ 708 } 709 710 // Visit each potential top-level cell except the last (handled below). 711 lastID := last.CellID().Parent(level) 712 for id := next.CellID().Parent(level); id != lastID; id = id.Next() { 713 // Skip any top-level cells that don't contain any index cells. 714 if id.RangeMax() < next.CellID() { 715 continue 716 } 717 718 // Find the range of index cells contained by this top-level cell and 719 // then shrink the cell if necessary so that it just covers them. 720 cellFirst := next.clone() 721 next.seek(id.RangeMax().Next()) 722 cellLast := next.clone() 723 cellLast.Prev() 724 e.addInitialRange(cellFirst, cellLast) 725 break 726 } 727 728 } 729 e.addInitialRange(next, last) 730 } 731 732 // addInitialRange adds an entry to the indexCovering and indexCells that covers the given 733 // inclusive range of cells. 734 // 735 // This requires that first and last cells have a common ancestor. 736 func (e *EdgeQuery) addInitialRange(first, last *ShapeIndexIterator) { 737 if first.CellID() == last.CellID() { 738 // The range consists of a single index cell. 739 e.indexCovering = append(e.indexCovering, first.CellID()) 740 e.indexCells = append(e.indexCells, first.IndexCell()) 741 } else { 742 // Add the lowest common ancestor of the given range. 743 level, _ := first.CellID().CommonAncestorLevel(last.CellID()) 744 e.indexCovering = append(e.indexCovering, first.CellID().Parent(level)) 745 e.indexCells = append(e.indexCells, nil) 746 } 747 } 748 749 // processEdges processes all the edges of the given index cell. 750 func (e *EdgeQuery) processEdges(entry *queryQueueEntry) { 751 for _, clipped := range entry.indexCell.shapes { 752 shape := e.index.Shape(clipped.shapeID) 753 for j := 0; j < clipped.numEdges(); j++ { 754 e.maybeAddResult(shape, int32(clipped.edges[j])) 755 } 756 } 757 } 758 759 // processOrEnqueue the given cell id and indexCell. 760 func (e *EdgeQuery) processOrEnqueue(id CellID, indexCell *ShapeIndexCell) { 761 if indexCell != nil { 762 // If this index cell has only a few edges, then it is faster to check 763 // them directly rather than computing the minimum distance to the Cell 764 // and inserting it into the queue. 765 const minEdgesToEnqueue = 10 766 numEdges := indexCell.numEdges() 767 if numEdges == 0 { 768 return 769 } 770 if numEdges < minEdgesToEnqueue { 771 // Set "distance" to zero to avoid the expense of computing it. 772 e.processEdges(&queryQueueEntry{ 773 distance: e.target.distance().zero(), 774 id: id, 775 indexCell: indexCell, 776 }) 777 return 778 } 779 } 780 781 // Otherwise compute the minimum distance to any point in the cell and add 782 // it to the priority queue. 783 cell := CellFromCellID(id) 784 dist := e.distanceLimit 785 var ok bool 786 if dist, ok = e.target.updateDistanceToCell(cell, dist); !ok { 787 return 788 } 789 if e.useConservativeCellDistance { 790 // Ensure that "distance" is a lower bound on the true distance to the cell. 791 dist = dist.sub(e.target.distance().fromChordAngle(e.opts.maxError)) 792 } 793 794 e.queue.push(&queryQueueEntry{ 795 distance: dist, 796 id: id, 797 indexCell: indexCell, 798 }) 799 } 800 801 // TODO(roberts): Remaining pieces 802 // GetEdge 803 // Project