/* * Copyright 2006 Google Inc. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package com.google.common.geometry; import com.google.common.base.Preconditions; import com.google.common.collect.Lists; import com.google.common.collect.Sets; import java.util.ArrayList; import java.util.Arrays; import java.util.Comparator; import java.util.HashSet; import java.util.List; import java.util.Set; public abstract strictfp class S2EdgeIndex { /** * Thicken the edge in all directions by roughly 1% of the edge length when * thickenEdge is true. */ private static final double THICKENING = 0.01; /** * Threshold for small angles, that help lenientCrossing to determine whether * two edges are likely to intersect. */ private static final double MAX_DET_ERROR = 1e-14; /** * The cell containing each edge, as given in the parallel array * edges. */ private long[] cells; /** * The edge contained by each cell, as given in the parallel array * cells. */ private int[] edges; /** * No cell strictly below this level appears in mapping. Initially leaf level, * that's the minimum level at which we will ever look for test edges. */ private int minimumS2LevelUsed; /** * Has the index been computed already? */ private boolean indexComputed; /** * Number of queries so far */ private int queryCount; /** * Empties the index in case it already contained something. */ public void reset() { minimumS2LevelUsed = S2CellId.MAX_LEVEL; indexComputed = false; queryCount = 0; cells = null; edges = null; } /** * Compares [cell1, edge1] to [cell2, edge2], by cell first and edge second. * * @return -1 if [cell1, edge1] is less than [cell2, edge2], 1 if [cell1, * edge1] is greater than [cell2, edge2], 0 otherwise. */ private static final int compare(long cell1, int edge1, long cell2, int edge2) { if (cell1 < cell2) { return -1; } else if (cell1 > cell2) { return 1; } else if (edge1 < edge2) { return -1; } else if (edge1 > edge2) { return 1; } else { return 0; } } /** Computes the index (if it has not been previously done). */ public final void computeIndex() { if (indexComputed) { return; } List cellList = Lists.newArrayList(); List edgeList = Lists.newArrayList(); for (int i = 0; i < getNumEdges(); ++i) { S2Point from = edgeFrom(i); S2Point to = edgeTo(i); ArrayList cover = Lists.newArrayList(); int level = getCovering(from, to, true, cover); minimumS2LevelUsed = Math.min(minimumS2LevelUsed, level); for (S2CellId cellId : cover) { cellList.add(cellId.id()); edgeList.add(i); } } cells = new long[cellList.size()]; edges = new int[edgeList.size()]; for (int i = 0; i < cells.length; i++) { cells[i] = cellList.get(i); edges[i] = edgeList.get(i); } sortIndex(); indexComputed = true; } /** Sorts the parallel cells and edges arrays. */ private void sortIndex() { // create an array of indices and sort based on the values in the parallel // arrays at each index Integer[] indices = new Integer[cells.length]; for (int i = 0; i < indices.length; i++) { indices[i] = i; } Arrays.sort(indices, new Comparator() { @Override public int compare(Integer index1, Integer index2) { return S2EdgeIndex.compare(cells[index1], edges[index1], cells[index2], edges[index2]); } }); // copy the cells and edges in the order given by the sorted list of indices long[] newCells = new long[cells.length]; int[] newEdges = new int[edges.length]; for (int i = 0; i < indices.length; i++) { newCells[i] = cells[indices[i]]; newEdges[i] = edges[indices[i]]; } // replace the cells and edges with the sorted arrays cells = newCells; edges = newEdges; } public final boolean isIndexComputed() { return indexComputed; } /** * Tell the index that we just received a new request for candidates. Useful * to compute when to switch to quad tree. */ protected final void incrementQueryCount() { ++queryCount; } /** * If the index hasn't been computed yet, looks at how much work has gone into * iterating using the brute force method, and how much more work is planned * as defined by 'cost'. If it were to have been cheaper to use a quad tree * from the beginning, then compute it now. This guarantees that we will never * use more than twice the time we would have used had we known in advance * exactly how many edges we would have wanted to test. It is the theoretical * best. * * The value 'n' is the number of iterators we expect to request from this * edge index. * * If we have m data edges and n query edges, then the brute force cost is m * * n * testCost where testCost is taken to be the cost of * EdgeCrosser.robustCrossing, measured to be about 30ns at the time of this * writing. * * If we compute the index, the cost becomes: m * costInsert + n * * costFind(m) * * - costInsert can be expected to be reasonably stable, and was measured at * 1200ns with the BM_QuadEdgeInsertionCost benchmark. * * - costFind depends on the length of the edge . For m=1000 edges, we got * timings ranging from 1ms (edge the length of the polygon) to 40ms. The * latter is for very long query edges, and needs to be optimized. We will * assume for the rest of the discussion that costFind is roughly 3ms. * * When doing one additional query, the differential cost is m * testCost - * costFind(m) With the numbers above, it is better to use the quad tree (if * we have it) if m >= 100. * * If m = 100, 30 queries will give m*n*testCost = m * costInsert = 100ms, * while the marginal cost to find is 3ms. Thus, this is a reasonable thing to * do. */ public final void predictAdditionalCalls(int n) { if (indexComputed) { return; } if (getNumEdges() > 100 && (queryCount + n) > 30) { computeIndex(); } } /** * Overwrite these functions to give access to the underlying data. The * function getNumEdges() returns the number of edges in the index, while * edgeFrom(index) and edgeTo(index) return the "from" and "to" endpoints of * the edge at the given index. */ protected abstract int getNumEdges(); protected abstract S2Point edgeFrom(int index); protected abstract S2Point edgeTo(int index); /** * Appends to "candidateCrossings" all edge references which may cross the * given edge. This is done by covering the edge and then finding all * references of edges whose coverings overlap this covering. Parent cells are * checked level by level. Child cells are checked all at once by taking * advantage of the natural ordering of S2CellIds. */ protected void findCandidateCrossings(S2Point a, S2Point b, List candidateCrossings) { Preconditions.checkState(indexComputed); ArrayList cover = Lists.newArrayList(); getCovering(a, b, false, cover); // Edge references are inserted into the map once for each covering cell, so // absorb duplicates here Set uniqueSet = new HashSet(); getEdgesInParentCells(cover, uniqueSet); // TODO(user): An important optimization for long query // edges (Contains queries): keep a bounding cap and clip the query // edge to the cap before starting the descent. getEdgesInChildrenCells(a, b, cover, uniqueSet); candidateCrossings.clear(); candidateCrossings.addAll(uniqueSet); } /** * Returns the smallest cell containing all four points, or * {@link S2CellId#sentinel()} if they are not all on the same face. The * points don't need to be normalized. */ private static S2CellId containingCell(S2Point pa, S2Point pb, S2Point pc, S2Point pd) { S2CellId a = S2CellId.fromPoint(pa); S2CellId b = S2CellId.fromPoint(pb); S2CellId c = S2CellId.fromPoint(pc); S2CellId d = S2CellId.fromPoint(pd); if (a.face() != b.face() || a.face() != c.face() || a.face() != d.face()) { return S2CellId.sentinel(); } while (!a.equals(b) || !a.equals(c) || !a.equals(d)) { a = a.parent(); b = b.parent(); c = c.parent(); d = d.parent(); } return a; } /** * Returns the smallest cell containing both points, or Sentinel if they are * not all on the same face. The points don't need to be normalized. */ private static S2CellId containingCell(S2Point pa, S2Point pb) { S2CellId a = S2CellId.fromPoint(pa); S2CellId b = S2CellId.fromPoint(pb); if (a.face() != b.face()) { return S2CellId.sentinel(); } while (!a.equals(b)) { a = a.parent(); b = b.parent(); } return a; } /** * Computes a cell covering of an edge. Clears edgeCovering and returns the * level of the s2 cells used in the covering (only one level is ever used for * each call). * * If thickenEdge is true, the edge is thickened and extended by 1% of its * length. * * It is guaranteed that no child of a covering cell will fully contain the * covered edge. */ private int getCovering( S2Point a, S2Point b, boolean thickenEdge, ArrayList edgeCovering) { edgeCovering.clear(); // Selects the ideal s2 level at which to cover the edge, this will be the // level whose S2 cells have a width roughly commensurate to the length of // the edge. We multiply the edge length by 2*THICKENING to guarantee the // thickening is honored (it's not a big deal if we honor it when we don't // request it) when doing the covering-by-cap trick. double edgeLength = a.angle(b); int idealLevel = S2Projections.MIN_WIDTH.getMaxLevel(edgeLength * (1 + 2 * THICKENING)); S2CellId containingCellId; if (!thickenEdge) { containingCellId = containingCell(a, b); } else { if (idealLevel == S2CellId.MAX_LEVEL) { // If the edge is tiny, instabilities are more likely, so we // want to limit the number of operations. // We pretend we are in a cell much larger so as to trigger the // 'needs covering' case, so we won't try to thicken the edge. containingCellId = (new S2CellId(0xFFF0)).parent(3); } else { S2Point pq = S2Point.mul(S2Point.minus(b, a), THICKENING); S2Point ortho = S2Point.mul(S2Point.normalize(S2Point.crossProd(pq, a)), edgeLength * THICKENING); S2Point p = S2Point.minus(a, pq); S2Point q = S2Point.add(b, pq); // If p and q were antipodal, the edge wouldn't be lengthened, // and it could even flip! This is not a problem because // idealLevel != 0 here. The farther p and q can be is roughly // a quarter Earth away from each other, so we remain // Theta(THICKENING). containingCellId = containingCell(S2Point.minus(p, ortho), S2Point.add(p, ortho), S2Point.minus(q, ortho), S2Point.add(q, ortho)); } } // Best case: edge is fully contained in a cell that's not too big. if (!containingCellId.equals(S2CellId.sentinel()) && containingCellId.level() >= idealLevel - 2) { edgeCovering.add(containingCellId); return containingCellId.level(); } if (idealLevel == 0) { // Edge is very long, maybe even longer than a face width, so the // trick below doesn't work. For now, we will add the whole S2 sphere. // TODO(user): Do something a tad smarter (and beware of the // antipodal case). for (S2CellId cellid = S2CellId.begin(0); !cellid.equals(S2CellId.end(0)); cellid = cellid.next()) { edgeCovering.add(cellid); } return 0; } // TODO(user): Check trick below works even when vertex is at // interface // between three faces. // Use trick as in S2PolygonBuilder.PointIndex.findNearbyPoint: // Cover the edge by a cap centered at the edge midpoint, then cover // the cap by four big-enough cells around the cell vertex closest to the // cap center. S2Point middle = S2Point.normalize(S2Point.div(S2Point.add(a, b), 2)); int actualLevel = Math.min(idealLevel, S2CellId.MAX_LEVEL - 1); S2CellId.fromPoint(middle).getVertexNeighbors(actualLevel, edgeCovering); return actualLevel; } /** * Filters a list of entries down to the inclusive range defined by the given * cells, in O(log N) time. * * @param cell1 One side of the inclusive query range. * @param cell2 The other side of the inclusive query range. * @return An array of length 2, containing the start/end indices. */ private int[] getEdges(long cell1, long cell2) { // ensure cell1 <= cell2 if (cell1 > cell2) { long temp = cell1; cell1 = cell2; cell2 = temp; } // The binary search returns -N-1 to indicate an insertion point at index N, // if an exact match cannot be found. Since the edge indices queried for are // not valid edge indices, we will always get -N-1, so we immediately // convert to N. return new int[]{ -1 - binarySearch(cell1, Integer.MIN_VALUE), -1 - binarySearch(cell2, Integer.MAX_VALUE)}; } private int binarySearch(long cell, int edge) { int low = 0; int high = cells.length - 1; while (low <= high) { int mid = (low + high) >>> 1; int cmp = compare(cells[mid], edges[mid], cell, edge); if (cmp < 0) { low = mid + 1; } else if (cmp > 0) { high = mid - 1; } else { return mid; } } return -(low + 1); } /** * Adds to candidateCrossings all the edges present in any ancestor of any * cell of cover, down to minimumS2LevelUsed. The cell->edge map is in the * variable mapping. */ private void getEdgesInParentCells(List cover, Set candidateCrossings) { // Find all parent cells of covering cells. Set parentCells = Sets.newHashSet(); for (S2CellId coverCell : cover) { for (int parentLevel = coverCell.level() - 1; parentLevel >= minimumS2LevelUsed; --parentLevel) { if (!parentCells.add(coverCell.parent(parentLevel))) { break; // cell is already in => parents are too. } } } // Put parent cell edge references into result. for (S2CellId parentCell : parentCells) { int[] bounds = getEdges(parentCell.id(), parentCell.id()); for (int i = bounds[0]; i < bounds[1]; i++) { candidateCrossings.add(edges[i]); } } } /** * Returns true if ab possibly crosses cd, by clipping tiny angles to zero. */ private static boolean lenientCrossing(S2Point a, S2Point b, S2Point c, S2Point d) { // assert (S2.isUnitLength(a)); // assert (S2.isUnitLength(b)); // assert (S2.isUnitLength(c)); double acb = S2Point.crossProd(a, c).dotProd(b); double bda = S2Point.crossProd(b, d).dotProd(a); if (Math.abs(acb) < MAX_DET_ERROR || Math.abs(bda) < MAX_DET_ERROR) { return true; } if (acb * bda < 0) { return false; } double cbd = S2Point.crossProd(c, b).dotProd(d); double dac = S2Point.crossProd(c, a).dotProd(c); if (Math.abs(cbd) < MAX_DET_ERROR || Math.abs(dac) < MAX_DET_ERROR) { return true; } return (acb * cbd >= 0) && (acb * dac >= 0); } /** * Returns true if the edge and the cell (including boundary) intersect. */ private static boolean edgeIntersectsCellBoundary(S2Point a, S2Point b, S2Cell cell) { S2Point[] vertices = new S2Point[4]; for (int i = 0; i < 4; ++i) { vertices[i] = cell.getVertex(i); } for (int i = 0; i < 4; ++i) { S2Point fromPoint = vertices[i]; S2Point toPoint = vertices[(i + 1) % 4]; if (lenientCrossing(a, b, fromPoint, toPoint)) { return true; } } return false; } /** * Appends to candidateCrossings the edges that are fully contained in an S2 * covering of edge. The covering of edge used is initially cover, but is * refined to eliminate quickly subcells that contain many edges but do not * intersect with edge. */ private void getEdgesInChildrenCells(S2Point a, S2Point b, List cover, Set candidateCrossings) { // Put all edge references of (covering cells + descendant cells) into // result. // This relies on the natural ordering of S2CellIds. S2Cell[] children = null; while (!cover.isEmpty()) { S2CellId cell = cover.remove(cover.size() - 1); int[] bounds = getEdges(cell.rangeMin().id(), cell.rangeMax().id()); if (bounds[1] - bounds[0] <= 16) { for (int i = bounds[0]; i < bounds[1]; i++) { candidateCrossings.add(edges[i]); } } else { // Add cells at this level bounds = getEdges(cell.id(), cell.id()); for (int i = bounds[0]; i < bounds[1]; i++) { candidateCrossings.add(edges[i]); } // Recurse on the children -- hopefully some will be empty. if (children == null) { children = new S2Cell[4]; for (int i = 0; i < 4; ++i) { children[i] = new S2Cell(); } } new S2Cell(cell).subdivide(children); for (S2Cell child : children) { // TODO(user): Do the check for the four cells at once, // as it is enough to check the four edges between the cells. At // this time, we are checking 16 edges, 4 times too many. // // Note that given the guarantee of AppendCovering, it is enough // to check that the edge intersect with the cell boundary as it // cannot be fully contained in a cell. if (edgeIntersectsCellBoundary(a, b, child)) { cover.add(child.id()); } } } } } /* * An iterator on data edges that may cross a query edge (a,b). Create the * iterator, call getCandidates(), then hasNext()/next() repeatedly. * * The current edge in the iteration has index index(), goes between from() * and to(). */ public static class DataEdgeIterator { /** * The structure containing the data edges. */ private final S2EdgeIndex edgeIndex; /** * Tells whether getCandidates() obtained the candidates through brute force * iteration or using the quad tree structure. */ private boolean isBruteForce; /** * Index of the current edge and of the edge before the last next() call. */ private int currentIndex; /** * Cache of edgeIndex.getNumEdges() so that hasNext() doesn't make an extra * call */ private int numEdges; /** * All the candidates obtained by getCandidates() when we are using a * quad-tree (i.e. isBruteForce = false). */ ArrayList candidates; /** * Index within array above. We have: currentIndex = * candidates.get(currentIndexInCandidates). */ private int currentIndexInCandidates; public DataEdgeIterator(S2EdgeIndex edgeIndex) { this.edgeIndex = edgeIndex; candidates = Lists.newArrayList(); } /** * Initializes the iterator to iterate over a set of candidates that may * cross the edge (a,b). */ public void getCandidates(S2Point a, S2Point b) { edgeIndex.predictAdditionalCalls(1); isBruteForce = !edgeIndex.isIndexComputed(); if (isBruteForce) { edgeIndex.incrementQueryCount(); currentIndex = 0; numEdges = edgeIndex.getNumEdges(); } else { candidates.clear(); edgeIndex.findCandidateCrossings(a, b, candidates); currentIndexInCandidates = 0; if (!candidates.isEmpty()) { currentIndex = candidates.get(0); } } } /** * Index of the current edge in the iteration. */ public int index() { Preconditions.checkState(hasNext()); return currentIndex; } /** * False if there are no more candidates; true otherwise. */ public boolean hasNext() { if (isBruteForce) { return (currentIndex < numEdges); } else { return currentIndexInCandidates < candidates.size(); } } /** * Iterate to the next available candidate. */ public void next() { Preconditions.checkState(hasNext()); if (isBruteForce) { ++currentIndex; } else { ++currentIndexInCandidates; if (currentIndexInCandidates < candidates.size()) { currentIndex = candidates.get(currentIndexInCandidates); } } } } }