diff --git a/s2/incident_edge_tracker.go b/s2/incident_edge_tracker.go new file mode 100644 index 0000000..3dbaaf8 --- /dev/null +++ b/s2/incident_edge_tracker.go @@ -0,0 +1,173 @@ +// Copyright 2025 The S2 Geometry Project Authors. All rights reserved. +// +// 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 s2 + +// incidentEdgeKey is a tuple of (shape id, vertex) that compares by shape id. +type incidentEdgeKey struct { + shapeID int32 + vertex Point +} + +// We need a strict ordering to be a valid key for an ordered container, but +// we don't actually care about the ordering of the vertices (as long as +// they're grouped by shape id). Vertices are 3D points so they don't have a +// natural ordering, so we'll just compare them lexicographically. +func (i incidentEdgeKey) Cmp(o incidentEdgeKey) int { + if i.shapeID < o.shapeID { + return -1 + } + if i.shapeID > o.shapeID { + return 1 + } + + return i.vertex.Cmp(o.vertex.Vector) +} + +// vertexEdge is a tuple of vertex and edgeID for processing incident edges. +type vertexEdge struct { + vertex Point + edgeID int32 +} + +// incidentEdgeTracker is a used for detecting and tracking shape edges that +// are incident on the same vertex. Edges of multiple shapes may be tracked, +// but lookup is by shape id and vertex: there is no facility to get all +// edges of all shapes at a vertex. Edge vertices must compare exactly equal +// to be considered the same vertex, no tolerance is applied as this isn't +// intended for e.g.: snapping shapes together, which Builder does better +// and more robustly. +// +// To use, instantiate and then add edges with one or more sequences of calls, +// where each sequence begins with startShape(), followed by addEdge() calls to +// add edges for that shape, and ends with finishShape(). Those sequences do +// not need to visit shapes or edges in order. Then, call incidentEdges() to get +// the resulting map from incidentEdgeKeys (which are shapeId, vertex pairs) to +// a set of edgeIds of the shape that are incident to that vertex.. +// +// This works on a block of edges at a time, meaning that to detect incident +// edges on a particular vertex, we must have at least three edges incident +// at that vertex when finishShape() is called. We don't maintain partial +// information between calls. However, subject to this constraint, a single +// shape's edges may be defined with multiple sequences of startShape(), +// addEdge()... , finishShape() calls. +// +// The reason for this is simple: most edges don't have more than two incident +// edges (the incoming and outgoing edge). If we had to maintain information +// between calls, we'd end up with a map that contains every vertex, to no +// benefit. Instead, when finishShape() is called, we discard vertices that +// contain two or fewer incident edges. +// +// In principle this isn't a real limitation because generally we process a +// ShapeIndex cell at a time, and, if a vertex has multiple edges, we'll see +// all the edges in the same cell as the vertex, and, in general, it's possible +// to aggregate edges before calling. +// +// The tracker maintains incident edges until it's cleared. If you call it with +// each cell from an ShapeIndex, then at the end you will have all the +// incident edge information for the whole index. If only a subset is needed, +// call reset() when you're done. +type incidentEdgeTracker struct { + currentShapeID int32 + + nursery []vertexEdge + + // We can and do encounter the same edges multiple times, so we need to + // deduplicate edges as they're inserted. + edgeMap map[incidentEdgeKey]map[int32]bool +} + +// newIncidentEdgeTracker returns a new tracker. +func newIncidentEdgeTracker() *incidentEdgeTracker { + return &incidentEdgeTracker{ + currentShapeID: -1, + nursery: []vertexEdge{}, + edgeMap: make(map[incidentEdgeKey]map[int32]bool), + } +} + +// startShape is used to start adding edges to the edge tracker. After calling, +// any vertices with multiple (> 2) incident edges will appear in the +// incident edge map. +func (t *incidentEdgeTracker) startShape(id int32) { + t.currentShapeID = id + t.nursery = t.nursery[:0] +} + +// addEdge adds the given edges start to the nursery, and if not degenerate, +// adds it second endpoint as well. +func (t *incidentEdgeTracker) addEdge(edgeID int32, e Edge) { + if t.currentShapeID < 0 { + return + } + + // Add non-degenerate edges to the nursery. + t.nursery = append(t.nursery, vertexEdge{vertex: e.V0, edgeID: edgeID}) + if !e.IsDegenerate() { + t.nursery = append(t.nursery, vertexEdge{vertex: e.V1, edgeID: edgeID}) + } +} + +func (t *incidentEdgeTracker) finishShape() { + // We want to keep any vertices with more than two incident edges. We could + // sort the array by vertex and remove any with fewer, but that would require + // shifting the array and could turn quadratic quickly. + // + // Instead we'll scan forward from each vertex, swapping entries with the same + // vertex into a contiguous range. Once we've done all the swapping we can + // just make sure that we have at least three edges in the range. + nurserySize := len(t.nursery) + for start := 0; start < nurserySize; { + end := start + 1 + + // Scan to the end of the array, swap entries so that entries with + // the same vertex as the start are adjacent. + next := start + currVertex := t.nursery[start].vertex + for next+1 < nurserySize { + next++ + if t.nursery[next].vertex == currVertex { + t.nursery[next], t.nursery[end] = t.nursery[end], t.nursery[next] + end++ + } + } + + // Most vertices will have two incident edges (the incoming edge and the + // outgoing edge), which aren't interesting, skip them. + numEdges := end - start + if numEdges <= 2 { + start = end + continue + } + + key := incidentEdgeKey{ + shapeID: t.currentShapeID, + vertex: t.nursery[start].vertex, + } + + // If we don't have this key yet, create it manually. + if _, ok := t.edgeMap[key]; !ok { + t.edgeMap[key] = map[int32]bool{} + } + + for ; start != end; start++ { + t.edgeMap[key][t.nursery[start].edgeID] = true + } + } +} + +// reset removes all incident edges from the tracker. +func (t *incidentEdgeTracker) reset() { + t.edgeMap = make(map[incidentEdgeKey]map[int32]bool) +} diff --git a/s2/incident_edge_tracker_test.go b/s2/incident_edge_tracker_test.go new file mode 100644 index 0000000..a80232b --- /dev/null +++ b/s2/incident_edge_tracker_test.go @@ -0,0 +1,70 @@ +// Copyright 2025 The S2 Geometry Project Authors. All rights reserved. +// +// 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 s2 + +import ( + "testing" +) + +func TestIncidentEdgeTrackerBasic(t *testing.T) { + tests := []struct { + index string + want int + }{ + // These shapeindex strings came from validation query's test + // corpus to determine which ones ended up actually getting + // tracked edges. + { + // Has 0 tracked edges + index: "## 0:0, 1:1", + want: 0, + }, + { + // Has 1 tracked edges + index: "## 2:0, 0:-2, -2:0, 0:2; 2:0, 0:-1, -1:0, 0:1", + want: 1, + }, + { + // Has 2 tracked edges + index: "## 2:0, 0:-2, -2:0, 0:2; 2:0, 0:-1, -2:0, 0:1", + want: 2, + }, + } + + for _, test := range tests { + index := makeShapeIndex(test.index) + index.Build() + + iter := index.Iterator() + celldata := newIndexCellData() + celldata.loadCell(index, iter.CellID(), iter.IndexCell()) + + tracker := newIncidentEdgeTracker() + + for _, clipped := range celldata.indexCell.shapes { + shapeID := clipped.shapeID + tracker.startShape(shapeID) + for _, e := range celldata.shapeEdges(shapeID) { + tracker.addEdge(e.ID, e.Edge) + } + tracker.finishShape() + } + + if got := len(tracker.edgeMap); got != test.want { + t.Errorf("incidentEdgeTracker should have %d edges, got %d", + test.want, got) + } + } +} diff --git a/s2/index_cell_data.go b/s2/index_cell_data.go new file mode 100644 index 0000000..8a41ede --- /dev/null +++ b/s2/index_cell_data.go @@ -0,0 +1,328 @@ +// Copyright 2025 The S2 Geometry Project Authors. All rights reserved. +// +// 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 s2 + +import ( + "sync" + "sync/atomic" +) + +// indexCellRegion represents a simple pair for defining an integer valued region. +type indexCellRegion struct { + start int + size int +} + +// shapeRegion is a mapping from shape id to the region of the edges array +// it's stored in. +type shapeRegion struct { + id int32 + region indexCellRegion +} + +// edgeAndIDChain is an extension of Edge with fields for the edge id, +// chain id, and offset. It's useful to bundle these together when decoding +// ShapeIndex cells because it allows us to avoid repetitive edge and chain +// lookups in many cases. +type edgeAndIDChain struct { + Edge // Embed the Edge type. + ID int32 // ID of the edge within its shape. + Chain int32 // ID of the chain the edge belongs to. + Offset int32 // Offset of the edge within the chain. +} + +// indexCellData provides methods for working with ShapeIndexCell data. For +// larger queries like validation, we often look up edges multiple times, +// and sometimes need to work with the edges themselves, their edge IDs, or +// their chain and offset. +// +// ShapeIndexCell and the ClippedShape APIs fundamentally work with edge IDs +// and can't be re-worked without significant effort and loss of efficiency. +// This class provides an alternative API than repeatedly querying through the +// shapes in the index. +// +// This is meant to support larger querying and validation operations such as +// ValidationQuery that have to proceed cell-by cell through an index. +// +// To use, simply call loadCell() to decode the contents of a cell. +// +// This type promises that the edges will be looked up once when loadCell() is +// called, and the edges, edgeIDs, chain, and chain offsets are loaded into a +// contiguous chunk of memory that we can serve requests from via slices. +// Since the chain and offset are computed anyways when looking up an edge via +// the shape.Edge() API, we simply cache those values so the cost is minimal. +// +// The memory layout looks like this: +// +// | 0D Shapes | 1D Shapes | 2D Shapes | Dimensions +// | 5 | 1 | 3 | 2 | 7 | 0 | 6 | 4 | 8 | Shapes +// [ ......................... Edges ..........................] Edges +// +// This allows us to look up individual shapes very quickly, as well as all +// shapes in a given dimension or contiguous range of dimensions: +// +// Edges() - Return slice over all edges. +// ShapeEdges() - Return slice over edges of a given shape. +// DimEdges() - Return slice over all edges of a given dimension. +// DimRangeEges() - Return slice over all edges of a range of dimensions. +// +// We use a stable sort, so similarly to ShapeIndexCell, we promise that +// shapes _within a dimension_ are in the same order they are in the index +// itself, and the edges _within a shape_ are similarly in the same order. +type indexCellData struct { + index *ShapeIndex + indexCell *ShapeIndexCell + cellID CellID + + // Computing the cell center and Cell can cost as much as looking up the + // edges themselves, so defer doing it until needed. + mu sync.Mutex + s2CellSet atomic.Bool + s2Cell Cell + centerSet atomic.Bool + cellCenter Point + + // Dimensions that we wish to decode, the default is all of them. + dimWanted [3]bool + + // Storage space for edges of the current cell. + edges []edgeAndIDChain + + // Mapping from shape id to the region of the edges array it's stored in. + shapeRegions []shapeRegion + + // Region for each dimension we might encounter. + dimRegions [3]indexCellRegion +} + +// newIndexCellData creates a new indexCellData with the expected defaults. +func newIndexCellData() *indexCellData { + return &indexCellData{ + dimWanted: [3]bool{true, true, true}, + } +} + +// newindexCellDataFromCell creates a new indexCellData and loads the given +// cell data. +func newIndexCellDataFromCell(index *ShapeIndex, id CellID, cell *ShapeIndexCell) *indexCellData { + d := newIndexCellData() + d.loadCell(index, id, cell) + return d +} + +// reset resets internal state to defaults. +// The next call to loadCell() will process the cell regardless of whether +// it was already loaded. Should also be called when processing a new index. +func (d *indexCellData) reset() { + d.index = nil + d.indexCell = nil + d.edges = d.edges[:0] + d.shapeRegions = d.shapeRegions[:0] + d.dimWanted = [3]bool{true, true, true} +} + +// setDimWanted configures whether a particular dimension of shape should be decoded. +func (d *indexCellData) setDimWanted(dim int, wanted bool) { + if dim < 0 || dim > 2 { + return + } + d.dimWanted[dim] = wanted +} + +// cell returns the S2 Cell for the current index cell, loading it if it +// was not already set. +func (d *indexCellData) cell() Cell { + if !d.s2CellSet.Load() { + d.mu.Lock() + if !d.s2CellSet.Load() { + d.s2Cell = CellFromCellID(d.cellID) + d.s2CellSet.Store(true) + } + d.mu.Unlock() + } + return d.s2Cell +} + +// center returns the center point of the current index cell, loading it +// if it was not already set. +func (d *indexCellData) center() Point { + if !d.centerSet.Load() { + d.mu.Lock() + if !d.centerSet.Load() { + d.cellCenter = d.cellID.Point() + d.centerSet.Store(true) + } + d.mu.Unlock() + } + return d.cellCenter +} + +// loadCell loads the data from the given cell, previous cell data is cleared. +// Both the index and cell lifetimes must span the lifetime until this +// indexCellData is destroyed or loadCell() is called again. +// +// If the index, id and cell pointer are the same as in the previous call to +// loadCell, loading is not performed since we already have the data decoded. +func (d *indexCellData) loadCell(index *ShapeIndex, id CellID, cell *ShapeIndexCell) { + // If this is still the same ShapeIndexCell as last time, then we are good. + if d.index == index && d.cellID == id { + return + } + + d.index = index + + // Cache cell information. + d.indexCell = cell + d.cellID = id + + // Reset atomic flags so we'll recompute cached values. + d.s2CellSet.Store(false) + d.centerSet.Store(false) + + // Clear previous edges + d.edges = d.edges[:0] + d.shapeRegions = d.shapeRegions[:0] + + // Reset per-dimension region information. + for i := range d.dimRegions { + d.dimRegions[i] = indexCellRegion{} + } + + minDim := 0 + for minDim <= 2 && !d.dimWanted[minDim] { + minDim++ + } + + maxDim := 2 + for maxDim >= 0 && !d.dimWanted[maxDim] { + maxDim-- + } + + // No dimensions wanted, we're done. + if minDim > 2 || maxDim < 0 { + return + } + + // For each dimension, load the edges from all shapes of that dimension + for dim := minDim; dim <= maxDim; dim++ { + dimStart := len(d.edges) + + for _, clipped := range cell.shapes { + shapeID := clipped.shapeID + shape := index.Shape(shapeID) + + // Only process the current dimension. + if shape.Dimension() != dim { + continue + } + + // In the event we wanted dimensions 0 and 2, but not 1. + // + // TODO(rsned): Should this be earlier in the loop? + // Why not skip this dim altogether if it's not wanted? + // Same question for C++. + if !d.dimWanted[dim] { + continue + } + + // Materialize clipped shape edges into the edges + // slice. Track where we start so we can add + // information about the region for this shape. + shapeStart := len(d.edges) + for _, edgeID := range clipped.edges { + // Looking up an edge requires looking up + // which chain it's in, which is often a binary + // search. So let's manually lookup the chain + // information and use that to find the edge, + // so we only have to do that search once. + position := shape.ChainPosition(edgeID) + edge := shape.ChainEdge(position.ChainID, position.Offset) + d.edges = append(d.edges, edgeAndIDChain{ + Edge: edge, + ID: int32(edgeID), + Chain: int32(position.ChainID), + Offset: int32(position.Offset), + }) + } + + // Note which block of edges belongs to the shape. + d.shapeRegions = append(d.shapeRegions, shapeRegion{ + id: shapeID, + region: indexCellRegion{ + start: shapeStart, + size: len(d.edges) - shapeStart, + }, + }) + } + + // Save region information for the current dimension. + d.dimRegions[dim] = indexCellRegion{ + start: dimStart, + size: len(d.edges) - dimStart, + } + } +} + +// shapeEdges returns a slice of the edges in the current cell for a given shape. +func (d *indexCellData) shapeEdges(shapeID int32) []edgeAndIDChain { + for _, sr := range d.shapeRegions { + if sr.id == shapeID { + region := sr.region + if region.start < len(d.edges) { + return d.edges[region.start : region.start+region.size] + } + return nil + } + } + return nil +} + +// dimEdges returns a slice of the edges in the current cell for the given +// dimension. +func (d *indexCellData) dimEdges(dim int) []edgeAndIDChain { + if dim < 0 || dim > 2 { + return nil + } + + region := d.dimRegions[dim] + if region.start < len(d.edges) { + return d.edges[region.start : region.start+region.size] + } + return nil +} + +// dimRangeEdges returns a slice of the edges in the current cell for a +// range of dimensions. +func (d *indexCellData) dimRangeEdges(dim0, dim1 int) []edgeAndIDChain { + if dim0 > dim1 || dim0 < 0 || dim0 > 2 || dim1 < 0 || dim1 > 2 { + return nil + } + + start := d.dimRegions[dim0].start + size := 0 + + for dim := dim0; dim <= dim1; dim++ { + start = minInt(start, d.dimRegions[dim].start) + size += d.dimRegions[dim].size + } + + if start < len(d.edges) { + return d.edges[start:size] + } + return nil +} + +// TODO(rsned): Differences from C++ +// shapeContains diff --git a/s2/index_cell_data_test.go b/s2/index_cell_data_test.go new file mode 100644 index 0000000..5e004b6 --- /dev/null +++ b/s2/index_cell_data_test.go @@ -0,0 +1,155 @@ +// Copyright 2025 The S2 Geometry Project Authors. All rights reserved. +// +// 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 s2 + +import ( + "testing" +) + +// TestIndexCellDataDimensionFilteringWorks verifies that we can filter +// shapes by dimension when loading cells. +func testIndexCellDataDimensionFilteringWorks(t *testing.T) { + tests := []struct { + dimWanted [3]bool + dimEmpty [3]bool + }{ + { + // Check that we get all dimensions by default. + dimWanted: [3]bool{true, true, true}, + dimEmpty: [3]bool{false, false, false}, + }, + { + // No dimensions should work too, we just don't decode edges. + dimWanted: [3]bool{false, false, false}, + dimEmpty: [3]bool{true, true, true}, + }, + // Should be able to get ranges even if a dimension is off. + { + dimWanted: [3]bool{false, true, true}, + dimEmpty: [3]bool{true, false, false}, + }, + { + dimWanted: [3]bool{false, true, false}, + dimEmpty: [3]bool{true, false, true}, + }, + { + dimWanted: [3]bool{true, false, true}, + dimEmpty: [3]bool{false, true, false}, + }, + { + dimWanted: [3]bool{true, false, false}, + dimEmpty: [3]bool{false, true, true}, + }, + { + dimWanted: [3]bool{false, false, true}, + dimEmpty: [3]bool{true, true, false}, + }, + } + + for i, test := range tests { + index := makeShapeIndex("0:0" + + "#1:1, 2:2" + + "#1:0, 0:1, -1:0, 0:-1") + + iter := index.Iterator() + iter.Begin() + data := newIndexCellData() + for j := 0; j < 3; j++ { + data.dimWanted[j] = test.dimWanted[j] + } + data.loadCell(index, iter.CellID(), iter.cell) + + for dim := 0; dim < 2; dim++ { + empty := len(data.dimEdges(dim)) == 0 + // Didn't get the answer we expected. + if empty != test.dimEmpty[dim] { + if test.dimEmpty[dim] { + t.Errorf("%d. Expected empty dimEdges(%d), but got non-empty", i, dim) + } else { + t.Errorf("%d. Expected non-empty dimEdges(%d), but got empty", i, dim) + } + } + } + } + + // Finally, test dimRangeEdges as well. + index := makeShapeIndex("0:0" + + "#1:1, 2:2" + + "#1:0, 0:1, -1:0, 0:-1") + + iter := index.Iterator() + iter.Begin() + data := newIndexCellData() + data.setDimWanted(0, false) + data.setDimWanted(1, true) + data.setDimWanted(2, true) + data.loadCell(index, iter.CellID(), iter.cell) + if len(data.dimRangeEdges(0, 0)) != 0 { + t.Error("Expected empty dimRangeEdges(0, 0)") + } + if len(data.dimRangeEdges(0, 2)) == 0 { + t.Error("Expected non-empty dimRangeEdges(0, 2)") + } +} + +// TestIndexCellDataCellAndCenterRecomputed verifies that when we load a +// new unique cell, the cached cell and center values are updated if we +// access them. +func testIndexCellDataCellAndCenterRecomputed(t *testing.T) { + // A line between two faces will guarantee we get at least two cells. + index := makeShapeIndex("# 0:0, 0:-90 #") + + // Get the first cell from the index + iter := index.Iterator() + iter.Begin() + + data := newIndexCellData() + data.loadCell(index, iter.CellID(), iter.cell) + + center0 := data.center() + cell0 := data.cell() + + // Move to the next cell + iter.Next() + if iter.Done() { + t.Errorf("Expected at least two cells in the index") + } + + // Load the second cell and check that the cached values change + data.loadCell(index, iter.CellID(), iter.cell) + center1 := data.center() + cell1 := data.cell() + + if cell0 == cell1 { + t.Error("Expected different cells") + } + if center0 == center1 { + t.Error("Expected different cell centers") + } + + // Load the same cell again, nothing should change + data.loadCell(index, iter.CellID(), iter.cell) + center2 := data.center() + cell2 := data.cell() + + if cell1 != cell2 { + t.Error("Expected same cells when reloading the same cell") + } + if center1 != center2 { + t.Error("Expected same cell centers when reloading the same cell") + } +} + +// TODO(rsned): Test shapeContains diff --git a/s2/lax_polygon.go b/s2/lax_polygon.go index ba9cf22..b6a25a4 100644 --- a/s2/lax_polygon.go +++ b/s2/lax_polygon.go @@ -14,12 +14,14 @@ package s2 +import "slices" + // Shape interface enforcement var _ Shape = (*LaxPolygon)(nil) // LaxPolygon represents a region defined by a collection of zero or more // closed loops. The interior is the region to the left of all loops. This -// is similar to Polygon except that this class supports polygons +// is similar to Polygon except that this type supports polygons // with degeneracies. Degeneracies are of two types: degenerate edges (from a // vertex to itself) and sibling edge pairs (consisting of two oppositely // oriented edges). Degeneracies can represent either "shells" or "holes" @@ -28,50 +30,50 @@ var _ Shape = (*LaxPolygon)(nil) // degenerate hole. Such edges form part of the boundary of the polygon. // // Loops with fewer than three vertices are interpreted as follows: -// - A loop with two vertices defines two edges (in opposite directions). -// - A loop with one vertex defines a single degenerate edge. -// - A loop with no vertices is interpreted as the "full loop" containing -// -// all points on the sphere. If this loop is present, then all other loops -// must form degeneracies (i.e., degenerate edges or sibling pairs). For -// example, two loops {} and {X} would be interpreted as the full polygon -// with a degenerate single-point hole at X. +// - A loop with two vertices defines two edges (in opposite directions). +// - A loop with one vertex defines a single degenerate edge. +// - A loop with no vertices is interpreted as the "full loop" containing +// all points on the sphere. If this loop is present, then all other loops +// must form degeneracies (i.e., degenerate edges or sibling pairs). For +// example, two loops {} and {X} would be interpreted as the full polygon +// with a degenerate single-point hole at X. // // LaxPolygon does not have any error checking, and it is perfectly fine to // create LaxPolygon objects that do not meet the requirements below (e.g., in // order to analyze or fix those problems). However, LaxPolygons must satisfy // some additional conditions in order to perform certain operations: // -// - In order to be valid for point containment tests, the polygon must -// -// satisfy the "interior is on the left" rule. This means that there must -// not be any crossing edges, and if there are duplicate edges then all but -// at most one of them must belong to a sibling pair (i.e., the number of -// edges in opposite directions must differ by at most one). -// -// - To be valid for polygon operations (BoundaryOperation), degenerate +// - In order to be valid for point containment tests, the polygon must +// satisfy the "interior is on the left" rule. This means that there must +// not be any crossing edges, and if there are duplicate edges then all but +// at most one of them must belong to a sibling pair (i.e., the number of +// edges in opposite directions must differ by at most one). // -// edges and sibling pairs cannot coincide with any other edges. For -// example, the following situations are not allowed: +// - To be valid for polygon operations (BooleanOperation), degenerate +// edges and sibling pairs cannot coincide with any other edges. For +// example, the following situations are not allowed: // -// {AA, AA} // degenerate edge coincides with another edge -// {AA, AB} // degenerate edge coincides with another edge -// {AB, BA, AB} // sibling pair coincides with another edge +// {AA, AA} // degenerate edge coincides with another edge +// {AA, AB} // degenerate edge coincides with another edge +// {AB, BA, AB} // sibling pair coincides with another edge // // Note that LaxPolygon is much faster to initialize and is more compact than // Polygon, but unlike Polygon it does not have any built-in operations. -// Instead you should use ShapeIndex based operations such as BoundaryOperation, +// Instead you should use ShapeIndex based operations such as BooleanOperation, // ClosestEdgeQuery, etc. type LaxPolygon struct { numLoops int vertices []Point - numVerts int - cumulativeVertices []int + numVerts int + // when numLoops > 1, store a list of size (numLoops+1) where "i" + // represents the total number of vertices in loops 0..i-1. + loopStarts []int +} - // TODO(roberts): C++ adds a prevLoop int field that claims to boost - // chain position lookups by 1.5-4.5x. Benchmark to see if this - // is useful here. +// LaxPolygonFromLoops creates a LaxPolygon from the given Polygon. +func LaxPolygonFromLoops(l []Loop) *LaxPolygon { + return nil } // LaxPolygonFromPolygon creates a LaxPolygon from the given Polygon. @@ -85,7 +87,17 @@ func LaxPolygonFromPolygon(p *Polygon) *LaxPolygon { copy(spans[i], loop.vertices) } } - return LaxPolygonFromPoints(spans) + lax := LaxPolygonFromPoints(spans) + + // Polygon and LaxPolygonShape holes are oriented oppositely, so we need + // to reverse the orientation of any loops representing holes. + for i := 0; i < p.NumLoops(); i++ { + if p.Loop(i).IsHole() { + v0 := lax.loopStarts[i] + slices.Reverse(lax.vertices[v0 : v0+lax.numLoopVertices(i)]) + } + } + return lax } // LaxPolygonFromPoints creates a LaxPolygon from the given points. @@ -99,19 +111,18 @@ func LaxPolygonFromPoints(loops [][]Point) *LaxPolygon { case 1: p.numVerts = len(loops[0]) p.vertices = make([]Point, p.numVerts) + p.loopStarts = []int{0, 0} copy(p.vertices, loops[0]) default: - p.cumulativeVertices = make([]int, p.numLoops+1) - numVertices := 0 + p.numVerts = 0 + p.loopStarts = make([]int, p.numLoops+1) for i, loop := range loops { - p.cumulativeVertices[i] = numVertices - numVertices += len(loop) + p.loopStarts[i] = p.numVerts + p.numVerts += len(loop) + p.vertices = append(p.vertices, loop...) } - p.cumulativeVertices[p.numLoops] = numVertices - for _, points := range loops { - p.vertices = append(p.vertices, points...) - } + p.loopStarts[p.numLoops] = p.numVerts } return p } @@ -121,7 +132,7 @@ func (p *LaxPolygon) numVertices() int { if p.numLoops <= 1 { return p.numVerts } - return p.cumulativeVertices[p.numLoops] + return p.loopStarts[p.numLoops] } // numLoopVertices reports the total number of vertices in the given loop. @@ -129,7 +140,7 @@ func (p *LaxPolygon) numLoopVertices(i int) int { if p.numLoops == 1 { return p.numVerts } - return p.cumulativeVertices[i+1] - p.cumulativeVertices[i] + return p.loopStarts[i+1] - p.loopStarts[i] } // loopVertex returns the vertex from loop i at index j. @@ -137,42 +148,20 @@ func (p *LaxPolygon) numLoopVertices(i int) int { // This requires: // // 0 <= i < len(loops) -// 0 <= j < len(loop[i].vertices) +// 0 <= j < numLoopVertices(i) func (p *LaxPolygon) loopVertex(i, j int) Point { if p.numLoops == 1 { return p.vertices[j] } - return p.vertices[p.cumulativeVertices[i]+j] + return p.vertices[p.loopStarts[i]+j] } func (p *LaxPolygon) NumEdges() int { return p.numVertices() } func (p *LaxPolygon) Edge(e int) Edge { - e1 := e + 1 - if p.numLoops == 1 { - // wrap the end vertex if this is the last edge. - if e1 == p.numVerts { - e1 = 0 - } - return Edge{p.vertices[e], p.vertices[e1]} - } - - // TODO(roberts): If this turns out to be performance critical in tests - // incorporate the maxLinearSearchLoops like in C++. - - // Check if e1 would cross a loop boundary in the set of all vertices. - nextLoop := 0 - for p.cumulativeVertices[nextLoop] <= e { - nextLoop++ - } - - // If so, wrap around to the first vertex of the loop. - if e1 == p.cumulativeVertices[nextLoop] { - e1 = p.cumulativeVertices[nextLoop-1] - } - - return Edge{p.vertices[e], p.vertices[e1]} + pos := p.ChainPosition(e) + return p.ChainEdge(pos.ChainID, pos.Offset) } func (p *LaxPolygon) Dimension() int { return 2 } @@ -184,10 +173,11 @@ func (p *LaxPolygon) ReferencePoint() ReferencePoint { return referencePointForS func (p *LaxPolygon) NumChains() int { return p.numLoops } func (p *LaxPolygon) Chain(i int) Chain { if p.numLoops == 1 { - return Chain{0, p.numVertices()} + return Chain{Start: 0, Length: p.numVertices()} } - start := p.cumulativeVertices[i] - return Chain{start, p.cumulativeVertices[i+1] - start} + + start := p.loopStarts[i] + return Chain{Start: start, Length: p.loopStarts[i+1] - start} } func (p *LaxPolygon) ChainEdge(i, j int) Edge { @@ -196,29 +186,33 @@ func (p *LaxPolygon) ChainEdge(i, j int) Edge { if j+1 != n { k = j + 1 } + if p.numLoops == 1 { - return Edge{p.vertices[j], p.vertices[k]} + return Edge{V0: p.vertices[j], V1: p.vertices[k]} } - base := p.cumulativeVertices[i] - return Edge{p.vertices[base+j], p.vertices[base+k]} + + start := p.loopStarts[i] + return Edge{V0: p.vertices[start+j], V1: p.vertices[start+k]} } -func (p *LaxPolygon) ChainPosition(e int) ChainPosition { +func (p *LaxPolygon) ChainPosition(edgeID int) ChainPosition { + if p.numLoops == 1 { - return ChainPosition{0, e} + return ChainPosition{ChainID: 0, Offset: edgeID} } - // TODO(roberts): If this turns out to be performance critical in tests - // incorporate the maxLinearSearchLoops like in C++. - - // Find the index of the first vertex of the loop following this one. - nextLoop := 1 - for p.cumulativeVertices[nextLoop] <= e { + // We need the loopStart that is less than or equal to the edgeID. + nextLoop := 0 + for p.loopStarts[nextLoop] <= edgeID { nextLoop++ } - return ChainPosition{p.cumulativeVertices[nextLoop] - p.cumulativeVertices[1], e - p.cumulativeVertices[nextLoop-1]} + return ChainPosition{ + ChainID: nextLoop - 1, + Offset: edgeID - p.loopStarts[nextLoop-1], + } } // TODO(roberts): Remaining to port from C++: -// EncodedLaxPolygon +// encode/decode +// Support for EncodedLaxPolygon diff --git a/s2/lax_polygon_test.go b/s2/lax_polygon_test.go index 70101cd..b0a0029 100644 --- a/s2/lax_polygon_test.go +++ b/s2/lax_polygon_test.go @@ -15,10 +15,13 @@ package s2 import ( + "math/rand" "testing" + + "github.com/golang/geo/s1" ) -func TestLaxPolygonShapeEmptyPolygon(t *testing.T) { +func TestLaxPolygonEmptyPolygon(t *testing.T) { shape := LaxPolygonFromPolygon((&Polygon{})) if got, want := shape.numLoops, 0; got != want { t.Errorf("shape.numLoops = %d, want %d", got, want) @@ -46,7 +49,7 @@ func TestLaxPolygonShapeEmptyPolygon(t *testing.T) { } } -func TestLaxPolygonFull(t *testing.T) { +func TestLaxPolygonFullPolygon(t *testing.T) { shape := LaxPolygonFromPolygon(PolygonFromLoops([]*Loop{makeLoop("full")})) if got, want := shape.numLoops, 1; got != want { t.Errorf("shape.numLoops = %d, want %d", got, want) @@ -124,7 +127,7 @@ func TestLaxPolygonSingleVertexPolygon(t *testing.T) { } } -func TestLaxPolygonShapeSingleLoopPolygon(t *testing.T) { +func TestLaxPolygonSingleLoopPolygon(t *testing.T) { vertices := parsePoints("0:0, 0:1, 1:1, 1:0") lenVerts := len(vertices) shape := LaxPolygonFromPolygon(PolygonFromLoops([]*Loop{LoopFromPoints(vertices)})) @@ -183,7 +186,7 @@ func TestLaxPolygonShapeSingleLoopPolygon(t *testing.T) { } } -func TestLaxPolygonShapeMultiLoopPolygon(t *testing.T) { +func TestLaxPolygonMultiLoopPolygon(t *testing.T) { // Test to make sure that the loops are oriented so that the interior of the // polygon is always on the left. loops := [][]Point{ @@ -245,7 +248,120 @@ func TestLaxPolygonShapeMultiLoopPolygon(t *testing.T) { } } -func TestLaxPolygonShapeDegenerateLoops(t *testing.T) { +// three ints is a tuple used in the many loop polygon test. +type threeInts struct { + e, i, j int +} + +func TestLaxPolygonManyLoopPolygon(t *testing.T) { + // Test a polygon with enough loops so that binary search is used to find + // the loop containing a given edge. + + const startingLoops = 100 + + loops := make([][]Point, startingLoops) + for i := 0; i < startingLoops; i++ { + center := PointFromLatLng(LatLngFromDegrees(0, float64(i))) + loops[i] = RegularLoop(center, s1.Angle(0.1)*s1.Degree, + randomUniformInt(3)).vertices + } + + shape := LaxPolygonFromPoints(loops) + + numLoops := len(loops) + + if shape.numLoops != numLoops { + t.Errorf("LaxPolygon num loops = %d, want %d", shape.numLoops, numLoops) + } + if shape.NumChains() != numLoops { + t.Errorf("LaxPolygon.NumChains() = %d, want %d", shape.NumChains(), numLoops) + } + + numVertices := 0 + for i := 0; i < numLoops; i++ { + loopLenI := len(loops[i]) + if loopLenI != shape.numLoopVertices(i) { + t.Errorf("loop[%d] num vertices = %d, want %d", i, shape.numLoopVertices(i), loopLenI) + } + if numVertices != shape.Chain(i).Start { + t.Errorf("LaxPolygon.Chain(%d).Start = %d, want %d", + i, shape.Chain(i).Start, numVertices) + } + if loopLenI != shape.Chain(i).Length { + t.Errorf("LaxPolygon.Chain(%d).Length = %d, want %d", + i, shape.Chain(i).Length, loopLenI) + } + for j := 0; j < loopLenI; j++ { + if loops[i][j] != shape.loopVertex(i, j) { + t.Errorf("loopVertex(%d, %d) = %v != original vertex %v", + i, j, shape.loopVertex(i, j), loops[i][j]) + } + e := numVertices + j + if shape.ChainPosition(e) != (ChainPosition{i, j}) { + t.Errorf("LaxPolygon.ChainPosition(%d) = %v,. want %v", + e, shape.ChainPosition(e), (ChainPosition{i, j})) + } + if loops[i][j] != shape.Edge(e).V0 { + t.Errorf("LaxPolygon.Edge(%d).V0 = %v, want %v", + e, shape.Edge(e).V0, loops[i][j]) + } + idx := (j + 1) % loopLenI + if loops[i][idx] != shape.Edge(e).V1 { + t.Errorf("LaxPolygon.Edge(%d).V1 = %v, want %v", + e, shape.Edge(e).V1, loops[i][idx]) + } + } + numVertices += loopLenI + } + if numVertices != shape.numVertices() { + t.Errorf("LaxPolygon.numVertices() = %d, want %d", shape.numVertices(), numVertices) + } + if numVertices != shape.NumEdges() { + t.Errorf("LaxPolygon.NumEdges() = %d, want %d", shape.NumEdges(), numVertices) + } + + // Now test all the edges in a random order in order. + edges := []threeInts{} + for i, e := 0, 0; i < numLoops; i++ { + for j := 0; j < len(loops[i]); j++ { + edges = append(edges, (threeInts{e, i, j})) + e++ + } + } + + // C++ uses the Mersienne Twister to shuffle the elements. For now just + // use rand.Shuffle unless it proves problematic. + rand.Shuffle(numVertices, func(i, j int) { + edges[i], edges[j] = edges[j], edges[i] + }) + + for _, edge := range edges { + if shape.ChainPosition(edge.e) != (ChainPosition{edge.i, edge.j}) { + t.Errorf("addasdaa") + } + v0 := loops[edge.i][edge.j] + v1 := loops[edge.i][(edge.j+1)%len(loops[edge.i])] + if shape.Edge(edge.e) != (Edge{v0, v1}) { + t.Errorf("sfsdaa") + } + } +} + +func TestLaxPolygonMultiLoopS2Polygon(t *testing.T) { + // Verify that the orientation of loops representing holes is reversed when + // converting from a Polygon to a LaxPolygonShape. + polygon := makePolygon("0:0, 0:3, 3:3; 1:1, 1:2, 2:2", true) + shape := LaxPolygonFromPolygon(polygon) + for i, loop := range polygon.Loops() { + for j := 0; j < loop.NumVertices(); j++ { + if loop.OrientedVertex(j) != shape.loopVertex(i, j) { + t.Errorf("LaxPolygon vertex %d in loop %d should equal the original loops oriented vertex but does not", j, i) + } + } + } +} + +func TestLaxPolygonDegenerateLoops(t *testing.T) { loops := [][]Point{ parsePoints("1:1, 1:2, 2:2, 1:2, 1:3, 1:2, 1:1"), parsePoints("0:0, 0:3, 0:6, 0:9, 0:6, 0:3, 0:0"), @@ -258,7 +374,7 @@ func TestLaxPolygonShapeDegenerateLoops(t *testing.T) { } } -func TestLaxPolygonShapeInvertedLoops(t *testing.T) { +func TestLaxPolygonInvertedLoops(t *testing.T) { loops := [][]Point{ parsePoints("1:2, 1:1, 2:2"), parsePoints("3:4, 3:3, 4:4"), @@ -271,10 +387,9 @@ func TestLaxPolygonShapeInvertedLoops(t *testing.T) { } // TODO(roberts): Remaining to port: -// LaxPolygonManyLoopPolygon -// LaxPolygonMultiLoopS2Polygon // LaxPolygonCompareToLoop once fractal testing is added. // LaxPolygonCoderWorks // LaxPolygonChainIteratorWorks // LaxPolygonChainVertexIteratorWorks // +// Add testLaxPolygonEncoding to all the above tests as well.