diff --git a/s2/builder.go b/s2/builder.go index 004dd5b..3331863 100644 --- a/s2/builder.go +++ b/s2/builder.go @@ -14,6 +14,12 @@ package s2 +import ( + "math" + + "github.com/golang/geo/s1" +) + const ( // maxEdgeDeviationRatio is set so that MaxEdgeDeviation will be large enough // compared to snapRadius such that edge splitting is rare. @@ -33,3 +39,345 @@ const ( // when MaxEdgeDeviation is exceeded. maxEdgeDeviationRatio = 1.1 ) + +// isFullPolygonPredicate is an interface for determining if Polygons are +// full or not. For output layers that represent polygons, there is an ambiguity +// inherent in spherical geometry that does not exist in planar geometry. +// Namely, if a polygon has no edges, does it represent the empty polygon +// (containing no points) or the full polygon (containing all points)? This +// ambiguity also occurs for polygons that consist only of degeneracies, e.g. +// a degenerate loop with only two edges could be either a degenerate shell in +// the empty polygon or a degenerate hole in the full polygon. +// +// To resolve this ambiguity, an IsFullPolygonPredicate may be specified for +// each output layer (see AddIsFullPolygonPredicate below). If the output +// after snapping consists only of degenerate edges and/or sibling pairs +// (including the case where there are no edges at all), then the layer +// implementation calls the given predicate to determine whether the polygon +// is empty or full except for those degeneracies. The predicate is given +// an S2Builder::Graph containing the output edges, but note that in general +// the predicate must also have knowledge of the input geometry in order to +// determine the correct result. +// +// This predicate is only needed by layers that are assembled into polygons. +// It is not used by other layer types. +type isFullPolygonPredicate func(g *graph) (bool, error) + +// builder is a tool for assembling polygonal geometry from edges. Here are +// some of the things it is designed for: +// +// 1. Building polygons, polylines, and polygon meshes from unsorted +// collections of edges. +// +// 2. Snapping geometry to discrete representations (such as CellID centers +// or E7 lat/lng coordinates) while preserving the input topology and with +// guaranteed error bounds. +// +// 3. Simplifying geometry (e.g. for indexing, display, or storage). +// +// 4. Importing geometry from other formats, including repairing geometry +// that has errors. +// +// 5. As a tool for implementing more complex operations such as polygon +// intersections and unions. +// +// The implementation is based on the framework of "snap rounding". Unlike +// most snap rounding implementations, Builder defines edges as geodesics on +// the sphere (straight lines) and uses the topology of the sphere (i.e., +// there are no "seams" at the poles or 180th meridian). The algorithm is +// designed to be 100% robust for arbitrary input geometry. It offers the +// following properties: +// +// - Guaranteed bounds on how far input vertices and edges can move during +// the snapping process (i.e., at most the given "snapRadius"). +// +// - Guaranteed minimum separation between edges and vertices other than +// their endpoints (similar to the goals of Iterated Snap Rounding). In +// other words, edges that do not intersect in the output are guaranteed +// to have a minimum separation between them. +// +// - Idempotency (similar to the goals of Stable Snap Rounding), i.e. if the +// input already meets the output criteria then it will not be modified. +// +// - Preservation of the input topology (up to the creation of +// degeneracies). This means that there exists a continuous deformation +// from the input to the output such that no vertex crosses an edge. In +// other words, self-intersections won't be created, loops won't change +// orientation, etc. +// +// - The ability to snap to arbitrary discrete point sets (such as CellID +// centers, E7 lat/lng points on the sphere, or simply a subset of the +// input vertices), rather than being limited to an integer grid. +// +// Here are some of its other features: +// +// - It can handle both directed and undirected edges. Undirected edges can +// be useful for importing data from other formats, e.g. where loops have +// unspecified orientations. +// +// - It can eliminate self-intersections by finding all edge pairs that cross +// and adding a new vertex at each intersection point. +// +// - It can simplify polygons to within a specified tolerance. For example, +// if two vertices are close enough they will be merged, and if an edge +// passes nearby a vertex then it will be rerouted through that vertex. +// Optionally, it can also detect nearly straight chains of short edges and +// replace them with a single long edge, while maintaining the same +// accuracy, separation, and topology guarantees ("simplifyEdgeChains"). +// +// - It supports many different output types through the concept of "layers" +// (polylines, polygons, polygon meshes, etc). You can build multiple +// layers at once in order to ensure that snapping does not create +// intersections between different objects (for example, you can simplify a +// set of contour lines without the risk of having them cross each other). +// +// - It supports edge labels, which allow you to attach arbitrary information +// to edges and have it preserved during the snapping process. (This can +// also be achieved using layers, at a coarser level of granularity.) +// +// Caveats: +// +// - Because Builder only works with edges, it cannot distinguish between +// the empty and full polygons. If your application can generate both the +// empty and full polygons, you must implement logic outside of this type. +// +// Example showing how to snap a polygon to E7 coordinates: +// +// builder := NewBuilder(newBuilderOptions(IntLatLngSnapFunction(7))); +// var output *Polygon +// builder.StartLayer(NewPolygonLayer(output)) +// builder.AddPolygon(input); +// if err := builder.Build(); err != nil { +// fmt.Printf("error building: %v\n"), err +// ... +// } +// +// TODO(rsned): Make the type public when Builder is ready. +type builder struct { + opts *builderOptions + + // The maximum distance (inclusive) that a vertex can move when snapped, + // equal to options.SnapFunction().SnapRadius()). + siteSnapRadiusCA s1.ChordAngle + + // The maximum distance (inclusive) that an edge can move when snapping to a + // snap site. It can be slightly larger than the site snap radius when + // edges are being split at crossings. + edgeSnapRadiusCA s1.ChordAngle + + // True if we need to check that snapping has not changed the input topology + // around any vertex (i.e. Voronoi site). Normally this is only necessary for + // forced vertices, but if the snap radius is very small (e.g., zero) and + // splitCrossingedges is true then we need to do this for all vertices. + // In all other situations, any snapped edge that crosses a vertex will also + // be closer than minEdgeVertexSeparation() to that vertex, which will + // cause us to add a separation site anyway. + checkAllSiteCrossings bool + + maxEdgeDeviation s1.Angle + edgeSiteQueryRadiusCA s1.ChordAngle + minEdgeLengthToSplitCA s1.ChordAngle + + minSiteSeparation s1.Angle + minSiteSeparationCA s1.ChordAngle + minEdgeSiteSeparationCA s1.ChordAngle + minEdgeSiteSeparationCALimit s1.ChordAngle + + maxAdjacentSiteSeparationCA s1.ChordAngle + + // The squared sine of the edge snap radius. This is equivalent to the snap + // radius (squared) for distances measured through the interior of the + // sphere to the plane containing an edge. This value is used only when + // interpolating new points along edges (see GetSeparationSite). + edgeSnapRadiusSin2 float64 + + // True if snapping was requested. This is true if either snapRadius() is + // positive, or splitCrossingEdges() is true (which implicitly requests + // snapping to ensure that both crossing edges are snapped to the + // intersection point). + snappingRequested bool + + // Initially false, and set to true when it is discovered that at least one + // input vertex or edge does not meet the output guarantees (e.g., that + // vertices are separated by at least snapFunction.minVertexSeparation). + snappingNeeded bool + + // A flag indicating whether labelSet has been modified since the last + // time labelSetID was computed. + labelSetModified bool + + inputVertices []Point + // inputEdges []builderInputEdge + + layers []*builderLayer + layerOptions []*graphOptions + layerBegins []int32 + layerIsFullPolygonPredicates []isFullPolygonPredicate + + // Each input edge has "label set id" (an int32) representing the set of + // labels attached to that edge. This vector is populated only if at least + // one label is used. + labelSetIDs []int32 + labelSetLexicon *idSetLexicon + + // The current set of labels (represented as a stack). + labelSet []int32 + + // The labelSetID corresponding to the current label set, computed on demand + // (by adding it to labelSetLexicon()). + labelSetID int32 + + // The remaining fields are used for snapping and simplifying. + + // The number of sites specified using forceVertex(). These sites are + // always at the beginning of the sites vector. + numForcedSites int32 + + // The set of snapped vertex locations ("sites"). + sites []Point + + // A map from each input edge to the set of sites "nearby" that edge, + // defined as the set of sites that are candidates for snapping and/or + // avoidance. Note that compactarray will inline up to two sites, which + // usually takes care of the vast majority of edges. Sites are kept sorted + // by increasing distance from the origin of the input edge. + // + // Once snapping is finished, this field is discarded unless edge chain + // simplification was requested, in which case instead the sites are + // filtered by removing the ones that each edge was snapped to, leaving only + // the "sites to avoid" (needed for simplification). + edgeSites [][]int32 + + // TODO(rsned): Add memoryTracker if it becomes available. +} + +// init initializes this instance with the given options. +func (b *builder) init(opts *builderOptions) { + b.opts = opts + + snapFunc := opts.snapFunction + sr := snapFunc.SnapRadius() + + // Cap the snap radius to the limit. + sr = min(sr, maxSnapRadius) + + // Convert the snap radius to an ChordAngle. This is the "true snap + // radius" used when evaluating exact predicates. + b.siteSnapRadiusCA = s1.ChordAngleFromAngle(sr) + + // When intersectionTolerance is non-zero we need to use a larger snap + // radius for edges than for vertices to ensure that both edges are snapped + // to the edge intersection location. This is because the computed + // intersection point is not exact; it may be up to intersectionTolerance + // away from its true position. The computed intersection point might then + // be snapped to some other vertex up to SnapRadius away. So to ensure + // that both edges are snapped to a common vertex, we need to increase the + // snap radius for edges to at least the sum of these two values (calculated + // conservatively). + edgeSnapRadius := opts.edgeSnapRadius() + b.edgeSnapRadiusCA = roundUp(edgeSnapRadius) + b.snappingRequested = (edgeSnapRadius > 0) + + // Compute the maximum distance that a vertex can be separated from an + // edge while still affecting how that edge is snapped. + b.maxEdgeDeviation = opts.maxEdgeDeviation() + b.edgeSiteQueryRadiusCA = s1.ChordAngleFromAngle(b.maxEdgeDeviation + + snapFunc.MinEdgeVertexSeparation()) + + // Compute the maximum edge length such that even if both endpoints move by + // the maximum distance allowed (i.e., edgeSnapRadius), the center of the + // edge will still move by less than maxEdgeDeviation. This saves us a + // lot of work since then we don't need to check the actual deviation. + if !b.snappingRequested { + b.minEdgeLengthToSplitCA = s1.InfChordAngle() + } else { + // This value varies between 30 and 50 degrees depending on + // the snap radius. + b.minEdgeLengthToSplitCA = s1.ChordAngleFromAngle(s1.Angle(2 * + math.Acos(math.Sin(edgeSnapRadius.Radians())/ + math.Sin(b.maxEdgeDeviation.Radians())))) + } + + // In rare cases we may need to explicitly check that the input topology is + // preserved, i.e. that edges do not cross vertices when snapped. This is + // only necessary (1) for vertices added using forceVertex, and (2) when the + // snap radius is smaller than intersectionTolerance (which is typically + // either zero or intersectionError, about 9e-16 radians). This + // condition arises because when a geodesic edge is snapped, the edge center + // can move further than its endpoints. This can cause an edge to pass on the + // wrong side of an input vertex. (Note that this could not happen in a + // planar version of this algorithm.) Usually we don't need to consider this + // possibility explicitly, because if the snapped edge passes on the wrong + // side of a vertex then it is also closer than minEdgeVertexSeparation + // to that vertex, which will cause a separation site to be added. + // + // If the condition below is true then we need to check all sites (i.e., + // snapped input vertices) for topology changes. However this is almost never + // the case because + // + // maxEdgeDeviation() == 1.1 * edgeSnapRadius + // and minEdgeVertexSeparation() >= 0.219 * SnapRadius + // + // for all currently implemented snap functions. The condition below is + // only true when intersectionTolerance() is non-zero (which causes + // edgeSnapRadius() to exceed SnapRadius() by intersectionError) and + // SnapRadius() is very small (at most intersectionError / 1.19). + b.checkAllSiteCrossings = (opts.maxEdgeDeviation() > + opts.edgeSnapRadius()+snapFunc.MinEdgeVertexSeparation()) + + // TODO(rsned): need to add check that b.checkAllSiteCrossings is false when tolerance is <= 0. + // if opts.intersectionTolerance <= 0 { + // } + + // To implement idempotency, we check whether the input geometry could + // possibly be the output of a previous Builder invocation. This involves + // testing whether any site/site or edge/site pairs are too close together. + // This is done using exact predicates, which require converting the minimum + // separation values to a ChordAngle. + b.minSiteSeparation = snapFunc.MinVertexSeparation() + b.minSiteSeparationCA = s1.ChordAngleFromAngle(b.minSiteSeparation) + b.minEdgeSiteSeparationCA = s1.ChordAngleFromAngle(snapFunc.MinEdgeVertexSeparation()) + + // This is an upper bound on the distance computed by ClosestPointQuery + // where the true distance might be less than minEdgeSiteSeparationCA. + b.minEdgeSiteSeparationCALimit = addPointToEdgeError(b.minEdgeSiteSeparationCA) + + // Compute the maximum possible distance between two sites whose Voronoi + // regions touch. (The maximum radius of each Voronoi region is + // edgeSnapRadius.) Then increase this bound to account for errors. + b.maxAdjacentSiteSeparationCA = addPointToPointError(roundUp(2 * opts.edgeSnapRadius())) + + // Finally, we also precompute sin^2(edgeSnapRadius), which is simply the + // squared distance between a vertex and an edge measured perpendicular to + // the plane containing the edge, and increase this value by the maximum + // error in the calculation to compare this distance against the bound. + d := math.Sin(opts.edgeSnapRadius().Radians()) + b.edgeSnapRadiusSin2 = d * d + b.edgeSnapRadiusSin2 += ((9.5*d+2.5+2*sqrt3)*d + 9*dblEpsilon) * dblEpsilon + + // Initialize the current label set. + b.labelSetID = emptySetID + b.labelSetModified = false + + // If snapping was requested, we try to determine whether the input geometry + // already meets the output requirements. This is necessary for + // idempotency, and can also save work. If we discover any reason that the + // input geometry needs to be modified, snappingNeeded is set to true. + b.snappingNeeded = false + + // TODO(rsned): Memory tracker init. +} + +// roundUp rounds the given angle up by the max error and returns it as a chord angle. +func roundUp(a s1.Angle) s1.ChordAngle { + ca := s1.ChordAngleFromAngle(a) + return ca.Expanded(ca.MaxAngleError()) +} + +func addPointToPointError(ca s1.ChordAngle) s1.ChordAngle { + return ca.Expanded(ca.MaxPointError()) +} + +func addPointToEdgeError(ca s1.ChordAngle) s1.ChordAngle { + return ca.Expanded(minUpdateDistanceMaxError(ca)) +} diff --git a/s2/builder_graph.go b/s2/builder_graph.go new file mode 100644 index 0000000..149d35c --- /dev/null +++ b/s2/builder_graph.go @@ -0,0 +1,40 @@ +// 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 + +// A Graph represents a collection of snapped edges that is passed +// to a Layer for assembly. (Example layers include polygons, polylines, and +// polygon meshes.) The Graph object does not own any of its underlying data; +// it is simply a view of data that is stored elsewhere. You will only +// need this interface if you want to implement a new Layer subtype. +// +// The graph consists of vertices and directed edges. Vertices are numbered +// sequentially starting from zero. An edge is represented as a pair of +// vertex ids. The edges are sorted in lexicographic order, therefore all of +// the outgoing edges from a particular vertex form a contiguous range. +// +// TODO(rsned): Consider pulling out the methods that are helper functions for +// Layer implementations (such as getDirectedLoops) into a builder_graph_util.go. +type graph struct { + opts *graphOptions + numVertices int + vertices []Point + edges []graphEdge + inputEdgeIDSetIDs []int32 + inputEdgeIDSetLexicon *idSetLexicon + labelSetIDs []int32 + labelSetLexicon *idSetLexicon + isFullPolygonPredicate isFullPolygonPredicate +} diff --git a/s2/builder_graph_edge_processor.go b/s2/builder_graph_edge_processor.go new file mode 100644 index 0000000..30bbb19 --- /dev/null +++ b/s2/builder_graph_edge_processor.go @@ -0,0 +1,333 @@ +// 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 ( + "errors" + "sort" +) + +// maxVertexID is the maximum possible vertex ID, used as a sentinel value. +const maxVertexID = int32(^uint32(0) >> 1) + +// graphEdge is a tuple of edge IDs. +type graphEdge struct { + first, second int32 +} + +// reverse returns a new graphEdge with the vertices in reverse order. +func (g graphEdge) reverse() graphEdge { + return graphEdge{first: g.second, second: g.first} +} + +// minGraphEdge returns the minimum of two edges in lexicographic order. +func minGraphEdge(a, b graphEdge) graphEdge { + if a.first < b.first || (a.first == b.first && a.second <= b.second) { + return a + } + return b +} + +// stableLessThan compares two graphEdges for stable sorting. +// It uses the graphEdge IDs as a tiebreaker to ensure a stable sort. +func stableLessThan(a, b graphEdge, aID, bID int32) bool { + if a.first != b.first { + return a.first < b.first + } + if a.second != b.second { + return a.second < b.second + } + return aID < bID +} + +// edgeProcessor processes edges in a Graph to handle duplicates, siblings, +// and degenerate edges according to the specified GraphOptions. +type edgeProcessor struct { + options *graphOptions + edges []graphEdge + inputIDs []int32 + idSetLexicon *idSetLexicon + outEdges []int32 + inEdges []int32 + newEdges []graphEdge + newInputIDs []int32 +} + +// newedgeProcessor creates a new edgeProcessor with the given options and data. +func newEdgeProcessor(opts *graphOptions, edges []graphEdge, inputIDs []int32, + idSetLexicon *idSetLexicon) *edgeProcessor { + // opts should not be nil at this point, but just in case. + if opts == nil { + opts = defaultGraphOptions() + } + ep := &edgeProcessor{ + options: opts, + edges: edges, + inputIDs: inputIDs, + idSetLexicon: idSetLexicon, + inEdges: make([]int32, len(edges)), + outEdges: make([]int32, len(edges)), + } + + // Sort the outgoing and incoming edges in lexicographic order. + // We use a stable sort to ensure that each undirected edge becomes a sibling pair, + // even if there are multiple identical input edges. + + // Fill the slice with a number sequence. + for i := range ep.outEdges { + ep.outEdges[i] = int32(i) + } + stableSortEdges(ep.outEdges, func(a, b int32) bool { + return stableLessThan(ep.edges[a], ep.edges[b], a, b) + }) + + // Fill the slice with a number sequence. + for i := range ep.inEdges { + ep.inEdges[i] = int32(i) + } + stableSortEdges(ep.inEdges, func(a, b int32) bool { + return stableLessThan(ep.edges[a].reverse(), ep.edges[b].reverse(), a, b) + }) + + ep.newEdges = make([]graphEdge, 0, len(edges)) + ep.newInputIDs = make([]int32, 0, len(edges)) + + return ep +} + +// stableSortEdges performs a stable sort on the given slice of EdgeIDs using the provided less function. +func stableSortEdges(edges []int32, less func(a, b int32) bool) { + sort.SliceStable(edges, func(i, j int) bool { + return less(edges[i], edges[j]) + }) +} + +// addEdge adds a single edge with its input edge ID set to the new edges. +func (ep *edgeProcessor) addEdge(edge graphEdge, inputEdgeIDSetID int32) { + ep.newEdges = append(ep.newEdges, edge) + ep.newInputIDs = append(ep.newInputIDs, inputEdgeIDSetID) +} + +// addEdges adds multiple copies of the same edge with the same input edge ID set. +func (ep *edgeProcessor) addEdges(numEdges int, edge graphEdge, inputEdgeIDSetID int32) { + for i := 0; i < numEdges; i++ { + ep.addEdge(edge, inputEdgeIDSetID) + } +} + +// copyEdges copies a range of edges from the input edges to the new edges. +func (ep *edgeProcessor) copyEdges(outBegin, outEnd int) { + for i := outBegin; i < outEnd; i++ { + ep.addEdge(ep.edges[ep.outEdges[i]], ep.inputIDs[ep.outEdges[i]]) + } +} + +// mergeInputIDs merges the input edge ID sets for a range of edges. +func (ep *edgeProcessor) mergeInputIDs(outBegin, outEnd int) int32 { + if outEnd-outBegin == 1 { + return ep.inputIDs[ep.outEdges[outBegin]] + } + + var tmpIDs []int32 + for i := outBegin; i < outEnd; i++ { + tmpIDs = append(tmpIDs, ep.idSetLexicon.idSet(ep.inputIDs[ep.outEdges[i]])...) + } + return ep.idSetLexicon.add(tmpIDs...) +} + +// Run processes the edges according to the specified options. +func (ep *edgeProcessor) Run() error { + numEdges := len(ep.edges) + if numEdges == 0 { + return nil + } + + // Walk through the two sorted arrays performing a merge join. For each + // edge, gather all the duplicate copies of the edge in both directions + // (outgoing and incoming). Then decide what to do based on options and + // how many copies of the edge there are in each direction. + out, in := 0, 0 + outEdge := ep.edges[ep.outEdges[out]] + inEdge := ep.edges[ep.inEdges[in]] + sentinel := graphEdge{first: maxVertexID, second: maxVertexID} + + for { + edge := minGraphEdge(outEdge, inEdge.reverse()) + if edge == sentinel { + break + } + + outBegin := out + inBegin := in + for outEdge == edge { + out++ + if out == numEdges { + outEdge = sentinel + } else { + outEdge = ep.edges[ep.outEdges[out]] + } + } + for inEdge.reverse() == edge { + in++ + if in == numEdges { + inEdge = sentinel + } else { + inEdge = ep.edges[ep.inEdges[in]] + } + } + nOut := out - outBegin + nIn := in - inBegin + + if edge.first == edge.second { + // This is a degenerate edge. + if err := ep.handleDegenerateEdge(edge, outBegin, out, nOut, nIn, inBegin, in); err != nil { + return err + } + } else if err := ep.handleNormalEdge(edge, outBegin, out, nOut, nIn); err != nil { + return err + } + } + + // Replace the old edges with the new ones. + ep.edges = ep.newEdges + ep.inputIDs = ep.newInputIDs + return nil +} + +// handleDegenerateEdge handles a degenerate edge (an edge from a vertex to itself). +func (ep *edgeProcessor) handleDegenerateEdge(edge graphEdge, outBegin, outEnd int, nOut, nIn, inBegin, in int) error { + // This is a degenerate edge. + if nOut != nIn { + return errors.New("inconsistent number of degenerate edges") + } + if ep.options.degenerateEdges == degenerateEdgesDiscard { + return nil + } + if ep.options.degenerateEdges == degenerateEdgesDiscardExcess { + // Check if there are any non-degenerate incident edges. + if (outBegin > 0 && ep.edges[ep.outEdges[outBegin-1]].first == edge.first) || + (outEnd < len(ep.edges) && ep.edges[ep.outEdges[outEnd]].first == edge.first) || + (inBegin > 0 && ep.edges[ep.inEdges[inBegin-1]].second == edge.first) || + (in < len(ep.edges) && ep.edges[ep.inEdges[in]].second == edge.first) { + return nil // There were non-degenerate incident edges, so discard. + } + } + + // degenerateEdgesDiscardExcess also merges degenerate edges. + merge := ep.options.duplicateEdges == duplicateEdgesMerge || + ep.options.degenerateEdges == degenerateEdgesDiscardExcess + + if ep.options.edgeType == edgeTypeUndirected && + (ep.options.siblingPairs == siblingPairsRequire || + ep.options.siblingPairs == siblingPairsCreate) { + // When we have undirected edges and are guaranteed to have siblings, + // we cut the number of edges in half (see Builder). + if nOut&1 != 0 { + return errors.New("odd number of undirected degenerate edges") + } + if merge { + ep.addEdges(1, edge, ep.mergeInputIDs(outBegin, outEnd)) + } else { + ep.addEdges(nOut/2, edge, ep.mergeInputIDs(outBegin, outEnd)) + } + } else if merge { + if ep.options.edgeType == edgeTypeUndirected { + ep.addEdges(2, edge, ep.mergeInputIDs(outBegin, outEnd)) + } else { + ep.addEdges(1, edge, ep.mergeInputIDs(outBegin, outEnd)) + } + } else if ep.options.siblingPairs == siblingPairsDiscard || + ep.options.siblingPairs == siblingPairsDiscardExcess { + // Any SiblingPair option that discards edges causes the labels of all + // duplicate edges to be merged together (see Builder). + ep.addEdges(nOut, edge, ep.mergeInputIDs(outBegin, outEnd)) + } else { + ep.copyEdges(outBegin, outEnd) + } + return nil +} + +// handleNormalEdge handles a non-degenerate edge. +func (ep *edgeProcessor) handleNormalEdge(edge graphEdge, outBegin, outEnd int, nOut, nIn int) error { + switch ep.options.siblingPairs { + case siblingPairsKeep: + if nOut > 1 && ep.options.duplicateEdges == duplicateEdgesMerge { + ep.addEdge(edge, ep.mergeInputIDs(outBegin, outEnd)) + } else { + ep.copyEdges(outBegin, outEnd) + } + case siblingPairsDiscard: + if ep.options.edgeType == edgeTypeDirected { + // If nOut == nIn: balanced sibling pairs + // If nOut < nIn: unbalanced siblings, in the form AB, BA, BA + // If nOut > nIn: unbalanced siblings, in the form AB, AB, BA + if nOut <= nIn { + return nil + } + // Any option that discards edges causes the labels of all duplicate + // edges to be merged together (see Builder). + if ep.options.duplicateEdges == duplicateEdgesMerge { + ep.addEdges(1, edge, ep.mergeInputIDs(outBegin, outEnd)) + } else { + ep.addEdges(nOut-nIn, edge, ep.mergeInputIDs(outBegin, outEnd)) + } + } else { + if nOut&1 == 0 { + return nil + } + ep.addEdge(edge, ep.mergeInputIDs(outBegin, outEnd)) + } + case siblingPairsDiscardExcess: + if ep.options.edgeType == edgeTypeDirected { + // See comments above. The only difference is that if there are + // balanced sibling pairs, we want to keep one such pair. + if nOut < nIn { + return nil + } + if ep.options.duplicateEdges == duplicateEdgesMerge { + ep.addEdges(1, edge, ep.mergeInputIDs(outBegin, outEnd)) + } else { + ep.addEdges(maxInt(1, nOut-nIn), edge, ep.mergeInputIDs(outBegin, outEnd)) + } + } else { + ep.addEdges((nOut&1)+1, edge, ep.mergeInputIDs(outBegin, outEnd)) + } + case siblingPairsCreate, siblingPairsRequire: + // In C++, this check also checked the state of the S2Error passed in + // to make sure no previous errors had occured before now. + if ep.options.siblingPairs == siblingPairsRequire && + (ep.options.edgeType == edgeTypeDirected && nOut != nIn || + ep.options.edgeType == edgeTypeUndirected && nOut&1 != 0) { + return errors.New("expected all input edges to have siblings but some were missing") + } + + if ep.options.duplicateEdges == duplicateEdgesMerge { + ep.addEdge(edge, ep.mergeInputIDs(outBegin, outEnd)) + } else if ep.options.edgeType == edgeTypeUndirected { + // Convert graph to use directed edges instead (see documentation of + // siblingPairsCreate/siblingPairsRequire for undirected edges). + ep.addEdges((nOut+1)/2, edge, ep.mergeInputIDs(outBegin, outEnd)) + } else { + ep.copyEdges(outBegin, outEnd) + if nIn > nOut { + // Automatically created edges have no input edge ids or labels. + ep.addEdges(nIn-nOut, edge, emptySetID) + } + } + default: + return errors.New("invalid sibling pairs option") + } + return nil +} diff --git a/s2/builder_graph_edge_processor_test.go b/s2/builder_graph_edge_processor_test.go new file mode 100644 index 0000000..603c4af --- /dev/null +++ b/s2/builder_graph_edge_processor_test.go @@ -0,0 +1,824 @@ +// 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 ( + "slices" + "testing" +) + +func TestGraphEdgeProcessorStableLessThan(t *testing.T) { + tests := []struct { + name string + a graphEdge + b graphEdge + aInputID int32 + bInputID int32 + want bool + }{ + { + name: "a < b lexicographically", + a: graphEdge{first: 1, second: 2}, + b: graphEdge{first: 2, second: 3}, + aInputID: 1, + bInputID: 2, + want: true, + }, + { + name: "a > b lexicographically", + a: graphEdge{first: 3, second: 4}, + b: graphEdge{first: 1, second: 2}, + aInputID: 1, + bInputID: 2, + want: false, + }, + { + name: "a == b lexicographically, a.inputID < b.inputID", + a: graphEdge{first: 1, second: 2}, + b: graphEdge{first: 1, second: 2}, + aInputID: 1, + bInputID: 2, + want: true, + }, + { + name: "a == b lexicographically, a.inputID > b.inputID", + a: graphEdge{first: 1, second: 2}, + b: graphEdge{first: 1, second: 2}, + aInputID: 3, + bInputID: 2, + want: false, + }, + { + name: "a == b lexicographically, a.inputID == b.inputID", + a: graphEdge{first: 1, second: 2}, + b: graphEdge{first: 1, second: 2}, + aInputID: 5, + bInputID: 5, + want: false, + }, + { + name: "first vertices equal, second vertices different", + a: graphEdge{first: 1, second: 2}, + b: graphEdge{first: 1, second: 3}, + aInputID: 1, + bInputID: 2, + want: true, + }, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + got := stableLessThan(test.a, test.b, test.aInputID, test.bInputID) + if got != test.want { + t.Errorf("stableLessThan() = %v, want %v", got, test.want) + } + }) + } +} + +func TestGraphEdgeProcessorAddEdge(t *testing.T) { + tests := []struct { + name string + edge graphEdge + inputEdgeIDSetID int32 + wantEdges int + wantInputIDs int + }{ + { + name: "add single edge", + edge: graphEdge{first: 1, second: 2}, + inputEdgeIDSetID: 1, + wantEdges: 1, + wantInputIDs: 1, + }, + { + name: "add second edge", + edge: graphEdge{first: 2, second: 3}, + inputEdgeIDSetID: 2, + wantEdges: 2, + wantInputIDs: 2, + }, + } + + ep := &edgeProcessor{ + newEdges: make([]graphEdge, 0), + newInputIDs: make([]int32, 0), + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + ep.addEdge(test.edge, test.inputEdgeIDSetID) + if len(ep.newEdges) != test.wantEdges { + t.Errorf("addEdge() edges = %v, want %v", len(ep.newEdges), test.wantEdges) + } + if len(ep.newInputIDs) != test.wantInputIDs { + t.Errorf("addEdge() inputIDs = %v, want %v", len(ep.newInputIDs), test.wantInputIDs) + } + }) + } +} + +func TestGraphEdgeProcessorAddEdges(t *testing.T) { + tests := []struct { + name string + numEdges int + edge graphEdge + inputEdgeIDSetID int32 + wantEdges int + wantInputIDs int + }{ + { + name: "add single edge", + numEdges: 1, + edge: graphEdge{first: 1, second: 2}, + inputEdgeIDSetID: 1, + wantEdges: 1, + wantInputIDs: 1, + }, + { + name: "add multiple edges", + numEdges: 3, + edge: graphEdge{first: 1, second: 2}, + inputEdgeIDSetID: 7, + wantEdges: 3, + wantInputIDs: 3, + }, + { + name: "add zero edges", + numEdges: 0, + edge: graphEdge{first: 1, second: 2}, + inputEdgeIDSetID: 8, + wantEdges: 0, // Should remain unchanged + wantInputIDs: 0, // Should remain unchanged + }, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + ep := &edgeProcessor{ + newEdges: make([]graphEdge, 0), + newInputIDs: make([]int32, 0), + } + + ep.addEdges(test.numEdges, test.edge, test.inputEdgeIDSetID) + if len(ep.newEdges) != test.wantEdges { + t.Errorf("addEdges() edges = %v, want %v", len(ep.newEdges), test.wantEdges) + } + + if len(ep.newInputIDs) != test.wantInputIDs { + t.Errorf("addEdges() inputIDs = %v, want %v", len(ep.newInputIDs), test.wantInputIDs) + } + + // addEdges uses the same inputEdgeIDSetID for each repeated edge. Ensure + // all the added ids match. + for k, v := range ep.newInputIDs { + if v != test.inputEdgeIDSetID { + t.Errorf("in addEdges, newInputIDs[%d] = %d, want %d", k, v, test.inputEdgeIDSetID) + } + } + }) + } +} + +func TestGraphEdgeProcessorHandleDegenerateEdge(t *testing.T) { + tests := []struct { + name string + edge graphEdge + options *graphOptions + outBegin int + outEnd int + nOut int + nIn int + inBegin int + in int + wantErr bool + wantEdges int + }{ + { + name: "discard degenerate edges", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesDiscard, + }, + outBegin: 0, + outEnd: 1, + nOut: 1, + nIn: 1, + inBegin: 0, + in: 1, + wantErr: false, + wantEdges: 0, + }, + { + name: "keep degenerate edges with merge", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesKeep, + duplicateEdges: duplicateEdgesMerge, + edgeType: edgeTypeDirected, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + inBegin: 0, + in: 2, + wantErr: false, + wantEdges: 1, + }, + { + name: "keep degenerate edges without merge", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesKeep, + duplicateEdges: duplicateEdgesKeep, + edgeType: edgeTypeDirected, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + inBegin: 0, + in: 2, + wantErr: false, + wantEdges: 2, + }, + { + name: "discard excess degenerate edges with incident edges", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesDiscardExcess, + edgeType: edgeTypeDirected, + }, + outBegin: 1, + outEnd: 2, + nOut: 1, + nIn: 1, + inBegin: 1, + in: 2, + wantErr: false, + wantEdges: 0, + }, + { + name: "discard excess degenerate edges without incident edges", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesDiscardExcess, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + inBegin: 0, + in: 2, + wantErr: false, + wantEdges: 1, + }, + { + name: "undirected degenerate edges with require siblings", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesKeep, + edgeType: edgeTypeUndirected, + siblingPairs: siblingPairsRequire, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + inBegin: 0, + in: 2, + wantErr: false, + wantEdges: 1, + }, + { + name: "undirected degenerate edges with create siblings", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesKeep, + edgeType: edgeTypeUndirected, + siblingPairs: siblingPairsCreate, + duplicateEdges: duplicateEdgesKeep, + }, + outBegin: 0, + outEnd: 4, + nOut: 4, + nIn: 4, + inBegin: 0, + in: 4, + wantErr: false, + wantEdges: 2, + }, + { + name: "inconsistent degenerate edges", + edge: graphEdge{first: 1, second: 1}, + options: &graphOptions{ + degenerateEdges: degenerateEdgesKeep, + }, + outBegin: 0, + outEnd: 1, + nOut: 1, + nIn: 2, // Mismatched counts + inBegin: 0, + in: 2, + wantErr: true, + }, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + // Create edges and inputIDs arrays with the correct size + edges := make([]graphEdge, test.outEnd) + inputIDs := make([]int32, test.outEnd) + for i := range edges { + edges[i] = test.edge + inputIDs[i] = int32(i + 1) + } + + ep := newEdgeProcessor(test.options, edges, inputIDs, newIDSetLexicon()) + // Add the input IDs to the lexicon. + for _, id := range inputIDs { + ep.idSetLexicon.add(id) + } + + gotErr := ep.handleDegenerateEdge(test.edge, test.outBegin, test.outEnd, + test.nOut, test.nIn, test.inBegin, test.in) + if (gotErr != nil) != test.wantErr { + t.Errorf("handleDegenerateEdge() error = %v, wantErr %v", + gotErr, test.wantErr) + } + if !test.wantErr && len(ep.newEdges) != test.wantEdges { + t.Errorf("handleDegenerateEdge() added %d edges, want %d", len(ep.newEdges), test.wantEdges) + } + }) + } +} + +func TestGraphEdgeProcessorHandleNormalEdge(t *testing.T) { + tests := []struct { + name string + edge graphEdge + options *graphOptions + outBegin int + outEnd int + nOut int + nIn int + wantErr bool + wantEdges int + }{ + { + name: "keep sibling pairs with merge", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsKeep, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 1, + wantErr: false, + wantEdges: 1, + }, + { + name: "keep sibling pairs without merge", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsKeep, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesKeep, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 1, + wantErr: false, + wantEdges: 2, + }, + { + name: "discard sibling pairs directed balanced", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsDiscard, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + wantErr: false, + wantEdges: 0, + }, + { + name: "discard sibling pairs directed unbalanced", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsDiscard, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 3, + nOut: 3, + nIn: 1, + wantErr: false, + wantEdges: 1, + }, + { + name: "discard sibling pairs undirected even", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsDiscard, + edgeType: edgeTypeUndirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + wantErr: false, + wantEdges: 0, + }, + { + name: "discard sibling pairs undirected odd", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsDiscard, + edgeType: edgeTypeUndirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 3, + nOut: 3, + nIn: 3, + wantErr: false, + wantEdges: 1, + }, + { + name: "discard excess sibling pairs directed balanced", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsDiscardExcess, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + wantErr: false, + wantEdges: 1, + }, + { + name: "discard excess sibling pairs directed unbalanced", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsDiscardExcess, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 3, + nOut: 3, + nIn: 1, + wantErr: false, + wantEdges: 1, + }, + { + name: "require sibling pairs directed balanced", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsRequire, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 2, + wantErr: false, + wantEdges: 1, + }, + { + name: "require sibling pairs directed unbalanced", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsRequire, + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + }, + outBegin: 0, + outEnd: 2, + nOut: 2, + nIn: 1, + wantErr: true, + }, + { + name: "create sibling pairs undirected", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairsCreate, + edgeType: edgeTypeUndirected, + duplicateEdges: duplicateEdgesKeep, + }, + outBegin: 0, + outEnd: 4, + nOut: 4, + nIn: 2, + wantErr: false, + wantEdges: 2, + }, + { + name: "invalid sibling pairs option", + edge: graphEdge{first: 1, second: 2}, + options: &graphOptions{ + siblingPairs: siblingPairs(255), // Use max uint8 value as invalid + edgeType: edgeTypeDirected, + }, + outBegin: 0, + outEnd: 1, + nOut: 1, + nIn: 1, + wantErr: true, + }, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + // Create edges and inputIDs arrays with the correct size + edges := make([]graphEdge, test.outEnd) + inputIDs := make([]int32, test.outEnd) + for i := range edges { + edges[i] = test.edge + inputIDs[i] = int32(i + 1) + } + + ep := newEdgeProcessor(test.options, edges, inputIDs, newIDSetLexicon()) + // Add the input IDs to the lexicon. + for _, id := range inputIDs { + ep.idSetLexicon.add(id) + } + + gotErr := ep.handleNormalEdge(test.edge, test.outBegin, test.outEnd, test.nOut, test.nIn) + if (gotErr != nil) != test.wantErr { + t.Errorf("handleNormalEdge() error = %v, wantErr %v", gotErr, test.wantErr) + } + if !test.wantErr && len(ep.newEdges) != test.wantEdges { + t.Errorf("handleNormalEdge() added %d edges, want %d", len(ep.newEdges), test.wantEdges) + } + }) + } +} + +func TestGraphEdgeProcessorMergeInputIDs(t *testing.T) { + tests := []struct { + name string + edges []graphEdge + inputIDs []int32 + outBegin int + outEnd int + wantInputIDSet []int32 + }{ + { + name: "single edge", + edges: []graphEdge{ + {first: 1, second: 2}, + }, + inputIDs: []int32{1}, + outBegin: 0, + outEnd: 1, + wantInputIDSet: []int32{1}, + }, + { + name: "multiple edges with same input ID should reduce to 1 output", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 1, second: 2}, + }, + inputIDs: []int32{1, 1}, + outBegin: 0, + outEnd: 2, + wantInputIDSet: []int32{1}, + }, + { + name: "multiple edges with different input IDs should keep distinct ids", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 1, second: 2}, + {first: 1, second: 2}, + }, + inputIDs: []int32{1, 2, 3}, + outBegin: 0, + outEnd: 3, + wantInputIDSet: []int32{1, 2, 3}, + }, + { + name: "subset of edges should return the smaller portion", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 1, second: 2}, + {first: 1, second: 2}, + }, + inputIDs: []int32{1, 2, 3}, + outBegin: 1, + outEnd: 3, + wantInputIDSet: []int32{2, 3}, + }, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + ep := newEdgeProcessor(defaultGraphOptions(), test.edges, test.inputIDs, newIDSetLexicon()) + // Add the input IDs to the lexicon + for _, id := range test.inputIDs { + ep.idSetLexicon.add(id) + } + + merged := ep.mergeInputIDs(test.outBegin, test.outEnd) + + // Get the actual set of IDs from the lexicon + got := ep.idSetLexicon.idSet(merged) + + // Sort both slices for comparison + slices.Sort(got) + slices.Sort(test.wantInputIDSet) + + if !slices.Equal(got, test.wantInputIDSet) { + t.Errorf("mergeInputIDs() = %v, want %v", got, test.wantInputIDSet) + } + }) + } +} + +func TestGraphEdgeProcessorRun(t *testing.T) { + tests := []struct { + name string + edges []graphEdge + inputIDs []int32 + options *graphOptions + wantErr bool + wantEdges int + }{ + { + name: "empty graph", + edges: []graphEdge{}, + inputIDs: []int32{}, + options: defaultGraphOptions(), + wantErr: false, + wantEdges: 0, + }, + { + name: "single edge", + edges: []graphEdge{ + {first: 1, second: 2}, + }, + inputIDs: []int32{1}, + options: defaultGraphOptions(), + wantErr: false, + wantEdges: 1, + }, + { + name: "degenerate edge with discard", + edges: []graphEdge{ + {first: 1, second: 1}, + }, + inputIDs: []int32{1}, + options: &graphOptions{ + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesKeep, + degenerateEdges: degenerateEdgesDiscard, + siblingPairs: siblingPairsKeep, + }, + wantErr: false, + wantEdges: 0, + }, + { + name: "duplicate edges with merge", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 1, second: 2}, + }, + inputIDs: []int32{1, 2}, + options: &graphOptions{ + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + degenerateEdges: degenerateEdgesKeep, + siblingPairs: siblingPairsKeep, + }, + wantErr: false, + wantEdges: 1, + }, + { + name: "sibling pairs with discard", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 2, second: 1}, + }, + inputIDs: []int32{1, 2}, + options: &graphOptions{ + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesKeep, + degenerateEdges: degenerateEdgesKeep, + siblingPairs: siblingPairsDiscard, + }, + wantErr: false, + wantEdges: 0, + }, + { + name: "undirected edges with require siblings", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 2, second: 1}, + }, + inputIDs: []int32{1, 2}, + options: &graphOptions{ + edgeType: edgeTypeUndirected, + duplicateEdges: duplicateEdgesKeep, + degenerateEdges: degenerateEdgesKeep, + siblingPairs: siblingPairsRequire, + }, + wantErr: true, + wantEdges: 0, + }, + { + name: "undirected edges with create siblings", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 2, second: 1}, + }, + inputIDs: []int32{1, 2}, + options: &graphOptions{ + edgeType: edgeTypeUndirected, + duplicateEdges: duplicateEdgesKeep, + degenerateEdges: degenerateEdgesKeep, + siblingPairs: siblingPairsCreate, + }, + wantErr: false, + wantEdges: 2, + }, + { + name: "require siblings with missing sibling", + edges: []graphEdge{ + {first: 1, second: 2}, + }, + inputIDs: []int32{1}, + options: &graphOptions{ + edgeType: edgeTypeUndirected, + duplicateEdges: duplicateEdgesKeep, + degenerateEdges: degenerateEdgesKeep, + siblingPairs: siblingPairsRequire, + }, + wantErr: true, + wantEdges: 0, + }, + { + name: "multiple edges with various options", + edges: []graphEdge{ + {first: 1, second: 2}, + {first: 2, second: 1}, + {first: 1, second: 2}, + {first: 3, second: 3}, + }, + inputIDs: []int32{1, 2, 3, 4}, + options: &graphOptions{ + edgeType: edgeTypeDirected, + duplicateEdges: duplicateEdgesMerge, + degenerateEdges: degenerateEdgesDiscard, + siblingPairs: siblingPairsDiscardExcess, + }, + wantErr: false, + wantEdges: 1, + }, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + lexicon := newIDSetLexicon() + ep := newEdgeProcessor(test.options, test.edges, test.inputIDs, lexicon) + err := ep.Run() + if (err != nil) != test.wantErr { + t.Errorf("Run() error = %v, wantErr %v", err, test.wantErr) + } + if !test.wantErr && len(ep.edges) != test.wantEdges { + t.Errorf("Run() produced %d edges, want %d", len(ep.edges), test.wantEdges) + } + }) + } +} diff --git a/s2/builder_graph_options.go b/s2/builder_graph_options.go new file mode 100644 index 0000000..360c9fb --- /dev/null +++ b/s2/builder_graph_options.go @@ -0,0 +1,179 @@ +// 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 + +// degenerateEdges controls how degenerate edges (i.e., an edge from a vertex to +// itself) are handled. Such edges may be present in the input, or they may be +// created when both endpoints of an edge are snapped to the same output vertex. +// The options available are: +type degenerateEdges uint8 + +const ( + // degenerateEdgesDiscard discards all degenerate edges. This is useful for + // layers that/do not support degeneracies, such as PolygonLayer. + degenerateEdgesDiscard degenerateEdges = iota + // degenerateEdgesDiscardExcess discards all degenerate edges that are + // connected to/non-degenerate edges and merges any remaining + // duplicate/degenerate edges. This is useful for simplifying/polygons + // while ensuring that loops that collapse to a/single point do not disappear. + degenerateEdgesDiscardExcess + // degenerateEdgesKeep: Keeps all degenerate edges. Be aware that this + // may create many redundant edges when simplifying geometry (e.g., a + // polyline of the form AABBBBBCCCCCCDDDD). degenerateEdgesKeep is mainly + // useful for algorithms that require an output edge for every input edge. + degenerateEdgesKeep +) + +// duplicateEdges controls how duplicate edges (i.e., edges that are present +// multiple times) are handled. Such edges may be present in the input, or they +// can be created when vertices are snapped together. When several edges are +// merged, the result is a single edge labelled with all of the original input +// edge ids. +type duplicateEdges uint8 + +const ( + duplicateEdgesMerge duplicateEdges = iota + duplicateEdgesKeep +) + +// siblingPairs controls how sibling edge pairs (i.e., pairs consisting +// of an edge and its reverse edge) are handled. Layer types that +// define an interior (e.g., polygons) normally discard such edge pairs +// since they do not affect the result (i.e., they define a "loop" with +// no interior). +// +// If edgeType is edgeTypeUndirected, a sibling edge pair is considered +// to consist of four edges (two duplicate edges and their siblings), since +// only two of these four edges will be used in the final output. +// +// Furthermore, since the options REQUIRE and CREATE guarantee that all +// edges will have siblings, Builder implements these options for +// undirected edges by discarding half of the edges in each direction and +// changing the edgeType to edgeTypeDirected. For example, two +// undirected input edges between vertices A and B would first be converted +// into two directed edges in each direction, and then one edge of each pair +// would be discarded leaving only one edge in each direction. +// +// Degenerate edges are considered not to have siblings. If such edges are +// present, they are passed through unchanged by siblingPairsDiscard. For +// siblingPairsRequire or siblingPairsCreate with undirected edges, the +// number of copies of each degenerate edge is reduced by a factor of two. +// Any of the options that discard edges (DISCARD, DISCARDEXCESS, and +// REQUIRE/CREATE in the case of undirected edges) have the side effect that +// when duplicate edges are present, all of the corresponding edge labels +// are merged together and assigned to the remaining edges. (This avoids +// the problem of having to decide which edges are discarded.) Note that +// this merging takes place even when all copies of an edge are kept. For +// example, consider the graph {AB1, AB2, AB3, BA4, CD5, CD6} (where XYn +// denotes an edge from X to Y with label "n"). With siblingPairsDiscard, +// we need to discard one of the copies of AB. But which one? Rather than +// choosing arbitrarily, instead we merge the labels of all duplicate edges +// (even ones where no sibling pairs were discarded), yielding {AB123, +// AB123, CD45, CD45} (assuming that duplicate edges are being kept). +// Notice that the labels of duplicate edges are merged even if no siblings +// were discarded (such as CD5, CD6 in this example), and that this would +// happen even with duplicate degenerate edges (e.g. the edges EE7, EE8). +type siblingPairs uint8 + +const ( + // siblingPairsDiscard discards all sibling edge pairs. + siblingPairsDiscard siblingPairs = iota + // siblingPairsDiscardExcess is like siblingPairsDiscard, except that a + // single sibling pair is kept if the result would otherwise be empty. + // This is useful for polygons with degeneracies (LaxPolygon), and for + // simplifying polylines while ensuring that they are not split into + // multiple disconnected pieces. + siblingPairsDiscardExcess + // siblingPairsKeep keeps sibling pairs. This can be used to create + // polylines that double back on themselves, or degenerate loops (with + // a layer type such as LaxPolygon). + siblingPairsKeep + // siblingPairsRequire requires that all edges have a sibling (and returns + // an error otherwise). This is useful with layer types that create a + // collection of adjacent polygons (a polygon mesh). + siblingPairsRequire + // siblingPairsCreate ensures that all edges have a sibling edge by + // creating them if necessary. This is useful with polygon meshes where + // the input polygons do not cover the entire sphere. Such edges always + // have an empty set of labels and do not have an associated InputEdgeID. + siblingPairsCreate +) + +// graphOptions is only needed by Layer implementations. A layer is +// responsible for assembling an Graph of snapped edges into the +// desired output format (e.g., an Polygon). The graphOptions allows +// each Layer type to specify requirements on its input graph: for example, if +// degenerateEdgesDiscard is specified, then Builder will ensure that all +// degenerate edges are removed before passing the graph to Layer's Build +// method. +type graphOptions struct { + // Specifies whether the Builder input edges should be treated as + // undirected. If true, then all input edges are duplicated into pairs + // consisting of an edge and a sibling (reverse) edge. Note that the + // automatically created sibling edge has an empty set of labels and does + // not have an associated InputEdgeId. + // + // The layer implementation is responsible for ensuring that exactly one + // edge from each pair is used in the output, i.e. *only half* of the graph + // edges will be used. (Note that some values of the siblingPairs option + // automatically take care of this issue by removing half of the edges and + // changing edgeType to Directed.) + // + // DEFAULT: edgeTypeDirected + edgeType edgeType + + // DEFAULT: degenerateEdgesKeep + degenerateEdges degenerateEdges + + // DEFAULT: duplicateEdgesKeep + duplicateEdges duplicateEdges + + // DEFAULT: siblingPairsKeep + siblingPairs siblingPairs + + // This is a specialized option that is only needed by clients that want to + // work with the graphs for multiple layers at the same time (e.g., in order + // to check whether the same edge is present in two different graphs). [Note + // that if you need to do this, usually it is easier just to build a single + // graph with suitable edge labels.] + // + // When there are a large number of layers, then by default Builder builds + // a minimal subgraph for each layer containing only the vertices needed by + // the edges in that layer. This ensures that layer types that iterate over + // the vertices run in time proportional to the size of that layer rather + // than the size of all layers combined. (For example, if there are a + // million layers with one edge each, then each layer would be passed a + // graph with 2 vertices rather than 2 million vertices.) + // + // If this option is set to false, this optimization is disabled. Instead + // the graph passed to this layer will contain the full set of vertices. + // (This is not recommended when the number of layers could be large.) + // + // DEFAULT: true + allowVertexFiltering bool +} + +// defaultGraphOptions returns a graphOptions that specify that all edges should +// be kept, since this produces the least surprising output and makes it easier +// to diagnose the problem when an option is left unspecified. +func defaultGraphOptions() *graphOptions { + return &graphOptions{ + edgeType: edgeTypeDirected, + degenerateEdges: degenerateEdgesKeep, + duplicateEdges: duplicateEdgesKeep, + siblingPairs: siblingPairsKeep, + allowVertexFiltering: true, + } +} diff --git a/s2/builder_layer.go b/s2/builder_layer.go new file mode 100644 index 0000000..7e4ee77 --- /dev/null +++ b/s2/builder_layer.go @@ -0,0 +1,32 @@ +// 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 + +// builderLayer defines the methods Layers must implement. +type builderLayer interface { + // GraphOptions returns the options defined by this layer. + GraphOptions() *graphOptions + + // Build assembles a graph of snapped edges into the geometry type + // implemented by this layer. If an error is encountered, error is + // set appropriately. + // + // Note that when there are multiple layers, the Graph object passed to all + // layers are guaranteed to be valid until the last Build() method returns. + // This makes it easier to write algorithms that gather the output graphs + // from several layers and process them all at once (such as + // closedSetNormalizer). + Build(g *graph) (bool, error) +} diff --git a/s2/builder_options.go b/s2/builder_options.go new file mode 100644 index 0000000..55e6c25 --- /dev/null +++ b/s2/builder_options.go @@ -0,0 +1,251 @@ +// 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 "github.com/golang/geo/s1" + +// polylineType Indicates whether polylines should be "paths" (which don't +// allow duplicate vertices, except possibly the first and last vertex) or +// "walks" (which allow duplicate vertices and edges). +type polylineType uint8 + +const ( + polylineTypePath polylineType = iota + polylineTypeWalk +) + +// edgeType indicates whether the input edges are undirected. Typically this is +// specified for each output layer (e.g., PolygonBuilderLayer). +// +// Directed edges are preferred, since otherwise the output is ambiguous. +// For example, output polygons may be the *inverse* of the intended result +// (e.g., a polygon intended to represent the world's oceans may instead +// represent the world's land masses). Directed edges are also somewhat +// more efficient. +// +// However even with undirected edges, most Builder layer types try to +// preserve the input edge direction whenever possible. Generally, edges +// are reversed only when it would yield a simpler output. For example, +// PolygonLayer assumes that polygons created from undirected edges should +// cover at most half of the sphere. Similarly, PolylineVectorBuilderLayer +// assembles edges into as few polylines as possible, even if this means +// reversing some of the "undirected" input edges. +// +// For shapes with interiors, directed edges should be oriented so that the +// interior is to the left of all edges. This means that for a polygon with +// holes, the outer loops ("shells") should be directed counter-clockwise +// while the inner loops ("holes") should be directed clockwise. Note that +// AddPolygon() follows this convention automatically. +type edgeType uint8 + +const ( + edgeTypeDirected edgeType = iota + edgeTypeUndirected +) + +type builderOptions struct { + // snapFunction holds the desired snap function. + // + // Note that if your input data includes vertices that were created using + // Intersection(), then you should use a "snapRadius" of + // at least intersectionMergeRadius, e.g. by calling + // + // options.setSnapFunction(IdentitySnapFunction(intersectionMergeRadius)); + // + // DEFAULT: IdentitySnapFunction(s1.Angle(0)) + // [This does no snapping and preserves all input vertices exactly.] + snapFunction Snapper + + // splitCrossingEdges determines how crossing edges are handled by Builder. + // If true, then detect all pairs of crossing edges and eliminate them by + // adding a new vertex at their intersection point. See also the + // AddIntersection() method which allows intersection points to be added + // selectively. + // + // When this option if true, intersectionTolerance is automatically set + // to a minimum of intersectionError (see intersectionTolerance + // for why this is necessary). Note that this means that edges can move + // by up to intersectionError even when the specified snap radius is + // zero. The exact distance that edges can move is always given by + // MaxEdgeDeviation(). + // + // Undirected edges should always be used when the output is a polygon, + // since splitting a directed loop at a self-intersection converts it into + // two loops that don't define a consistent interior according to the + // "interior is on the left" rule. (On the other hand, it is fine to use + // directed edges when defining a polygon *mesh* because in that case the + // input consists of sibling edge pairs.) + // + // Self-intersections can also arise when importing data from a 2D + // projection. You can minimize this problem by subdividing the input + // edges so that the S2 edges (which are geodesics) stay close to the + // original projected edges (which are curves on the sphere). This can + // be done using EdgeTessellator, for example. + // + // DEFAULT: false + splitCrossingEdges bool + + // intersectionTolerance specifies the maximum allowable distance between + // a vertex added by AddIntersection() and the edge(s) that it is intended + // to snap to. This method must be called before AddIntersection() can be + // used. It has the effect of increasing the snap radius for edges (but not + // vertices) by the given distance. + // + // The intersection tolerance should be set to the maximum error in the + // intersection calculation used. For example, if Intersection() + // is used then the error should be set to intersectionError. If + // PointOnLine is used then the error should be set to PointOnLineError. + // If Project is used then the error should be set to + // projectPerpendicularError. If more than one method is used then the + // intersection tolerance should be set to the maximum such error. + // + // The reason this option is necessary is that computed intersection + // points are not exact. For example, Intersection(a, b, c, d) + // returns a point up to intersectionError away from the true + // mathematical intersection of the edges AB and CD. Furthermore such + // intersection points are subject to further snapping in order to ensure + // that no pair of vertices is closer than the specified snap radius. For + // example, suppose the computed intersection point X of edges AB and CD + // is 1 nanonmeter away from both edges, and the snap radius is 1 meter. + // In that case X might snap to another vertex Y exactly 1 meter away, + // which would leave us with a vertex Y that could be up to 1.000000001 + // meters from the edges AB and/or CD. This means that AB and/or CD might + // not snap to Y leaving us with two edges that still cross each other. + // + // However if the intersection tolerance is set to 1 nanometer then the + // snap radius for edges is increased to 1.000000001 meters ensuring that + // both edges snap to a common vertex even in this worst case. (Tthis + // technique does not work if the vertex snap radius is increased as well; + // it requires edges and vertices to be handled differently.) + // + // Note that this option allows edges to move by up to the given + // intersection tolerance even when the snap radius is zero. The exact + // distance that edges can move is always given by maxEdgeDeviation() + // defined above. + // + // When splitCrossingEdges is true, the intersection tolerance is + // automatically set to a minimum of intersectionError. A larger + // value can be specified by calling this method explicitly. + // + // DEFAULT: s1.Angle(0) + intersectionTolerance s1.Angle + + // simplifyEdgeChains determines if the output geometry should be simplified + // by replacing nearly straight chains of short edges with a single long edge. + // + // The combined effect of snapping and simplifying will not change the + // input by more than the guaranteed tolerances (see the list documented + // with the SnapFunction class). For example, simplified edges are + // guaranteed to pass within snapRadius() of the *original* positions of + // all vertices that were removed from that edge. This is a much tighter + // guarantee than can be achieved by snapping and simplifying separately. + // + // However, note that this option does not guarantee idempotency. In + // other words, simplifying geometry that has already been simplified once + // may simplify it further. (This is unavoidable, since tolerances are + // measured with respect to the original geometry, which is no longer + // available when the geometry is simplified a second time.) + // + // When the output consists of multiple layers, simplification is + // guaranteed to be consistent: for example, edge chains are simplified in + // the same way across layers, and simplification preserves topological + // relationships between layers (e.g., no crossing edges will be created). + // Note that edge chains in different layers do not need to be identical + // (or even have the same number of vertices, etc) in order to be + // simplified together. All that is required is that they are close + // enough together so that the same simplified edge can meet all of their + // individual snapping guarantees. + // + // Note that edge chains are approximated as parametric curves rather than + // point sets. This means that if an edge chain backtracks on itself (for + // example, ABCDEFEDCDEFGH) then such backtracking will be preserved to + // within snapRadius() (for example, if the preceding point were all in a + // straight line then the edge chain would be simplified to ACFCFH, noting + // that C and F have degree > 2 and therefore can't be simplified away). + // + // Simplified edges are assigned all labels associated with the edges of + // the simplified chain. + // + // For this option to have any effect, a SnapFunction with a non-zero + // snapRadius() must be specified. Also note that vertices specified + // using ForceVertex are never simplified away. + // + // DEFAULT: false + simplifyEdgeChains bool + + // idempotent determines if snapping occurs only when the input geometry + // does not already meet the Builder output guarantees (see the Snapper + // type description for details). This means that if all input vertices + // are at snapped locations, all vertex pairs are separated by at least + // MinVertexSeparation(), and all edge-vertex pairs are separated by at + // least MinEdgeVertexSeparation(), then no snapping is done. + // + // If false, then all vertex pairs and edge-vertex pairs closer than + // "SnapRadius" will be considered for snapping. This can be useful, for + // example, if you know that your geometry contains errors and you want to + // make sure that features closer together than "SnapRadius" are merged. + // + // This option is automatically turned off when simplifyEdgeChains is true + // since simplifying edge chains is never guaranteed to be idempotent. + // + // DEFAULT: true + idempotent bool +} + +// defaultBuilderOptions returns a new instance with the proper defaults. +func defaultBuilderOptions() *builderOptions { + return &builderOptions{ + snapFunction: NewIdentitySnapper(0), + splitCrossingEdges: false, + intersectionTolerance: s1.Angle(0), + simplifyEdgeChains: false, + idempotent: true, + } +} + +// edgeSnapRadius reports the maximum distance from snapped edge vertices to +// the original edge. This is the same as SnapFunction().SnapRadius() except +// when splitCrossingEdges is true (see below), in which case the edge snap +// radius is increased by intersectionError. +func (o builderOptions) edgeSnapRadius() s1.Angle { + return o.snapFunction.SnapRadius() + o.intersectionTolerance +} + +// maxEdgeDeviation returns maximum distance that any point along an edge can +// move when snapped. It is slightly larger than edgeSnapRadius() because when +// a geodesic edge is snapped, the edge center moves further than its endpoints. +// Builder ensures that this distance is at most 10% larger than +// edgeSnapRadius(). +func (o builderOptions) maxEdgeDeviation() s1.Angle { + // We want maxEdgeDeviation to be large enough compared to SnapRadius() + // such that edge splitting is rare. + // + // Using spherical trigonometry, if the endpoints of an edge of length L + // move by at most a distance R, the center of the edge moves by at most + // asin(sin(R) / cos(L / 2)). Thus the (maxEdgeDeviation / SnapRadius) + // ratio increases with both the snap radius R and the edge length L. + // + // We arbitrarily limit the edge deviation to be at most 10% more than the + // snap radius. With the maximum allowed snap radius of 70 degrees, this + // means that edges up to 30.6 degrees long are never split. For smaller + // snap radii, edges up to 49 degrees long are never split. (Edges of any + // length are not split unless their endpoints move far enough so that the + // actual edge deviation exceeds the limit; in practice, splitting is rare + // even with long edges.) Note that it is always possible to split edges + // when maxEdgeDeviation() is exceeded; see maybeAddExtraSites(). + // + // TODO(rsned): What should we do when snapFunction.SnapRadius() > maxSnapRadius); + return maxEdgeDeviationRatio * o.edgeSnapRadius() +}