crossing_edge_query.go (13363B)
1 // Copyright 2017 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/r2" 21 ) 22 23 // CrossingEdgeQuery is used to find the Edge IDs of Shapes that are crossed by 24 // a given edge(s). 25 // 26 // Note that if you need to query many edges, it is more efficient to declare 27 // a single CrossingEdgeQuery instance and reuse it. 28 // 29 // If you want to find *all* the pairs of crossing edges, it is more efficient to 30 // use the not yet implemented VisitCrossings in shapeutil. 31 type CrossingEdgeQuery struct { 32 index *ShapeIndex 33 34 // temporary values used while processing a query. 35 a, b r2.Point 36 iter *ShapeIndexIterator 37 38 // candidate cells generated when finding crossings. 39 cells []*ShapeIndexCell 40 } 41 42 // NewCrossingEdgeQuery creates a CrossingEdgeQuery for the given index. 43 func NewCrossingEdgeQuery(index *ShapeIndex) *CrossingEdgeQuery { 44 c := &CrossingEdgeQuery{ 45 index: index, 46 iter: index.Iterator(), 47 } 48 return c 49 } 50 51 // Crossings returns the set of edge of the shape S that intersect the given edge AB. 52 // If the CrossingType is Interior, then only intersections at a point interior to both 53 // edges are reported, while if it is CrossingTypeAll then edges that share a vertex 54 // are also reported. 55 func (c *CrossingEdgeQuery) Crossings(a, b Point, shape Shape, crossType CrossingType) []int { 56 edges := c.candidates(a, b, shape) 57 if len(edges) == 0 { 58 return nil 59 } 60 61 crosser := NewEdgeCrosser(a, b) 62 out := 0 63 n := len(edges) 64 65 for in := 0; in < n; in++ { 66 b := shape.Edge(edges[in]) 67 sign := crosser.CrossingSign(b.V0, b.V1) 68 if crossType == CrossingTypeAll && (sign == MaybeCross || sign == Cross) || crossType != CrossingTypeAll && sign == Cross { 69 edges[out] = edges[in] 70 out++ 71 } 72 } 73 74 if out < n { 75 edges = edges[0:out] 76 } 77 return edges 78 } 79 80 // EdgeMap stores a sorted set of edge ids for each shape. 81 type EdgeMap map[Shape][]int 82 83 // CrossingsEdgeMap returns the set of all edges in the index that intersect the given 84 // edge AB. If crossType is CrossingTypeInterior, then only intersections at a 85 // point interior to both edges are reported, while if it is CrossingTypeAll 86 // then edges that share a vertex are also reported. 87 // 88 // The edges are returned as a mapping from shape to the edges of that shape 89 // that intersect AB. Every returned shape has at least one crossing edge. 90 func (c *CrossingEdgeQuery) CrossingsEdgeMap(a, b Point, crossType CrossingType) EdgeMap { 91 edgeMap := c.candidatesEdgeMap(a, b) 92 if len(edgeMap) == 0 { 93 return nil 94 } 95 96 crosser := NewEdgeCrosser(a, b) 97 for shape, edges := range edgeMap { 98 out := 0 99 n := len(edges) 100 for in := 0; in < n; in++ { 101 edge := shape.Edge(edges[in]) 102 sign := crosser.CrossingSign(edge.V0, edge.V1) 103 if (crossType == CrossingTypeAll && (sign == MaybeCross || sign == Cross)) || (crossType != CrossingTypeAll && sign == Cross) { 104 edgeMap[shape][out] = edges[in] 105 out++ 106 } 107 } 108 109 if out == 0 { 110 delete(edgeMap, shape) 111 } else { 112 if out < n { 113 edgeMap[shape] = edgeMap[shape][0:out] 114 } 115 } 116 } 117 return edgeMap 118 } 119 120 // candidates returns a superset of the edges of the given shape that intersect 121 // the edge AB. 122 func (c *CrossingEdgeQuery) candidates(a, b Point, shape Shape) []int { 123 var edges []int 124 125 // For small loops it is faster to use brute force. The threshold below was 126 // determined using benchmarks. 127 const maxBruteForceEdges = 27 128 maxEdges := shape.NumEdges() 129 if maxEdges <= maxBruteForceEdges { 130 edges = make([]int, maxEdges) 131 for i := 0; i < maxEdges; i++ { 132 edges[i] = i 133 } 134 return edges 135 } 136 137 // Compute the set of index cells intersected by the query edge. 138 c.getCellsForEdge(a, b) 139 if len(c.cells) == 0 { 140 return nil 141 } 142 143 // Gather all the edges that intersect those cells and sort them. 144 // TODO(roberts): Shapes don't track their ID, so we need to range over 145 // the index to find the ID manually. 146 var shapeID int32 147 for k, v := range c.index.shapes { 148 if v == shape { 149 shapeID = k 150 } 151 } 152 153 for _, cell := range c.cells { 154 if cell == nil { 155 continue 156 } 157 clipped := cell.findByShapeID(shapeID) 158 if clipped == nil { 159 continue 160 } 161 edges = append(edges, clipped.edges...) 162 } 163 164 if len(c.cells) > 1 { 165 edges = uniqueInts(edges) 166 } 167 168 return edges 169 } 170 171 // uniqueInts returns the sorted uniqued values from the given input. 172 func uniqueInts(in []int) []int { 173 var edges []int 174 m := make(map[int]bool) 175 for _, i := range in { 176 if m[i] { 177 continue 178 } 179 m[i] = true 180 edges = append(edges, i) 181 } 182 sort.Ints(edges) 183 return edges 184 } 185 186 // candidatesEdgeMap returns a map from shapes to the superse of edges for that 187 // shape that intersect the edge AB. 188 // 189 // CAVEAT: This method may return shapes that have an empty set of candidate edges. 190 // However the return value is non-empty only if at least one shape has a candidate edge. 191 func (c *CrossingEdgeQuery) candidatesEdgeMap(a, b Point) EdgeMap { 192 edgeMap := make(EdgeMap) 193 194 // If there are only a few edges then it's faster to use brute force. We 195 // only bother with this optimization when there is a single shape. 196 if len(c.index.shapes) == 1 { 197 // Typically this method is called many times, so it is worth checking 198 // whether the edge map is empty or already consists of a single entry for 199 // this shape, and skip clearing edge map in that case. 200 shape := c.index.Shape(0) 201 202 // Note that we leave the edge map non-empty even if there are no candidates 203 // (i.e., there is a single entry with an empty set of edges). 204 edgeMap[shape] = c.candidates(a, b, shape) 205 return edgeMap 206 } 207 208 // Compute the set of index cells intersected by the query edge. 209 c.getCellsForEdge(a, b) 210 if len(c.cells) == 0 { 211 return edgeMap 212 } 213 214 // Gather all the edges that intersect those cells and sort them. 215 for _, cell := range c.cells { 216 for _, clipped := range cell.shapes { 217 s := c.index.Shape(clipped.shapeID) 218 for j := 0; j < clipped.numEdges(); j++ { 219 edgeMap[s] = append(edgeMap[s], clipped.edges[j]) 220 } 221 } 222 } 223 224 if len(c.cells) > 1 { 225 for s, edges := range edgeMap { 226 edgeMap[s] = uniqueInts(edges) 227 } 228 } 229 230 return edgeMap 231 } 232 233 // getCells returns the set of ShapeIndexCells that might contain edges intersecting 234 // the edge AB in the given cell root. This method is used primarily by loop and shapeutil. 235 func (c *CrossingEdgeQuery) getCells(a, b Point, root *PaddedCell) []*ShapeIndexCell { 236 aUV, bUV, ok := ClipToFace(a, b, root.id.Face()) 237 if ok { 238 c.a = aUV 239 c.b = bUV 240 edgeBound := r2.RectFromPoints(c.a, c.b) 241 if root.Bound().Intersects(edgeBound) { 242 c.computeCellsIntersected(root, edgeBound) 243 } 244 } 245 246 if len(c.cells) == 0 { 247 return nil 248 } 249 250 return c.cells 251 } 252 253 // getCellsForEdge populates the cells field to the set of index cells intersected by an edge AB. 254 func (c *CrossingEdgeQuery) getCellsForEdge(a, b Point) { 255 c.cells = nil 256 257 segments := FaceSegments(a, b) 258 for _, segment := range segments { 259 c.a = segment.a 260 c.b = segment.b 261 262 // Optimization: rather than always starting the recursive subdivision at 263 // the top level face cell, instead we start at the smallest S2CellId that 264 // contains the edge (the edge root cell). This typically lets us skip 265 // quite a few levels of recursion since most edges are short. 266 edgeBound := r2.RectFromPoints(c.a, c.b) 267 pcell := PaddedCellFromCellID(CellIDFromFace(segment.face), 0) 268 edgeRoot := pcell.ShrinkToFit(edgeBound) 269 270 // Now we need to determine how the edge root cell is related to the cells 271 // in the spatial index (cellMap). There are three cases: 272 // 273 // 1. edgeRoot is an index cell or is contained within an index cell. 274 // In this case we only need to look at the contents of that cell. 275 // 2. edgeRoot is subdivided into one or more index cells. In this case 276 // we recursively subdivide to find the cells intersected by AB. 277 // 3. edgeRoot does not intersect any index cells. In this case there 278 // is nothing to do. 279 relation := c.iter.LocateCellID(edgeRoot) 280 if relation == Indexed { 281 // edgeRoot is an index cell or is contained by an index cell (case 1). 282 c.cells = append(c.cells, c.iter.IndexCell()) 283 } else if relation == Subdivided { 284 // edgeRoot is subdivided into one or more index cells (case 2). We 285 // find the cells intersected by AB using recursive subdivision. 286 if !edgeRoot.isFace() { 287 pcell = PaddedCellFromCellID(edgeRoot, 0) 288 } 289 c.computeCellsIntersected(pcell, edgeBound) 290 } 291 } 292 } 293 294 // computeCellsIntersected computes the index cells intersected by the current 295 // edge that are descendants of pcell and adds them to this queries set of cells. 296 func (c *CrossingEdgeQuery) computeCellsIntersected(pcell *PaddedCell, edgeBound r2.Rect) { 297 298 c.iter.seek(pcell.id.RangeMin()) 299 if c.iter.Done() || c.iter.CellID() > pcell.id.RangeMax() { 300 // The index does not contain pcell or any of its descendants. 301 return 302 } 303 if c.iter.CellID() == pcell.id { 304 // The index contains this cell exactly. 305 c.cells = append(c.cells, c.iter.IndexCell()) 306 return 307 } 308 309 // Otherwise, split the edge among the four children of pcell. 310 center := pcell.Middle().Lo() 311 312 if edgeBound.X.Hi < center.X { 313 // Edge is entirely contained in the two left children. 314 c.clipVAxis(edgeBound, center.Y, 0, pcell) 315 return 316 } else if edgeBound.X.Lo >= center.X { 317 // Edge is entirely contained in the two right children. 318 c.clipVAxis(edgeBound, center.Y, 1, pcell) 319 return 320 } 321 322 childBounds := c.splitUBound(edgeBound, center.X) 323 if edgeBound.Y.Hi < center.Y { 324 // Edge is entirely contained in the two lower children. 325 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, 0, 0), childBounds[0]) 326 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, 1, 0), childBounds[1]) 327 } else if edgeBound.Y.Lo >= center.Y { 328 // Edge is entirely contained in the two upper children. 329 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, 0, 1), childBounds[0]) 330 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, 1, 1), childBounds[1]) 331 } else { 332 // The edge bound spans all four children. The edge itself intersects 333 // at most three children (since no padding is being used). 334 c.clipVAxis(childBounds[0], center.Y, 0, pcell) 335 c.clipVAxis(childBounds[1], center.Y, 1, pcell) 336 } 337 } 338 339 // clipVAxis computes the intersected cells recursively for a given padded cell. 340 // Given either the left (i=0) or right (i=1) side of a padded cell pcell, 341 // determine whether the current edge intersects the lower child, upper child, 342 // or both children, and call c.computeCellsIntersected recursively on those children. 343 // The center is the v-coordinate at the center of pcell. 344 func (c *CrossingEdgeQuery) clipVAxis(edgeBound r2.Rect, center float64, i int, pcell *PaddedCell) { 345 if edgeBound.Y.Hi < center { 346 // Edge is entirely contained in the lower child. 347 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, i, 0), edgeBound) 348 } else if edgeBound.Y.Lo >= center { 349 // Edge is entirely contained in the upper child. 350 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, i, 1), edgeBound) 351 } else { 352 // The edge intersects both children. 353 childBounds := c.splitVBound(edgeBound, center) 354 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, i, 0), childBounds[0]) 355 c.computeCellsIntersected(PaddedCellFromParentIJ(pcell, i, 1), childBounds[1]) 356 } 357 } 358 359 // splitUBound returns the bound for two children as a result of spliting the 360 // current edge at the given value U. 361 func (c *CrossingEdgeQuery) splitUBound(edgeBound r2.Rect, u float64) [2]r2.Rect { 362 v := edgeBound.Y.ClampPoint(interpolateFloat64(u, c.a.X, c.b.X, c.a.Y, c.b.Y)) 363 // diag indicates which diagonal of the bounding box is spanned by AB: 364 // it is 0 if AB has positive slope, and 1 if AB has negative slope. 365 var diag int 366 if (c.a.X > c.b.X) != (c.a.Y > c.b.Y) { 367 diag = 1 368 } 369 return splitBound(edgeBound, 0, diag, u, v) 370 } 371 372 // splitVBound returns the bound for two children as a result of spliting the 373 // current edge into two child edges at the given value V. 374 func (c *CrossingEdgeQuery) splitVBound(edgeBound r2.Rect, v float64) [2]r2.Rect { 375 u := edgeBound.X.ClampPoint(interpolateFloat64(v, c.a.Y, c.b.Y, c.a.X, c.b.X)) 376 var diag int 377 if (c.a.X > c.b.X) != (c.a.Y > c.b.Y) { 378 diag = 1 379 } 380 return splitBound(edgeBound, diag, 0, u, v) 381 } 382 383 // splitBound returns the bounds for the two childrenn as a result of spliting 384 // the current edge into two child edges at the given point (u,v). uEnd and vEnd 385 // indicate which bound endpoints of the first child will be updated. 386 func splitBound(edgeBound r2.Rect, uEnd, vEnd int, u, v float64) [2]r2.Rect { 387 var childBounds = [2]r2.Rect{ 388 edgeBound, 389 edgeBound, 390 } 391 392 if uEnd == 1 { 393 childBounds[0].X.Lo = u 394 childBounds[1].X.Hi = u 395 } else { 396 childBounds[0].X.Hi = u 397 childBounds[1].X.Lo = u 398 } 399 400 if vEnd == 1 { 401 childBounds[0].Y.Lo = v 402 childBounds[1].Y.Hi = v 403 } else { 404 childBounds[0].Y.Hi = v 405 childBounds[1].Y.Lo = v 406 } 407 408 return childBounds 409 }