crossing_edge_query.go

Documentation: github.com/golang/geo/s2

     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  }
   410
View as plain text