diff --git a/go.mod b/go.mod index 9ce5f35..4ed1d55 100644 --- a/go.mod +++ b/go.mod @@ -1,3 +1,7 @@ module github.com/golang/geo -go 1.18 +go 1.22.0 + +toolchain go1.22.5 + +require golang.org/x/exp v0.0.0-20240909161429-701f63a606c0 diff --git a/go.sum b/go.sum new file mode 100644 index 0000000..b314e3d --- /dev/null +++ b/go.sum @@ -0,0 +1,2 @@ +golang.org/x/exp v0.0.0-20240909161429-701f63a606c0 h1:e66Fs6Z+fZTbFBAxKfP3PALWBtpfqks2bwGcexMxgtk= +golang.org/x/exp v0.0.0-20240909161429-701f63a606c0/go.mod h1:2TbTHSBQa924w8M6Xs1QcRcFwyucIwBGpK1p2f1YFFY= diff --git a/s2/chain_interpolation_query.go b/s2/chain_interpolation_query.go new file mode 100644 index 0000000..f6b4a64 --- /dev/null +++ b/s2/chain_interpolation_query.go @@ -0,0 +1,336 @@ +package s2 + +import ( + "errors" + "slices" + + "github.com/golang/geo/s1" +) + +var ( + // ErrEmptyChain is returned by ChainInterpolationQuery when the query + // contains no edges. + ErrEmptyChain = errors.New("empty chain") + + // ErrInvalidDivisionsCount is returned by ChainInterpolationQuery when + // divisionsCount is less than the number of edges in the shape. + ErrInvalidDivisionsCount = errors.New("invalid divisions count") + + // ErrInvalidIndexes is returned by ChainInterpolationQuery when + // start or end indexes are invalid. + ErrInvalidIndexes = errors.New("invalid indexes") +) + +// ChainInterpolationQuery is a helper struct for querying points on Shape's +// edges by spherical distance. The distance is computed cumulatively along the +// edges contained in the shape, using the order in which the edges are stored +// by the Shape object. +type ChainInterpolationQuery struct { + Shape Shape + ChainID int + cumulativeValues []s1.Angle + firstEdgeID int + lastEdgeID int +} + +// InitChainInterpolationQuery initializes and conctructs a ChainInterpolationQuery. +// If a particular chain id is specified at the query initialization, then the +// distance values are computed along that single chain, which allows per-chain +// interpolation. If no chain is specified, then the interpolated point as a +// function of distance is discontinuous at chain boundaries. Using multiple +// chains can be used in such algorithms as quasi-random sampling along the +// total span of a multiline. +// +// Once the query object is initialized, the complexity of each subsequent query +// is O( log(number of edges) ). The complexity of the initialization and the +// memory footprint of the query object are both O(number of edges). +func InitChainInterpolationQuery(shape Shape, chainID int) ChainInterpolationQuery { + if shape == nil || chainID >= shape.NumChains() { + return ChainInterpolationQuery{nil, 0, nil, 0, 0} + } + + var firstEdgeID, lastEdgeID int + var cumulativeValues []s1.Angle + + if chainID >= 0 { + // If a valid chain id was provided, then the range of edge ids is defined + // by the start and the length of the chain. + chain := shape.Chain(chainID) + firstEdgeID = chain.Start + lastEdgeID = firstEdgeID + chain.Length - 1 + } else { + // If no valid chain id was provided then we use the whole range of shape's + // edge ids. + firstEdgeID = 0 + lastEdgeID = shape.NumEdges() - 1 + } + + var cumulativeAngle s1.Angle + + for i := firstEdgeID; i <= lastEdgeID; i++ { + cumulativeValues = append(cumulativeValues, cumulativeAngle) + edge := shape.Edge(i) + edgeAngle := edge.V0.Angle(edge.V1.Vector) + cumulativeAngle += edgeAngle + } + + if len(cumulativeValues) != 0 { + cumulativeValues = append(cumulativeValues, cumulativeAngle) + } + return ChainInterpolationQuery{shape, chainID, cumulativeValues, firstEdgeID, lastEdgeID} +} + +// Gets the total length of the chain(s), which corresponds to the distance at +// the end vertex of the last edge of the chain(s). Returns zero length for +// shapes containing no edges. +func (s ChainInterpolationQuery) GetLength() (s1.Angle, error) { + // The total length equals the cumulative value at the end of the last + // edge, if there is at least one edge in the shape. + if len(s.cumulativeValues) == 0 { + return 0, ErrEmptyChain + } + return s.cumulativeValues[len(s.cumulativeValues)-1], nil +} + +// Returns the cumulative length along the edges being interpolated up to the +// end of the given edge ID. Returns s1.InfAngle() if the edge +// ID does not lie within the set of edges being interpolated. Returns +// ErrEmptyChain if the ChainInterpolationQuery is empty. +func (s ChainInterpolationQuery) GetLengthAtEdgeEnd(edgeID int) (s1.Angle, error) { + if len(s.cumulativeValues) == 0 { + return 0, ErrEmptyChain + } + + if edgeID < s.firstEdgeID || edgeID > s.lastEdgeID { + return s1.InfAngle(), nil + } + + return s.cumulativeValues[edgeID-s.firstEdgeID+1], nil +} + +// Computes the Point located at the given distance along the edges from the +// first vertex of the first edge. Also computes the edge id and the actual +// distance corresponding to the resulting point. +// +// This method returns a valid result if the query has been initialized with +// at least one edge. +// +// If the input distance exceeds the total length, then the resulting point is +// the end vertex of the last edge, and the resulting distance is set to the +// total length. +// +// If there are one or more degenerate (zero-length) edges corresponding to +// the given distance, then the resulting point is located on the first of +// these edges. +func (s ChainInterpolationQuery) AtDistance(inputDistance s1.Angle) (point Point, edgeID int, distance s1.Angle, err error) { + if len(s.cumulativeValues) == 0 { + return point, 0, 0, ErrEmptyChain + } + + distance = inputDistance + + position, found := slices.BinarySearch(s.cumulativeValues, inputDistance) + + if position <= 0 { + // Corner case: the first vertex of the shape at distance = 0. + return s.Shape.Edge(s.firstEdgeID).V0, s.firstEdgeID, s.cumulativeValues[0], nil + } else if (found && position == len(s.cumulativeValues)-1) || (!found && position >= len(s.cumulativeValues)) { + // Corner case: the input distance is greater than the total length, hence + // we snap the result to the last vertex of the shape at distance = total + // length. + return s.Shape.Edge(s.lastEdgeID).V1, s.lastEdgeID, s.cumulativeValues[len(s.cumulativeValues)-1], nil + } else { + // Obtain the edge index and compute the interpolated result from the edge + // vertices. + edgeID = max(position+s.firstEdgeID-1, 0) + edge := s.Shape.Edge(edgeID) + point = GetPointOnLine(edge.V0, edge.V1, inputDistance-s.cumulativeValues[max(0, position-1)]) + } + + return point, edgeID, distance, nil +} + +// Similar to the above function, but takes the normalized fraction of the +// distance as input, with inputFraction = 0 corresponding to the beginning of the +// shape or chain and inputFraction = 1 to the end. Forwards the call to +// AtDistance(). A small precision loss may occur due to converting the +// fraction to distance by multiplying it by the total length. +func (s ChainInterpolationQuery) AtFraction(inputFraction float64) (point Point, edgeID int, distance s1.Angle, err error) { + length, error := s.GetLength() + if error != nil { + return point, 0, 0, error + } + + return s.AtDistance(s1.Angle(inputFraction * float64(length))) +} + +// Returns the vector of points that is a slice of the chain from +// beginFraction to endFraction. If beginFraction is greater than +// endFraction, then the points are returned in reverse order. +// +// For example, Slice(0,1) returns the entire chain, Slice(0, 0.5) returns the +// first half of the chain, and Slice(1, 0.5) returns the second half of the +// chain in reverse. +// +// The endpoints of the slice are interpolated (except when coinciding with an +// existing vertex of the chain), and all the internal points are copied from +// the chain as is. +// +// If the query is either uninitialized, or initialized with a shape +// containing no edges, then an empty vector is returned. +func (s ChainInterpolationQuery) Slice(beginFraction, endFraction float64) []Point { + var points []Point + s.AddSlice(beginFraction, endFraction, &points) + return points +} + +// Returns the vector of points that is a slice of the chain from +// beginFraction to endFraction. If beginFraction is greater than +// endFraction, then the points are returned in reverse order. +// +// For example, Slice(0,1) returns the entire chain, Slice(0, 0.5) returns the +// first half of the chain, and Slice(1, 0.5) returns the second half of the +// chain in reverse. +// +// The endpoints of the slice are interpolated (except when coinciding with an +// existing vertex of the chain), and all the internal points are copied from +// the chain as is. +// +// divisions is the number of segments to divide the polyline into. +// divisions must be >= len(Slice(beginFraction, endFraction)). +// +// If the query is either uninitialized, or initialized with a shape +// containing no edges, then an empty vector is returned. +func (s ChainInterpolationQuery) SliceDivided(beginFraction, endFraction float64, divisions int) []Point { + var points []Point + s.AddDividedSlice(beginFraction, endFraction, &points, divisions) + return points +} + +// Appends the chain slice from beginFraction to endFraction to the given +// slice. If beginFraction is greater than endFraction, then the points are +// appended in reverse order. If the query is either uninitialized, or +// initialized with a shape containing no edges, then no points are appended. +func (s ChainInterpolationQuery) AddSlice(beginFraction, endFraction float64, points *[]Point) { + if len(s.cumulativeValues) == 0 { + return + } + + reverse := beginFraction > endFraction + if reverse { + // Swap the begin and end fractions so that we can iterate in ascending order. + beginFraction, endFraction = endFraction, beginFraction + } + + atBegin, beginEdgeID, _, err := s.AtFraction(beginFraction) + if err != nil { + return + } + *points = append(*points, atBegin) + lastPoint := atBegin + + atEnd, endEdgeID, _, err := s.AtFraction(endFraction) + if err != nil { + return + } + + // Copy the internal points from the chain. + for edgeID := beginEdgeID; edgeID < endEdgeID; edgeID++ { + edge := s.Shape.Edge(edgeID) + if lastPoint != edge.V1 { + lastPoint = edge.V1 + *points = append(*points, lastPoint) + } + } + *points = append(*points, atEnd) + + // Reverse the slice if necessary. + if reverse { + slices.Reverse(*points) + } +} + +// Appends the slice from beginFraction to endFraction to the given +// slice. If beginFraction is greater than endFraction, then the points are +// appended in reverse order. If the query is either uninitialized, or +// initialized with a shape containing no edges, then no points are appended. +// divisions is the number of segments to divide the polyline into. +// divisions must be greater or equal of NumEdges of Shape. +// A polyline is divided into segments of equal length, and then edges are added to the slice. +func (s ChainInterpolationQuery) AddDividedSlice(beginFraction, endFraction float64, points *[]Point, pointsNum int) { + if len(s.cumulativeValues) == 0 { + return + } + + pointsLength := len(*points) + + *points = append(*points, s.Slice(beginFraction, endFraction)...) + + if len(*points) > pointsNum { + *points = (*points)[0:pointsLength] + return + } else if len(*points) == pointsNum { + return + } + + *points = (*points)[0:pointsLength] + + reverse := beginFraction > endFraction + if reverse { + // Swap the begin and end fractions so that we can iterate in ascending order. + beginFraction, endFraction = endFraction, beginFraction + } + + atBegin, currentEdgeID, _, err := s.AtFraction(beginFraction) + if err != nil { + return + } + + atEnd, _, _, err := s.AtFraction(endFraction) + if err != nil { + return + } + + // divisionsExcludingEdges := pointsNum - len(slice) + + *points = append(*points, atBegin) + + // // Copy the internal points from the chain. + for fraction := beginFraction + (endFraction-beginFraction)/float64(pointsNum-1); fraction < endFraction; fraction += (endFraction - beginFraction) / float64(pointsNum-1) { + atFraction, edgeID, _, err := s.AtFraction(fraction) + if err != nil { + return + } + + // If the current edge is the same as the previous edge, then skip it. + // Otherwise, append all edges in between. + if currentEdgeID != edgeID { + for i := currentEdgeID; i < edgeID; i++ { + edge := s.Shape.Edge(i) + if edge.V1 != atFraction { + if len(*points) == pointsNum-1 { + break + } + *points = append(*points, edge.V1) + } + } + currentEdgeID = edgeID + continue + } else if edge := s.Shape.Edge(edgeID); edge.V1.approxEqual(atFraction, epsilon) { + *points = append(*points, edge.V1) + currentEdgeID++ + continue + } + if len(*points) == pointsNum-1 { + break + } + *points = append(*points, atFraction) + } + // Append last edge + *points = append(*points, atEnd) + + // Reverse the slice if necessary. + if reverse { + slices.Reverse(*points) + } +} diff --git a/s2/chain_interpolation_query_test.go b/s2/chain_interpolation_query_test.go new file mode 100644 index 0000000..266ed16 --- /dev/null +++ b/s2/chain_interpolation_query_test.go @@ -0,0 +1,826 @@ +package s2 + +import ( + "testing" + + "github.com/golang/geo/s1" +) + +const ( + latitudeB = 1. + latitudeC = 2.5 + totalLengthAbc = latitudeC + kEpsilon = 1e-8 + kEpsilonAngle = s1.Angle(kEpsilon) +) + +type result struct { + point Point + edgeID int + distance s1.Angle + err error +} + +func TestSimplePolylines(t *testing.T) { + a := PointFromLatLng(LatLngFromDegrees(0, 0)) + b := PointFromLatLng(LatLngFromDegrees(latitudeB, 0)) + c := PointFromLatLng(LatLngFromDegrees(latitudeC, 0)) + + emptyLaxPolyline := Shape(&LaxPolyline{}) + acPolyline := Shape(&LaxPolyline{vertices: []Point{a, c}}) + abcPolyline := Shape(&LaxPolyline{vertices: []Point{a, b, c}}) + bbPolyline := Shape(&LaxPolyline{vertices: []Point{b, b}}) + ccPolyline := Shape(&LaxPolyline{vertices: []Point{c}}) + + uninitializedQuery := ChainInterpolationQuery{} + emptyQuery := InitChainInterpolationQuery(emptyLaxPolyline, 0) + acQuery := InitChainInterpolationQuery(acPolyline, 0) + abcQuery := InitChainInterpolationQuery(abcPolyline, 0) + bbQuery := InitChainInterpolationQuery(bbPolyline, 0) + ccQuery := InitChainInterpolationQuery(ccPolyline, 0) + + distances := []float64{ + -1., + 0., + 1.0e-8, + latitudeB / 2, + latitudeB - 1.0e-7, + latitudeB, + latitudeB + 1.0e-5, + latitudeB + 0.5, + latitudeC - 10.e-7, + latitudeC, + latitudeC + 10.e-16, + 1.e6, + } + + groundTruth := make([]result, len(distances)) + for i, distance := range distances { + lat := max(.0, min(totalLengthAbc, distance)) + point := PointFromLatLng(LatLngFromDegrees(lat, 0)) + var edgeID int + if distance < latitudeB { + edgeID = 0 + } else { + edgeID = 1 + } + groundTruth[i] = result{point: point, edgeID: edgeID, distance: s1.Angle(distance)} + } + + for _, args := range []struct { + query ChainInterpolationQuery + want float64 + wantErr bool + }{ + {query: uninitializedQuery, want: 0, wantErr: true}, + {query: emptyQuery, want: 0, wantErr: true}, + {query: acQuery, want: totalLengthAbc, wantErr: false}, + {query: abcQuery, want: totalLengthAbc, wantErr: false}, + {query: bbQuery, want: 0, wantErr: false}, + {query: ccQuery, want: 0, wantErr: true}, + } { + length, err := args.query.GetLength() + if args.wantErr != (err != nil) { + t.Fatalf("got %v, want %v", err, args.wantErr) + } + if !float64Near(length.Degrees(), args.want, kEpsilon) { + t.Errorf("got %v, want %v", length.Degrees(), args.want) + } + } + + for _, args := range []struct { + query ChainInterpolationQuery + totalFraction float64 + wantPoint Point + wantEdgeID int + wantDistance s1.Angle + wantErr bool + }{ + {query: uninitializedQuery, totalFraction: 0, wantPoint: Point{}, wantEdgeID: 0, wantDistance: 0, wantErr: true}, + {query: emptyQuery, totalFraction: 0, wantPoint: Point{}, wantEdgeID: 0, wantDistance: 0, wantErr: true}, + {query: acQuery, totalFraction: 0, wantPoint: a, wantEdgeID: 0, wantDistance: s1.Angle(0), wantErr: false}, + {query: abcQuery, totalFraction: 0, wantPoint: a, wantEdgeID: 0, wantDistance: s1.Angle(0), wantErr: false}, + {query: bbQuery, totalFraction: 0, wantPoint: b, wantEdgeID: 0, wantDistance: s1.Angle(0), wantErr: false}, + {query: ccQuery, totalFraction: 0, wantPoint: c, wantEdgeID: 0, wantDistance: 0, wantErr: true}, + } { + distancePoint, distanceEdgeID, newDistance, err := args.query.AtFraction(args.totalFraction) + if args.wantErr != (err != nil) { + t.Fatalf("got %v, want %v", err, args.wantErr) + } + if distancePoint.Angle(args.wantPoint.Vector) >= kEpsilonAngle { + t.Errorf("got %v, want %v", distancePoint, args.wantPoint.Vector) + } + if distanceEdgeID != args.wantEdgeID { + t.Errorf("got %v, want %v", distanceEdgeID, args.wantEdgeID) + } + if !float64Near(newDistance.Radians(), args.wantDistance.Radians(), kEpsilon) { + t.Errorf("got %v, want %v", newDistance, args.wantDistance) + } + } +} +func TestDistances(t *testing.T) { + // Initialize test data + distances := []float64{ + -1.0, -1.0e-8, 0.0, 1.0e-8, 0.2, 0.5, + 1.0 - 1.0e-8, 1.0, 1.0 + 1.e-8, 1.2, 1.2, 1.2 + 1.0e-10, + 1.5, 1.999999, 2.0, 2.00000001, 1.e6, + } + + vertices := parsePoints( + `0:0, 0:0, 1.0e-7:0, 0.1:0, 0.2:0, 0.2:0, 0.6:0, 0.999999:0, 0.999999:0, + 1:0, 1:0, 1.000001:0, 1.000001:0, 1.1:0, 1.2:0, 1.2000001:0, 1.7:0, + 1.99999999:0, 2:0`, + ) + + totalLength := vertices[0].Angle(vertices[len(vertices)-1].Vector).Degrees() + + shape := Polyline(vertices) + query := InitChainInterpolationQuery(&shape, 0) + + angle, err := query.GetLength() + + if err != nil { + t.Errorf("got %v, want %v", err, nil) + } + length := angle.Degrees() + + results := make([]result, len(distances)) + for i := 0; i < len(distances); i++ { + point, edgeID, distance, err := query.AtDistance(s1.Angle(distances[i] * float64(s1.Degree))) + + results[i] = result{point, edgeID, distance, err} + } + + if !float64Near(length, totalLength, kEpsilon) { + t.Errorf("got %v, want %v", length, totalLength) + } + + // Run tests + + for i := 0; i < len(distances); i++ { + if results[i].err != nil { + t.Errorf("got %v, want %v", results[i].err, nil) + } + + d := distances[i] + lat := LatLngFromPoint(results[i].point).Lat.Degrees() + edgeID := results[i].edgeID + distance := results[i].distance + + if d < 0 { + if !float64Eq(lat, LatLngFromPoint(shape.Edge(0).V0).Lat.Degrees()) { + t.Errorf("got %v, want %v", lat, 0) + } + + if edgeID != 0 { + t.Errorf("got %v, want %v", edgeID, 0) + } + + if !float64Eq(distance.Degrees(), 0) { + t.Errorf("got %v, want %v", distance, 0) + } + } else if d > 2 { + if !float64Near(lat, 2, kEpsilon) { + t.Errorf("got %v, want %v", lat, 2) + } + + if edgeID != shape.NumEdges()-1 { + t.Errorf("got %v, want %v", edgeID, shape.NumEdges()-1) + } + + if !float64Eq(distance.Degrees(), totalLength) { + t.Errorf("got %v, want %v", distance, totalLength) + } + } else { + if !float64Near(lat, d, kEpsilon) { + t.Errorf("got %v, want %v", lat, d) + } + + if edgeID < 0 { + t.Errorf("got %v, want %v", edgeID, 0) + } + + if edgeID >= shape.NumEdges() { + t.Errorf("got %v, want %v", edgeID, shape.NumEdges()-1) + } + + edge := shape.Edge(edgeID) + + if lat < LatLngFromPoint(edge.V0).Lat.Degrees() { + t.Errorf("got %v, want %v", lat, LatLngFromPoint(edge.V0).Lat.Degrees()) + } + + if lat > LatLngFromPoint(edge.V1).Lat.Degrees() { + t.Errorf("got %v, want %v", lat, LatLngFromPoint(edge.V1).Lat.Degrees()) + } + + if !float64Near(distance.Degrees(), d, kEpsilon) { + t.Errorf("got %v, want %v", distance, d) + } + } + } +} +func TestChains(t *testing.T) { + loops := [][]Point{ + parsePoints(`0:0, 1:0`), + parsePoints(`2:0, 3:0`), + } + + laxPolygon := LaxPolygonFromPoints(loops) + + tests := []struct { + query ChainInterpolationQuery + want result + args float64 + }{ + { + query: InitChainInterpolationQuery(laxPolygon, -1), + want: result{ + point: PointFromLatLng(LatLngFromDegrees(1, 0)), + edgeID: 1, + distance: s1.Angle(1 * s1.Degree), + err: nil, + }, + args: 0.25, + }, + { + query: InitChainInterpolationQuery(laxPolygon, 0), + want: result{ + point: PointFromLatLng(LatLngFromDegrees(0.5, 0)), + edgeID: 0, + distance: s1.Angle(0.5 * s1.Degree), + err: nil, + }, + args: 0.25, + }, + { + query: InitChainInterpolationQuery(laxPolygon, 1), + want: result{ + point: PointFromLatLng(LatLngFromDegrees(2.5, 0)), + edgeID: 2, + distance: s1.Angle(2.5 * s1.Degree), + err: nil, + }, + args: 0.25, + }, + } + + for i, tt := range tests { + point, edgeID, distance, err := tt.query.AtFraction(tt.args) + + got := result{ + point: point, + edgeID: edgeID, + distance: distance, + err: err, + } + + if !float64Near(LatLngFromPoint(got.point).Lat.Degrees(), LatLngFromPoint(tt.want.point).Lat.Degrees(), kEpsilon) { + t.Errorf("%d. got %v, want %v", i, LatLngFromPoint(got.point).Lat.Degrees(), LatLngFromPoint(tt.want.point).Lat.Degrees()) + } + + if got.edgeID != tt.want.edgeID { + t.Errorf("%d. got %v, want %v", i, got.edgeID, tt.want.edgeID) + } + if got.err != tt.want.err { + t.Errorf("%d. got %v, want %v", i, got.err, tt.want.err) + } + } +} + +func TestGetLengthAtEdgeEmpty(t *testing.T) { + query := InitChainInterpolationQuery(&laxPolyline{}, 0) + + angle, err := query.GetLengthAtEdgeEnd(0) + + if err == nil { + t.Errorf("got %v, want %v", err, nil) + } + + if !float64Eq(angle.Degrees(), 0) { + t.Errorf("got %v, want %v", angle, 0) + } +} +func TestGetLengthAtEdgePolyline(t *testing.T) { + points := []Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 3)), + PointFromLatLng(LatLngFromDegrees(0, 6)), + } + + polyline := laxPolyline{points} + + query := InitChainInterpolationQuery(&polyline, 0) + + tests := []struct { + edgeID int + want s1.Angle + }{ + {-100, s1.InfAngle()}, + {0, s1.Degree}, + {1, s1.Degree * 3}, + {2, s1.Degree * 6}, + {100, s1.InfAngle()}, + } + + for _, tt := range tests { + got, err := query.GetLengthAtEdgeEnd(tt.edgeID) + + if err != nil { + t.Errorf("edgeID %d: got %v, want %v", tt.edgeID, err, nil) + } + + if tt.edgeID <= polyline.NumEdges() && tt.edgeID >= 0 { + if !float64Near(got.Degrees(), tt.want.Degrees(), kEpsilon) { + t.Errorf("edgeID %d: got %v, want %v", tt.edgeID, got, tt.want) + } + } else if got != tt.want { + t.Errorf("edgeID %d: got %v, want %v", tt.edgeID, got, tt.want) + } + } +} +func TestGetLengthAtEdgePolygon(t *testing.T) { + polygon := laxPolygonFromPoints([][]Point{ + { + PointFromLatLng(LatLngFromDegrees(1, 1)), + PointFromLatLng(LatLngFromDegrees(2, 1)), + PointFromLatLng(LatLngFromDegrees(2, 3)), + PointFromLatLng(LatLngFromDegrees(1, 3)), + }, + { + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 4)), + PointFromLatLng(LatLngFromDegrees(3, 4)), + PointFromLatLng(LatLngFromDegrees(3, 0)), + }}) + + tolerance := .01 + + tests := []struct { + name string + args struct { + shape Shape + edge int + chainID int + } + want s1.Angle + }{ + { + name: "edge before first edge of first loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: -100, + chainID: 0, + }, + want: s1.InfAngle(), + }, + { + name: "first edge of first loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 0, + chainID: 0, + }, + want: s1.Degree, + }, + { + name: "second edge of first loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 1, + chainID: 0, + }, + want: s1.Degree * 3, + }, + { + name: "last edge of first loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 3, + chainID: 0, + }, + want: s1.Degree * 6, + }, + { + name: "edge after last edge of first loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 4, + chainID: 0, + }, + want: s1.InfAngle(), + }, + { + name: "edge before first edge of second loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 3, + chainID: 1, + }, + want: s1.InfAngle(), + }, + { + name: "first edge of second loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 4, + chainID: 1, + }, + want: s1.Degree * 4, + }, + { + name: "second edge of second loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 5, + chainID: 1, + }, + want: s1.Degree * 7, + }, + { + name: "midlle edge of second loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 6, + chainID: 1, + }, + want: s1.Degree * 11, + }, + { + name: "last edge of second loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 7, + chainID: 1, + }, + want: s1.Degree * 14, + }, + { + name: "edge after last edge of second loop", + args: struct { + shape Shape + edge int + chainID int + }{ + shape: polygon, + edge: 8, + chainID: 1, + }, + want: s1.InfAngle(), + }, + } + + for _, tt := range tests { + query := InitChainInterpolationQuery(tt.args.shape, tt.args.chainID) + + got, err := query.GetLengthAtEdgeEnd(tt.args.edge) + + if err != nil { + t.Errorf("%d. got %v, want %v", tt.args.edge, err, nil) + } + + if tt.want == s1.InfAngle() { + if got != tt.want { + t.Errorf("%d. got %v, want %v", tt.args.edge, got, tt.want) + } + } else if !float64Near(got.Degrees(), tt.want.Degrees(), float64(tolerance)) { + t.Errorf("%d. got %v, want %v", tt.args.edge, got.Degrees(), tt.want.Degrees()) + } + } +} + +func TestSlice(t *testing.T) { + tests := []struct { + name string + args struct { + shape Shape + startSliceFraction float64 + endSliceFraction float64 + } + want string + }{ + { + name: "empty shape", + args: struct { + shape Shape + startSliceFraction float64 + endSliceFraction float64 + }{nil, 0, 1}, + want: ``, + }, + { + name: "full polyline", + args: struct { + shape Shape + startSliceFraction float64 + endSliceFraction float64 + }{laxPolylineFromPoints(parsePoints(`0:0, 0:1, 0:2`)), 0, 1}, + want: `0:0, 0:1, 0:2`, + }, + { + name: "first half of polyline", + args: struct { + shape Shape + startSliceFraction float64 + endSliceFraction float64 + }{laxPolylineFromPoints(parsePoints(`0:0, 0:1, 0:2`)), 0, 0.5}, + want: `0:0, 0:1`, + }, + { + name: "second half of polyline", + args: struct { + shape Shape + startSliceFraction float64 + endSliceFraction float64 + }{laxPolylineFromPoints(parsePoints(`0:0, 0:1, 0:2`)), 1, 0.5}, + want: `0:2, 0:1`, + }, + { + name: "middle of polyline", + args: struct { + shape Shape + startSliceFraction float64 + endSliceFraction float64 + }{laxPolylineFromPoints(parsePoints(`0:0, 0:1, 0:2`)), 0.25, 0.75}, + want: `0:0.5, 0:1, 0:1.5`, + }, + } + + for _, test := range tests { + query := InitChainInterpolationQuery(test.args.shape, 0) + if got := pointsToString(query.Slice(test.args.startSliceFraction, test.args.endSliceFraction)); got != test.want { + t.Errorf("%v: got %v, want %v", test.name, got, test.want) + } + } +} + +func TestSliceDivided(t *testing.T) { + type args struct { + shape Shape + startSliceFraction float64 + endSliceFraction float64 + divisions int + } + tests := []struct { + name string + args args + want string + }{ + { + name: "empty shape", + args: args{nil, 0, 1., 1}, + want: ``, + }, + {name: "full polyline", args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0, + endSliceFraction: 1, + divisions: 3, + }, want: `0:0, 0:1, 0:2`}, + { + name: "first half of polyline", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0, + endSliceFraction: 0.5, + divisions: 2, + }, + want: `0:0, 0:1`, + }, + { + name: "second half of polyline", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 1, + endSliceFraction: 0.5, + divisions: 2, + }, + want: `0:2, 0:1`, + }, + { + name: "middle of polyline", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.25, + endSliceFraction: 0.75, + divisions: 3, + }, + want: `0:0.5, 0:1, 0:1.5`, + }, + { + name: "middle of polyline; divisions = 5", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.25, + endSliceFraction: 0.75, + divisions: 5, + }, + want: `0:0.5, 0:0.75, 0:1, 0:1.25, 0:1.5`, + }, + { + name: "middle of polyline; divisions = 11", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.25, + endSliceFraction: 0.75, + divisions: 11, + }, + want: `0:0.5, 0:0.6, 0:0.7, 0:0.8, 0:0.9, 0:1, 0:1.1, 0:1.2, 0:1.3, 0:1.4, 0:1.5`, + }, + { + name: "corner case: divisions = s.NumEdges()+1", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.3, + endSliceFraction: 0.6, + divisions: 4, + }, + want: `0:0.6, 0:0.8, 0:1, 0:1.2`, + }, + { + name: "divisions = s.NumEdges()+2", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.25, + endSliceFraction: 0.75, + divisions: 5, + }, + want: `0:0.5, 0:0.75, 0:1, 0:1.25, 0:1.5`, + }, + { + name: "divisions = s.NumEdges()+3", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.45, + endSliceFraction: 0.75, + divisions: 6, + }, + want: `0:0.9, 0:1, 0:1.14, 0:1.26, 0:1.38, 0:1.5`, + }, + { + name: "divisions = s.NumEdges()+10", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.105, + endSliceFraction: 0.605, + divisions: 11, + }, + want: `0:0.21, 0:0.31, 0:0.41, 0:0.51, 0:0.61, 0:0.71, 0:0.81, 0:0.91, 0:1, 0:1.11, 0:1.21`, + }, + { + name: "divisions = 10, 0 edges inside resulting points", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.05, + endSliceFraction: 0.1, + divisions: 11, + }, + want: `0:0.1, 0:0.11, 0:0.12, 0:0.13, 0:0.14, 0:0.15, 0:0.16, 0:0.17, 0:0.18, 0:0.19, 0:0.2`, + }, + { + name: "divisions = s.NumEdges()+1", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + PointFromLatLng(LatLngFromDegrees(0, 3)), + PointFromLatLng(LatLngFromDegrees(0, 4)), + PointFromLatLng(LatLngFromDegrees(0, 5)), + }, + ), + startSliceFraction: 0.3, + endSliceFraction: 0.84, + divisions: 5, + }, + want: `0:1.5, 0:2, 0:3, 0:4, 0:4.2`, + }, + { + name: "divisions = s.NumEdges()+1", + args: args{ + shape: laxPolylineFromPoints([]Point{ + PointFromLatLng(LatLngFromDegrees(0, 0)), + PointFromLatLng(LatLngFromDegrees(0, 1)), + PointFromLatLng(LatLngFromDegrees(0, 2)), + }, + ), + startSliceFraction: 0.3, + endSliceFraction: 0.99999999999995, + divisions: 3, + }, + want: `0:0.6, 0:1, 0:1.9999999999999`, + }, + } + + for _, test := range tests { + query := InitChainInterpolationQuery(test.args.shape, 0) + got := query.SliceDivided( + test.args.startSliceFraction, + test.args.endSliceFraction, + test.args.divisions, + ) + want := parsePoints(test.want) + if len(got) != test.args.divisions && len(got) != len(want) { + t.Errorf("length mismatch: got %d, want %d", len(got), test.args.divisions) + } + if !pointSlicesApproxEqual(got, want, kEpsilon) { + t.Errorf("%v: got %v, want %v", test.name, got, want) + } + } +} diff --git a/s2/edge_distances.go b/s2/edge_distances.go index 3067cd7..4ba0bfc 100644 --- a/s2/edge_distances.go +++ b/s2/edge_distances.go @@ -406,3 +406,14 @@ func EdgePairClosestPoints(a0, a1, b0, b1 Point) (Point, Point) { panic("illegal case reached") } } + +func GetPointOnRay(origin Point, dir Point, r s1.Angle) Point { + vector := (origin.Vector.Mul(math.Cos(r.Radians()))).Add(dir.Vector.Mul(math.Sin(r.Radians()))).Normalize() + return PointFromCoords(vector.X, vector.Y, vector.Z) +} + +func GetPointOnLine(startPoint Point, endPoint Point, angle s1.Angle) Point { + vector, _ := robustNormalWithLength(startPoint.Vector, endPoint.Vector) + dir := vector.Cross(startPoint.Vector).Normalize() + return GetPointOnRay(startPoint, PointFromCoords(dir.X, dir.Y, dir.Z), angle) +} diff --git a/s2/lax_loop.go b/s2/lax_loop.go index c544f03..1fb9b98 100644 --- a/s2/lax_loop.go +++ b/s2/lax_loop.go @@ -70,7 +70,13 @@ func (l *LaxLoop) Edge(e int) Edge { func (l *LaxLoop) Dimension() int { return 2 } func (l *LaxLoop) ReferencePoint() ReferencePoint { return referencePointForShape(l) } func (l *LaxLoop) NumChains() int { return minInt(1, l.numVertices) } -func (l *LaxLoop) Chain(i int) Chain { return Chain{0, l.numVertices} } +func (l *LaxLoop) Chain(i int) Chain { + if l.numVertices == 1 { + return Chain{0, l.numVertices} + } + return Chain{i, l.numVertices - i} +} + func (l *LaxLoop) ChainEdge(i, j int) Edge { var k int if j+1 == l.numVertices { diff --git a/s2/lax_loop_test.go b/s2/lax_loop_test.go index 9e3e5be..01578bb 100644 --- a/s2/lax_loop_test.go +++ b/s2/lax_loop_test.go @@ -68,7 +68,13 @@ func (l *laxLoop) Edge(e int) Edge { func (l *laxLoop) Dimension() int { return 2 } func (l *laxLoop) ReferencePoint() ReferencePoint { return referencePointForShape(l) } func (l *laxLoop) NumChains() int { return minInt(1, l.numVertices) } -func (l *laxLoop) Chain(i int) Chain { return Chain{0, l.numVertices} } +func (l *laxLoop) Chain(i int) Chain { + if l.numVertices == 1 { + return Chain{0, l.numVertices} + } + return Chain{i, l.numVertices - i} +} + func (l *laxLoop) ChainEdge(i, j int) Edge { var k int if j+1 == l.numVertices { diff --git a/s2/lax_polyline.go b/s2/lax_polyline.go index 87ad967..c5c731b 100644 --- a/s2/lax_polyline.go +++ b/s2/lax_polyline.go @@ -40,11 +40,16 @@ func LaxPolylineFromPolyline(p Polyline) *LaxPolyline { return LaxPolylineFromPoints(p) } -func (l *LaxPolyline) NumEdges() int { return maxInt(0, len(l.vertices)-1) } -func (l *LaxPolyline) Edge(e int) Edge { return Edge{l.vertices[e], l.vertices[e+1]} } -func (l *LaxPolyline) ReferencePoint() ReferencePoint { return OriginReferencePoint(false) } -func (l *LaxPolyline) NumChains() int { return minInt(1, l.NumEdges()) } -func (l *LaxPolyline) Chain(i int) Chain { return Chain{0, l.NumEdges()} } +func (l *LaxPolyline) NumEdges() int { return maxInt(0, len(l.vertices)-1) } +func (l *LaxPolyline) Edge(e int) Edge { return Edge{l.vertices[e], l.vertices[e+1]} } +func (l *LaxPolyline) ReferencePoint() ReferencePoint { return OriginReferencePoint(false) } +func (l *LaxPolyline) NumChains() int { return minInt(1, l.NumEdges()) } +func (l *LaxPolyline) Chain(i int) Chain { + if l.NumEdges() == 1 { + return Chain{0, l.NumEdges()} + } + return Chain{i, l.NumEdges() - i} +} func (l *LaxPolyline) ChainEdge(i, j int) Edge { return Edge{l.vertices[j], l.vertices[j+1]} } func (l *LaxPolyline) ChainPosition(e int) ChainPosition { return ChainPosition{0, e} } func (l *LaxPolyline) Dimension() int { return 1 } diff --git a/s2/lax_polyline_test.go b/s2/lax_polyline_test.go index 90d8a77..144f1f1 100644 --- a/s2/lax_polyline_test.go +++ b/s2/lax_polyline_test.go @@ -36,11 +36,16 @@ func laxPolylineFromPolyline(p Polyline) *laxPolyline { return laxPolylineFromPoints(p) } -func (l *laxPolyline) NumEdges() int { return maxInt(0, len(l.vertices)-1) } -func (l *laxPolyline) Edge(e int) Edge { return Edge{l.vertices[e], l.vertices[e+1]} } -func (l *laxPolyline) ReferencePoint() ReferencePoint { return OriginReferencePoint(false) } -func (l *laxPolyline) NumChains() int { return minInt(1, l.NumEdges()) } -func (l *laxPolyline) Chain(i int) Chain { return Chain{0, l.NumEdges()} } +func (l *laxPolyline) NumEdges() int { return maxInt(0, len(l.vertices)-1) } +func (l *laxPolyline) Edge(e int) Edge { return Edge{l.vertices[e], l.vertices[e+1]} } +func (l *laxPolyline) ReferencePoint() ReferencePoint { return OriginReferencePoint(false) } +func (l *laxPolyline) NumChains() int { return minInt(1, l.NumEdges()) } +func (l *laxPolyline) Chain(i int) Chain { + if l.NumEdges() == 1 { + return Chain{0, l.NumEdges()} + } + return Chain{i, l.NumEdges() - i} +} func (l *laxPolyline) ChainEdge(i, j int) Edge { return Edge{l.vertices[j], l.vertices[j+1]} } func (l *laxPolyline) ChainPosition(e int) ChainPosition { return ChainPosition{0, e} } func (l *laxPolyline) Dimension() int { return 1 } diff --git a/s2/loop.go b/s2/loop.go index 73cb67c..88f7b0f 100644 --- a/s2/loop.go +++ b/s2/loop.go @@ -501,7 +501,10 @@ func (l *Loop) NumChains() int { // Chain returns the i-th edge chain in the Shape. func (l *Loop) Chain(chainID int) Chain { - return Chain{0, l.NumEdges()} + if l.NumEdges() == 1 { + return Chain{0, l.NumEdges()} + } + return Chain{chainID, l.NumEdges() - chainID} } // ChainEdge returns the j-th edge of the i-th edge chain. diff --git a/s2/point_vector.go b/s2/point_vector.go index f8e6f65..b591b75 100644 --- a/s2/point_vector.go +++ b/s2/point_vector.go @@ -28,11 +28,16 @@ var ( // Its methods are on *PointVector due to implementation details of ShapeIndex. type PointVector []Point -func (p *PointVector) NumEdges() int { return len(*p) } -func (p *PointVector) Edge(i int) Edge { return Edge{(*p)[i], (*p)[i]} } -func (p *PointVector) ReferencePoint() ReferencePoint { return OriginReferencePoint(false) } -func (p *PointVector) NumChains() int { return len(*p) } -func (p *PointVector) Chain(i int) Chain { return Chain{i, 1} } +func (p *PointVector) NumEdges() int { return len(*p) } +func (p *PointVector) Edge(i int) Edge { return Edge{(*p)[i], (*p)[i]} } +func (p *PointVector) ReferencePoint() ReferencePoint { return OriginReferencePoint(false) } +func (p *PointVector) NumChains() int { return len(*p) } +func (p *PointVector) Chain(i int) Chain { + if len(*p) == 1 { + return Chain{0, 1} + } + return Chain{i, 1} +} func (p *PointVector) ChainEdge(i, j int) Edge { return Edge{(*p)[i], (*p)[j]} } func (p *PointVector) ChainPosition(e int) ChainPosition { return ChainPosition{e, 0} } func (p *PointVector) Dimension() int { return 0 } diff --git a/s2/point_vector_test.go b/s2/point_vector_test.go index 3fef472..8fd6b8e 100644 --- a/s2/point_vector_test.go +++ b/s2/point_vector_test.go @@ -72,7 +72,7 @@ func TestPointVectorBasics(t *testing.T) { rand.Seed(seed) for i := 0; i < numPoints; i++ { - if got, want := shape.Chain(i).Start, i; got != want { + if got, want := shape.Chain(i).Start, i; i != got { t.Errorf("shape.Chain(%d).Start = %d, want %d", i, got, want) } if got, want := shape.Chain(i).Length, 1; got != want { diff --git a/s2/polyline.go b/s2/polyline.go index a5d9a64..b1ec2e7 100644 --- a/s2/polyline.go +++ b/s2/polyline.go @@ -195,7 +195,10 @@ func (p *Polyline) NumChains() int { // Chain returns the i-th edge Chain in the Shape. func (p *Polyline) Chain(chainID int) Chain { - return Chain{0, p.NumEdges()} + if p.NumEdges() == 1 { + return Chain{0, p.NumEdges()} + } + return Chain{chainID, p.NumEdges() - chainID} } // ChainEdge returns the j-th edge of the i-th edge Chain.