diff --git a/go.mod b/go.mod index 044ca35..910987c 100644 --- a/go.mod +++ b/go.mod @@ -1,6 +1,6 @@ module github.com/pdok/texel -go 1.21 +go 1.21.6 require ( github.com/cpuguy83/go-md2man/v2 v2.0.3 // indirect @@ -16,6 +16,7 @@ require ( github.com/muesli/reflow v0.3.0 github.com/perimeterx/marshmallow v1.1.5 github.com/stretchr/testify v1.8.4 + github.com/tobshub/go-sortedmap v1.0.3 github.com/wk8/go-ordered-map/v2 v2.1.8 golang.org/x/exp v0.0.0-20231110203233-9a3e6036ecaa ) @@ -44,6 +45,5 @@ require ( github.com/gdey/errors v0.0.0-20190426172550-8ebd5bc891fb // indirect github.com/pmezard/go-difflib v1.0.0 // indirect github.com/russross/blackfriday/v2 v2.1.0 // indirect - github.com/umpc/go-sortedmap v0.0.0-20180422175548-64ab94c482f4 github.com/xrash/smetrics v0.0.0-20201216005158-039620a65673 // indirect ) diff --git a/go.sum b/go.sum index 331a59e..9db1eb2 100644 --- a/go.sum +++ b/go.sum @@ -67,10 +67,10 @@ github.com/stretchr/testify v1.8.0/go.mod h1:yNjHg4UonilssWZ8iaSj1OCr/vHnekPRkoO github.com/stretchr/testify v1.8.2/go.mod h1:w2LPCIKwWwSfY2zedu0+kehJoqGctiVI29o6fzry7u4= github.com/stretchr/testify v1.8.4 h1:CcVxjf3Q8PM0mHUKJCdn+eZZtm5yQwehR5yeSVQQcUk= github.com/stretchr/testify v1.8.4/go.mod h1:sz/lmYIOXD/1dqDmKjjqLyZ2RngseejIcXlSw2iwfAo= +github.com/tobshub/go-sortedmap v1.0.3 h1:oUhj/5tqzjTX4bhWqB1ZFTDtMULJ1ZYUnS8WAugSfjY= +github.com/tobshub/go-sortedmap v1.0.3/go.mod h1:JLxyU94+lfKuCgelxXpwRr29ei6SqLbaeVuNVMvENbE= github.com/ugorji/go/codec v1.2.7 h1:YPXUKf7fYbp/y8xloBqZOw2qaVggbfwMlI8WM3wZUJ0= github.com/ugorji/go/codec v1.2.7/go.mod h1:WGN1fab3R1fzQlVQTkfxVtIBhWDRqOviHU95kRgeqEY= -github.com/umpc/go-sortedmap v0.0.0-20180422175548-64ab94c482f4 h1:qk1XyC6UGfPa51PGmsTQJavyhfMLScqw97pEV3sFClI= -github.com/umpc/go-sortedmap v0.0.0-20180422175548-64ab94c482f4/go.mod h1:X6iKjXCleSyo/LZzKZ9zDF/ZB2L9gC36I5gLMf32w3M= github.com/urfave/cli/v2 v2.25.7 h1:VAzn5oq403l5pHjc4OhD54+XGO9cdKVL/7lDjF+iKUs= github.com/urfave/cli/v2 v2.25.7/go.mod h1:8qnjx1vcq5s2/wpsqoZFndg2CE5tNFyrTvS6SinrnYQ= github.com/wk8/go-ordered-map/v2 v2.1.8 h1:5h/BUHu93oj4gIdvHHHGsScSTMijfx5PeYkE/fJgbpc= diff --git a/snap/geom.go b/snap/geom.go new file mode 100644 index 0000000..383b1ef --- /dev/null +++ b/snap/geom.go @@ -0,0 +1,18 @@ +package snap + +import "math" + +// https://en.wikipedia.org/wiki/Shoelace_formula +func shoelace(pts [][2]float64) float64 { + sum := 0. + if len(pts) == 0 { + return 0. + } + + p0 := pts[len(pts)-1] + for _, p1 := range pts { + sum += p0[1]*p1[0] - p0[0]*p1[1] + p0 = p1 + } + return math.Abs(sum / 2) +} diff --git a/snap/pointindex_test.go b/snap/pointindex_test.go index 1f419ee..e82adb5 100644 --- a/snap/pointindex_test.go +++ b/snap/pointindex_test.go @@ -450,31 +450,6 @@ func newSimplePointIndex(maxDepth Level, cellSize float64) *PointIndex { return &ix } -//nolint:unparam -func newSimpleTileMatrixSet(maxDepth Level, cellSize float64) tms20.TileMatrixSet { - zeroZero := tms20.TwoDPoint([2]float64{0.0, 0.0}) - tms := tms20.TileMatrixSet{ - CRS: fakeCRS{}, - OrderedAxes: []string{"X", "Y"}, - TileMatrices: make(map[tms20.TMID]tms20.TileMatrix), - } - for tmID := 0; tmID <= int(maxDepth); tmID++ { - tmCellSize := cellSize * float64(pow2(maxDepth-uint(tmID))) - tms.TileMatrices[tmID] = tms20.TileMatrix{ - ID: "0", - ScaleDenominator: tmCellSize / tms20.StandardizedRenderingPixelSize, - CellSize: tmCellSize, - CornerOfOrigin: tms20.BottomLeft, - PointOfOrigin: &zeroZero, - TileWidth: 1, - TileHeight: 1, - MatrixWidth: 1, - MatrixHeight: 1, - } - } - return tms -} - func loadEmbeddedTileMatrixSet(t *testing.T, tmsID string) tms20.TileMatrixSet { tms, err := tms20.LoadEmbeddedTileMatrixSet(tmsID) require.NoError(t, err) diff --git a/snap/snap.go b/snap/snap.go index 4b39b0e..a9c7a47 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -2,17 +2,20 @@ package snap import ( "fmt" + "github.com/tobshub/go-sortedmap" + "log" "math" "slices" "sort" + "github.com/go-spatial/geom/winding" + "github.com/go-spatial/geom/encoding/wkt" "github.com/muesli/reflow/truncate" "github.com/go-spatial/geom" "github.com/pdok/texel/intgeom" "github.com/pdok/texel/tms20" - "github.com/umpc/go-sortedmap" orderedmap "github.com/wk8/go-ordered-map/v2" "golang.org/x/exp/constraints" "golang.org/x/exp/maps" @@ -23,6 +26,8 @@ const ( internalPixelResolution = 16 ) +type IsOuter = bool + // SnapPolygon snaps polygons' points to a tile's internal pixel grid // and adds points to lines to prevent intersections. // @@ -89,7 +94,8 @@ func tileMatrixIDsByLevels(tms tms20.TileMatrixSet, tmIDs []tms20.TMID) map[Leve //nolint:cyclop func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[Level][]geom.Polygon { levelMap := asKeys(levels) - newPolygons := make(map[Level][][][][2]float64, len(levels)) + newOuters := make(map[Level][][][2]float64, len(levels)) + newInners := make(map[Level][][][2]float64, len(levels)) newPointsAndLines := make(map[Level][][][2]float64, len(levels)) // Could use polygon.AsSegments(), but it skips rings with <3 segments and starts with the last segment. @@ -99,7 +105,7 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ } isOuter := ringIdx == 0 // winding order is reversed if incorrect - ensureCorrectWindingOrder(ring, !isOuter) + ring = ensureCorrectWindingOrder(ring, !isOuter) ringLen := len(ring) newRing := make(map[Level][][2]float64, len(levelMap)) for level := range levelMap { @@ -122,23 +128,30 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ // walk through the new ring and append to the polygon (on all levels) for level := range levelMap { outerRings, innerRings, pointsAndLines := cleanupNewRing(newRing[level], isOuter, ix.hitMultiple[level], ringIdx) + // Check if outer ring has become too small if isOuter && len(outerRings) == 0 && (!keepPointsAndLines || len(pointsAndLines) == 0) { - delete(levelMap, level) // outer ring has become too small + delete(levelMap, level) // If too small, delete it continue } - for _, outerRing := range outerRings { // (there are only outer rings if isOuter) - newPolygons[level] = append(newPolygons[level], [][][2]float64{outerRing}) + for _, outerRing := range outerRings { + newOuters[level] = append(newOuters[level], outerRing) } - if len(newPolygons[level]) == 0 && len(innerRings) > 0 { // should never happen - panicInnerRingsButNoOuterRings(level, polygon, innerRings) - continue - } - newPolygons[level] = matchInnersToPolygon(newPolygons[level], innerRings) + newInners[level] = append(newInners[level], innerRings...) if keepPointsAndLines { newPointsAndLines[level] = append(newPointsAndLines[level], pointsAndLines...) } } } + + newPolygons := make(map[Level][][][][2]float64, len(levels)) + for l := range levelMap { + newOuters[l], newInners[l] = dedupeInnersOuters(newOuters[l], newInners[l]) + newPolygonsForLevel := matchInnersToPolygons(outersToPolygons(newOuters[l]), newInners[l], len(polygon) > 1) + if len(newPolygonsForLevel) > 0 { + newPolygons[l] = newPolygonsForLevel + } + } + // points and lines at the end, as outer rings for level, pointsAndLines := range newPointsAndLines { for _, pointOrLine := range pointsAndLines { @@ -148,46 +161,174 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ return floatPolygonsToGeomPolygonsForAllLevels(newPolygons) } -func matchInnersToPolygon(polygons [][][][2]float64, innerRings [][][2]float64) [][][][2]float64 { +func outersToPolygons(outers [][][2]float64) [][][][2]float64 { + polygons := make([][][][2]float64, len(outers)) + for i := 0; i < len(outers); i++ { + polygons[i] = [][][2]float64{outers[i]} + } + return polygons +} + +func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]float64, [][][2]float64) { + // ToDo: optimize by deleting rings from allRings on the fly + lenOuters := len(outers) + lenInners := len(inners) + lenAll := lenOuters + lenInners + processedIndexes := make(map[int]IsOuter) + indexesToDelete := make(map[int]IsOuter) + for i := 0; i < lenAll; i++ { + if _, processed := processedIndexes[i]; processed { + continue + } + iIsOuter := i < lenOuters + equalIndexes := orderedmap.New[int, IsOuter]( // ordered for predictable outcome when ranging + orderedmap.WithInitialData(orderedmap.Pair[int, IsOuter]{Key: i, Value: iIsOuter})) + for j := i + 1; j < lenAll; j++ { + if _, processed := processedIndexes[j]; processed { + continue + } + jIsOuter := j < lenOuters + var ringI [][2]float64 + if iIsOuter { + ringI = outers[i] + } else { + ringI = inners[i-lenOuters] + } + var ringJ [][2]float64 + if jIsOuter { + ringJ = outers[j] + } else { + ringJ = inners[j-lenOuters] + } + if !ringsAreEqual(ringI, ringJ, iIsOuter, jIsOuter) { + continue + } + equalIndexes.Set(j, jIsOuter) + } + if equalIndexes.Len() <= 1 { + continue + } + + lenEqualOuters := countVals(equalIndexes, true) + lenEqualInners := countVals(equalIndexes, false) + difference := int(math.Abs(float64(lenEqualOuters) - float64(lenEqualInners))) + var numOutersToDelete, numInnersToDelete int + if difference == 0 { + // delete all but one + numOutersToDelete = lenEqualOuters - 1 + numInnersToDelete = lenEqualInners - 1 + } else { // difference > 0 + // delete the surplus + numOutersToDelete = min(lenEqualOuters, lenEqualInners) + numInnersToDelete = numOutersToDelete + } + for p := equalIndexes.Oldest(); p != nil; p = p.Next() { + equalI, isOuter := p.Key, p.Value + processedIndexes[equalI] = isOuter + if isOuter && numOutersToDelete > 0 { + indexesToDelete[equalI] = isOuter + numOutersToDelete-- + } else if !isOuter && numInnersToDelete > 0 { + indexesToDelete[equalI] = isOuter + numInnersToDelete-- + } + } + } + + if len(indexesToDelete) == 0 { + return outers, inners + } + newOuters := deleteFromSliceByIndex(outers, indexesToDelete, 0) + newInners := deleteFromSliceByIndex(inners, indexesToDelete, lenOuters) + return newOuters, newInners +} + +func ringsAreEqual(ringI, ringJ [][2]float64, iIsOuter, jIsOuter bool) bool { + // check length + ringLen := len(ringI) + if ringLen != len(ringJ) { + return false + } + idx := slices.Index(ringJ, ringI[0]) // empty ring panics, but that's ok + if idx < 0 { + return false + } + differentWindingOrder := iIsOuter && !jIsOuter + // Check if rings are equal + for k := 0; k < ringLen; k++ { + if !differentWindingOrder && ringI[k] != ringJ[(idx+k)%ringLen] { + return false + } + if differentWindingOrder && ringI[k] != ringJ[(idx+ringLen-k)%ringLen] { // "+iLen" to ensure the modulus returns a positive number + return false + } + } + return true +} + +func matchInnersToPolygons(polygons [][][][2]float64, innerRings [][][2]float64, hasInners bool) [][][][2]float64 { lenPolygons := len(polygons) - if lenPolygons == 0 { - return polygons - } else if lenPolygons == 1 { - polygons[0] = append(polygons[0], innerRings...) + if len(innerRings) == 0 { return polygons } + + var polyISortedByOuterAreaDesc []int + var innersTurnedOuters [][][2]float64 matchInners: for _, innerRing := range innerRings { - containsPerPolyI := make(map[int]uint, lenPolygons) + containsPerPolyI := orderedmap.New[int, uint](orderedmap.WithCapacity[int, uint](lenPolygons)) // TODO don't need ordered map anymore? // this is pretty nested, but usually breaks early for _, vertex := range innerRing { for polyI := range polygons { - contains, onBoundary := ringContains(polygons[polyI][0], vertex) + contains, _ := ringContains(polygons[polyI][0], vertex) + // it doesn't matter if on boundary or not, if not on boundary there could still be multiple (nested) matching polygons if contains { - if !onBoundary { - polygons[polyI] = append(polygons[polyI], innerRing) - continue matchInners - } - containsPerPolyI[polyI]++ + containsPerPolyI.Set(polyI, containsPerPolyI.Value(polyI)+1) } } - matchingPolyI, _, matchCount := findKeyWithMaxValue(containsPerPolyI) + matchingPolyI, _, matchCount := findLastKeyWithMaxValue(containsPerPolyI) if matchCount == 1 { polygons[matchingPolyI] = append(polygons[matchingPolyI], innerRing) continue matchInners } } - // no (single) matching outer ring was found (should never happen) - if len(containsPerPolyI) == 0 { - panicNoMatchingOuterForInnerRing(polygons, innerRing) + if containsPerPolyI.Len() == 0 { + // no (single) matching outer ring was found + // presumably because the inner ring's winding order is incorrect and it should have been an outer + // TODO is that presumption correct and is this really never a panic? // panicNoMatchingOuterForInnerRing(polygons, innerRing) + // TODO should it be a candidate for other the other inner rings? + log.Printf("no matching outer for inner ring found, turned inner into outer. original has inners: %v", hasInners) + innersTurnedOuters = append(innersTurnedOuters, reverseClone(innerRing)) continue } - panicMoreThanOneMatchingOuterRing(polygons, innerRing) - continue + // multiple matching outer rings were found. use the smallest one + // TODO dedupe poly outers (not here) + if polyISortedByOuterAreaDesc == nil { + polyISortedByOuterAreaDesc = sortPolyIdxsByOuterAreaDesc(polygons) + } + smallestMatchingPolyI := lastMatch(polyISortedByOuterAreaDesc, orderedMapKeys(containsPerPolyI)) + polygons[smallestMatchingPolyI] = append(polygons[smallestMatchingPolyI], innerRing) + } + for i := range innersTurnedOuters { + polygons = append(polygons, [][][2]float64{innersTurnedOuters[i]}) } return polygons } +func sortPolyIdxsByOuterAreaDesc(polygons [][][][2]float64) []int { + areas := sortedmap.New[int, float64](len(polygons), func(i, j float64) bool { + return i > j // desc + }) + for i := range polygons { + if len(polygons[i]) == 0 { + areas.Insert(i, 0.0) + } else { + areas.Insert(i, shoelace(polygons[i][0])) + } + } + return areas.Keys() +} + // from paulmach/orb, modified to also return whether it's on the boundary // ringContains returns true if the point is inside the ring. // Points on the boundary are also considered in. In which case the second returned var is true too. @@ -217,58 +358,58 @@ func ringContains(ring [][2]float64, point [2]float64) (contains, onBoundary boo // Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go // //nolint:cyclop,nestif -func rayIntersect(p, s, e [2]float64) (intersects, on bool) { - if s[0] > e[0] { - s, e = e, s +func rayIntersect(pt, start, end [2]float64) (intersects, on bool) { + if start[xAx] > end[xAx] { + start, end = end, start } - if p[0] == s[0] { - if p[1] == s[1] { - // p == start + if pt[xAx] == start[xAx] { + if pt[yAx] == start[yAx] { + // pt == start return false, true - } else if s[0] == e[0] { - // vertical segment (s -> e) + } else if start[xAx] == end[xAx] { + // vertical segment (start -> end) // return true if within the line, check to see if start or end is greater. - if s[1] > e[1] && s[1] >= p[1] && p[1] >= e[1] { + if start[yAx] > end[yAx] && start[yAx] >= pt[yAx] && pt[yAx] >= end[yAx] { return false, true } - if e[1] > s[1] && e[1] >= p[1] && p[1] >= s[1] { + if end[yAx] > start[yAx] && end[yAx] >= pt[yAx] && pt[yAx] >= start[yAx] { return false, true } } // Move the y coordinate to deal with degenerate case - p[0] = math.Nextafter(p[0], math.Inf(1)) - } else if p[0] == e[0] { - if p[1] == e[1] { + pt[xAx] = math.Nextafter(pt[xAx], math.Inf(1)) + } else if pt[xAx] == end[xAx] { + if pt[yAx] == end[yAx] { // matching the end point return false, true } - p[0] = math.Nextafter(p[0], math.Inf(1)) + pt[xAx] = math.Nextafter(pt[xAx], math.Inf(1)) } - if p[0] < s[0] || p[0] > e[0] { + if pt[xAx] < start[xAx] || pt[xAx] > end[xAx] { return false, false } - if s[1] > e[1] { - if p[1] > s[1] { + if start[yAx] > end[yAx] { + if pt[yAx] > start[yAx] { return false, false - } else if p[1] < e[1] { + } else if pt[yAx] < end[yAx] { return true, false } } else { - if p[1] > e[1] { + if pt[yAx] > end[yAx] { return false, false - } else if p[1] < s[1] { + } else if pt[yAx] < start[yAx] { return true, false } } - rs := (p[1] - s[1]) / (p[0] - s[0]) - ds := (e[1] - s[1]) / (e[0] - s[0]) + rs := (pt[yAx] - start[yAx]) / (pt[xAx] - start[xAx]) + ds := (end[yAx] - start[yAx]) / (end[xAx] - start[xAx]) if rs == ds { return false, true @@ -318,32 +459,27 @@ func cleanupNewRing(newRing [][2]float64, isOuter bool, hitMultiple map[intgeom. } // if winding order is incorrect, ring is reversed to correct winding order -func ensureCorrectWindingOrder(ring [][2]float64, shouldBeClockwise bool) { +func ensureCorrectWindingOrder(ring [][2]float64, shouldBeClockwise bool) [][2]float64 { if !windingOrderIsCorrect(ring, shouldBeClockwise) { - slices.Reverse(ring) + return reverseClone(ring) } + return ring } // validate winding order (CCW for outer rings, CW for inner rings) func windingOrderIsCorrect(ring [][2]float64, shouldBeClockwise bool) bool { - // modulo function that returns least positive remainder (i.e., mod(-1, 5) returns 4) - mod := func(a, b int) int { - return (a%b + b) % b - } - i, _ := findRightmostLowestPoint(ring) - // check angle between the vectors goint into and coming out of the rightmost lowest point - points := [3][2]float64{ring[mod(i-1, len(ring))], ring[i], ring[mod(i+1, len(ring))]} - return isClockwise(points) == shouldBeClockwise + wo := winding.Order{}.OfPoints(ring...) + return wo.IsClockwise() && shouldBeClockwise || wo.IsCounterClockwise() && !shouldBeClockwise || wo.IsColinear() } // TODO: rewrite by using intgeoms for as long as possible func isHitMultiple(hitMultiple map[intgeom.Point][]int, vertex [2]float64, ringIdx int) bool { intVertex := intgeom.FromGeomPoint(vertex) return slices.Contains(hitMultiple[intVertex], ringIdx) || // exact match - slices.Contains(hitMultiple[intgeom.Point{intVertex[0] + 1, intVertex[1]}], ringIdx) || // fuzzy search - slices.Contains(hitMultiple[intgeom.Point{intVertex[0] - 1, intVertex[1]}], ringIdx) || - slices.Contains(hitMultiple[intgeom.Point{intVertex[0], intVertex[1] + 1}], ringIdx) || - slices.Contains(hitMultiple[intgeom.Point{intVertex[0], intVertex[1] - 1}], ringIdx) + slices.Contains(hitMultiple[intgeom.Point{intVertex[xAx] + 1, intVertex[yAx]}], ringIdx) || // fuzzy search + slices.Contains(hitMultiple[intgeom.Point{intVertex[xAx] - 1, intVertex[yAx]}], ringIdx) || + slices.Contains(hitMultiple[intgeom.Point{intVertex[xAx], intVertex[yAx] + 1}], ringIdx) || + slices.Contains(hitMultiple[intgeom.Point{intVertex[xAx], intVertex[yAx] - 1}], ringIdx) } // split ring into multiple rings at any point where the ring goes through the point more than once @@ -428,38 +564,36 @@ func splitRing(ring [][2]float64, isOuter bool, hitMultiple map[intgeom.Point][] } } } - return outerRings, innerRings, pointsAndLines -} - -// determines whether a pair of vectors turns clockwise by examining the relationship of their relative angle to the angles of the vectors -func isClockwise(points [3][2]float64) bool { - // modulo function that returns least positive remainder (i.e., mod(-1, 5) returns 4) - mod := func(a, b float64) float64 { - return math.Mod(math.Mod(a, b)+b, b) - } - vector1 := vector2d{x: (points[1][0] - points[0][0]), y: (points[1][1] - points[0][1])} - vector2 := vector2d{x: (points[2][0] - points[1][0]), y: (points[2][1] - points[1][1])} - relativeAngle := math.Acos(vector1.dot(vector2)/(vector1.magnitude()*vector2.magnitude())) * (180 / math.Pi) - if vector1.angle() == 0.0 { - return relativeAngle > 180 + // outer ring(s) incorrectly saved as inner ring(s) or vice versa due to winding order, swap + if isOuter && len(outerRings) == 0 && len(innerRings) > 0 { + for _, innerRing := range innerRings { + slices.Reverse(innerRing) // in place, not used elsewhere + outerRings = append(outerRings, innerRing) + } + innerRings = make([][][2]float64, 0) + } else if !isOuter && len(innerRings) == 0 && len(outerRings) > 0 { + for _, outerRing := range outerRings { + slices.Reverse(outerRing) // in place, not used elsewhere + innerRings = append(innerRings, outerRing) + } + outerRings = make([][][2]float64, 0) } - return math.Round(mod((vector2.angle()-relativeAngle), 360)) != math.Round(vector1.angle()) + return outerRings, innerRings, pointsAndLines } // deduplication using an implementation of the Knuth-Morris-Pratt algorithm // //nolint:cyclop,funlen -func kmpDeduplicate(newRing [][2]float64) [][2]float64 { - deduplicatedRing := make([][2]float64, len(newRing)) - copy(deduplicatedRing, newRing) - // map of indices to remove, sorted by starting index of each sequence to remove - indicesToRemove := sortedmap.New(len(newRing), func(x, y interface{}) bool { - return x.([2]int)[0] < y.([2]int)[0] +func kmpDeduplicate(ring [][2]float64) [][2]float64 { + ringLen := len(ring) + // sequences (from index uptoandincluding index) to remove, sorted by starting index, mapped to prevent dupes + sequencesToRemove := sortedmap.New[string, [2]int](ringLen, func(a, b [2]int) bool { + return a[xAx] < b[xAx] }) - // walk through newRing until a step back is taken, then identify how many steps back are taken and search for repeats + // walk through ring until a step back is taken, then identify how many steps back are taken and search for repeats visitedPoints := [][2]float64{} - for i := 0; i < len(newRing); { - vertex := newRing[i] + for i := 0; i < ringLen; { + vertex := ring[i] // not a step back, continue if len(visitedPoints) <= 1 || visitedPoints[len(visitedPoints)-2] != vertex { visitedPoints = append(visitedPoints, vertex) @@ -470,7 +604,7 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { reverseSegment := [][2]float64{visitedPoints[len(visitedPoints)-1], visitedPoints[len(visitedPoints)-2]} for j := 3; j <= len(visitedPoints); j++ { nextI := i + (j - 2) - if nextI <= len(newRing)-1 && visitedPoints[len(visitedPoints)-j] == newRing[nextI] { + if nextI <= ringLen-1 && visitedPoints[len(visitedPoints)-j] == ring[nextI] { reverseSegment = append(reverseSegment, visitedPoints[len(visitedPoints)-j]) } else { // end of segment @@ -480,13 +614,13 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { // create segment from reverse segment segment := make([][2]float64, len(reverseSegment)) copy(segment, reverseSegment) - slices.Reverse(segment) + slices.Reverse(segment) // in place, not used elsewhere // create search corpus: initialise with section of 3*segment length, then continuously // add 2*segment length until corpus contains a point that is not in segment start := i - len(segment) end := start + (3 * len(segment)) k := 0 - corpus := newRing[start:min(end, len(newRing))] + corpus := ring[start:min(end, ringLen)] for { stop := false // check if (additional) corpus contains a point that is not in segment @@ -496,8 +630,8 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { break } } - // corpus already runs until the end of newRing - if end > len(newRing) { + // corpus already runs until the end of ring + if end > ringLen { stop = true } if stop { @@ -505,7 +639,7 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { } // expand corpus k = len(corpus) - corpus = append(corpus, newRing[end:min(end+(2*len(segment)), len(newRing))]...) + corpus = append(corpus, ring[end:min(end+(2*len(segment)), ringLen)]...) end += 2 * len(segment) } // search corpus for all matches of segment and reverseSegment @@ -517,11 +651,7 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { // mark all but one occurrance of segment for removal sequenceStart := start + len(segment) sequenceEnd := start + matches[len(matches)-1] + len(segment) - segmentRec := sortedmap.Record{ - Key: fmt.Sprintf("%v", segment), - Val: [2]int{sequenceStart, sequenceEnd}, - } - indicesToRemove.Insert(segmentRec.Key, segmentRec.Val) + sequencesToRemove.Insert(fmt.Sprint(segment), [2]int{sequenceStart, sequenceEnd}) // skip past matched section and reset visitedPoints i = sequenceEnd visitedPoints = [][2]float64{} @@ -530,11 +660,7 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { // mark all but one occurrance of segment and one occurrance of its reverse for removal sequenceStart := start + (2 * len(segment)) - 1 sequenceEnd := start + matches[len(matches)-1] + len(segment) - segmentRec := sortedmap.Record{ - Key: fmt.Sprintf("%v", segment), - Val: [2]int{sequenceStart, sequenceEnd}, - } - indicesToRemove.Insert(segmentRec.Key, segmentRec.Val) + sequencesToRemove.Insert(fmt.Sprint(segment), [2]int{sequenceStart, sequenceEnd}) // skip past matched section and reset visitedPoints i = sequenceEnd visitedPoints = [][2]float64{} @@ -557,67 +683,28 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { sequenceEnd = start + 2*(len(segment)-1)*len(reverseMatches) endPointIdx = start + matches[len(matches)-1] + len(segment) } - segmentRec := sortedmap.Record{ - Key: fmt.Sprintf("%v", segment), - Val: [2]int{sequenceStart, sequenceEnd}, - } - indicesToRemove.Insert(segmentRec.Key, segmentRec.Val) - // check if remaining points are on a straight line, retain only start and end points if so - startPointX := newRing[sequenceEnd][0] - startPointY := newRing[sequenceEnd][1] - onLine := true - furthestDistance := 0.0 - furthestPointIdx := sequenceEnd - for n := sequenceEnd + 1; n < endPointIdx; n++ { - if newRing[n][0] != startPointX && newRing[n][1] != startPointY { - onLine = false - break - } - pointDistance := math.Sqrt(math.Abs(newRing[n][0]-startPointX) + math.Abs(newRing[n][1]-startPointY)) - if pointDistance > furthestDistance { - furthestDistance = pointDistance - furthestPointIdx = n - } - } - if onLine { - // remove any intermediate points on the return from the end point to the start point - sequenceStart := furthestPointIdx + 1 - sequenceEnd := endPointIdx - 1 - segmentRec := sortedmap.Record{ - Key: fmt.Sprintf("%v", newRing[furthestPointIdx:furthestPointIdx+1]), - Val: [2]int{sequenceStart, sequenceEnd}, - } - indicesToRemove.Insert(segmentRec.Key, segmentRec.Val) - } + sequencesToRemove.Insert(fmt.Sprint(segment), [2]int{sequenceStart, sequenceEnd}) + // (checking if remaining points are on a straight line could be done here + // but is not necessary because snap always inserts it) // skip past matched section and reset visitedPoints i = endPointIdx - 1 visitedPoints = [][2]float64{} } } - offset := 0 - for _, key := range indicesToRemove.Keys() { - sequence := indicesToRemove.Map()[key].([2]int) - deduplicatedRing = append(deduplicatedRing[:sequence[0]-offset], deduplicatedRing[sequence[1]-offset:]...) - offset += sequence[1] - sequence[0] - } - return deduplicatedRing + return removeSequences(ring, sequencesToRemove) } -func findRightmostLowestPoint(ring [][2]float64) (int, [2]float64) { - rightmostLowestPoint := [2]float64{math.MinInt, math.MaxInt} - j := 0 - for i, vertex := range ring { - // check if vertex is rightmost lowest point (either y is lower, or y is equal and x is higher) - if vertex[1] < rightmostLowestPoint[1] { - rightmostLowestPoint[0] = vertex[0] - rightmostLowestPoint[1] = vertex[1] - j = i - } else if vertex[1] == rightmostLowestPoint[1] && vertex[0] > rightmostLowestPoint[0] { - rightmostLowestPoint[0] = vertex[0] - j = i - } +func removeSequences(ring [][2]float64, sequencesToRemove *sortedmap.SortedMap[string, [2]int]) (newRing [][2]float64) { + mmap := sequencesToRemove.Map() + keepFrom := 0 + for _, key := range sequencesToRemove.Keys() { + sequenceToRemove := mmap[key] + keepTo := sequenceToRemove[0] + newRing = append(newRing, ring[keepFrom:keepTo]...) + keepFrom = sequenceToRemove[1] } - return j, rightmostLowestPoint + newRing = append(newRing, ring[keepFrom:]...) + return newRing } // repeatedly calls kmpSearch, returning all starting indexes of 'find' in 'corpus' @@ -700,23 +787,64 @@ func asKeys[T constraints.Ordered](elements []T) map[T]any { return mapped } -func findKeyWithMaxValue[K comparable, V constraints.Ordered](m map[K]V) (maxK K, maxV V, winners uint) { +func findLastKeyWithMaxValue[K comparable, V constraints.Ordered](m *orderedmap.OrderedMap[K, V]) (maxK K, maxV V, numWinners uint) { first := true - for k, v := range m { - if first || v > maxV { - maxK = k - maxV = v - winners = 1 + for p := m.Newest(); p != nil; p = p.Prev() { + if first || p.Value > maxV { + maxK = p.Key + maxV = p.Value + numWinners = 1 first = false continue } - if v == maxV { - winners++ + if p.Value == maxV { + numWinners++ } } return } +func orderedMapKeys[K comparable, V any](m *orderedmap.OrderedMap[K, V]) []K { + l := make([]K, m.Len()) + i := 0 + for p := m.Oldest(); p != nil; p = p.Next() { + l[i] = p.Key + i++ + } + return l +} + +func lastMatch[T comparable](haystack, needle []T) T { + for i := len(haystack) - 1; i >= 0; i-- { + if slices.Contains(needle, haystack[i]) { + return haystack[i] + } + } + var empty T + return empty +} + +func deleteFromSliceByIndex[V any, X any](s []V, indexesToDelete map[int]X, indexOffset int) []V { + r := make([]V, 0, len(s)) + for i := range s { + if _, skip := indexesToDelete[i+indexOffset]; skip { + continue + } + r = append(r, s[i]) + } + return r +} + +func countVals[K, V comparable](m *orderedmap.OrderedMap[K, V], v V) int { + n := 0 + for p := m.Oldest(); p != nil; p = p.Next() { + if p.Value == v { + n++ + } + } + return n +} + func floatPolygonsToGeomPolygonsForAllLevels(floatersPerLevel map[Level][][][][2]float64) map[Level][]geom.Polygon { geomsPerLevel := make(map[Level][]geom.Polygon, len(floatersPerLevel)) for l := range floatersPerLevel { @@ -778,6 +906,29 @@ func panicMoreThanOneMatchingOuterRing(polygons [][][][2]float64, innerRing [][2 panic(panicMsg) } +func panicDeletedPolygonHasInnerRing(polygon [][][2]float64) { + panicMsg := fmt.Sprintf("a deleted dupe polygon had more than one ring, %v", + truncatedWkt(floatPolygonToGeomPolygon(polygon), 100)) + panic(panicMsg) +} + func truncatedWkt(geom geom.Geometry, width uint) string { return truncate.StringWithTail(wkt.MustEncode(geom), width, "...") } + +func reverseClone[S ~[]E, E any](s S) S { + if s == nil { + return nil + } + l := len(s) + c := make(S, l) + for i := 0; i < l; i++ { + c[l-1-i] = s[i] + } + return c +} + +func roundFloat(f float64, p uint) float64 { + r := math.Pow(10, float64(p)) + return math.Round(f*r) / r +} diff --git a/snap/snap_test.go b/snap/snap_test.go index 5781461..faa0994 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -3,7 +3,7 @@ package snap import ( "testing" - "log" + "github.com/go-spatial/geom/encoding/wkt" "github.com/pdok/texel/tms20" @@ -549,7 +549,7 @@ func TestSnap_snapPolygon(t *testing.T) { }}, }, { - name: "outer ring with horizontal rightmostlowest", + name: "inner but no outer error because of not reversing because of horizontal rightmostlowest", tms: loadEmbeddedTileMatrixSet(t, "NetherlandsRDNewQuad"), tmIDs: []tms20.TMID{5}, polygon: geom.Polygon{{ @@ -564,16 +564,96 @@ func TestSnap_snapPolygon(t *testing.T) { geom.Polygon{{{139520.48, 527777.44}, {139513.76, 527777.44}}}, }}, // want no panicInnerRingsButNoOuterRings }, + { + name: "inner but no outer error because of not reversing because of a very sharp leg/extension", + tms: loadEmbeddedTileMatrixSet(t, "NetherlandsRDNewQuad"), + tmIDs: []tms20.TMID{0}, + polygon: geom.Polygon{{ + {48158.204, 392310.062}, + {47753.125, 391885.44}, {48565.4, 391515.876}, {47751.195, 391884.821}, // a very sharp leg/extension + {47677.592, 392079.403}, + }}, + want: map[tms20.TMID][]geom.Polygon{0: { + {{{47587.52, 392144.32}, {47802.56, 391929.28}, {48232.64, 392359.36}}}, // turned counterclockwise + {{{47802.56, 391929.28}, {48662.72, 391499.2}}}, + }}, // want no panicInnerRingsButNoOuterRings + }, + { + name: "split ring from outer is cw, should be ccw", + tms: loadEmbeddedTileMatrixSet(t, "NetherlandsRDNewQuad"), + tmIDs: []tms20.TMID{0}, + polygon: geom.Polygon{{{179334.089, 408229.072}, {179121.631, 408528.181}, {179328.228, 408231.924}, {178889.903, 408431.167}, {178531.386, 408106.618}, {178497.492, 407886.329}, {178535.353, 408103.574}, {178862.244, 408226.852}, {178891.816, 408426.547}, {179173.349, 408187.199}, {178893.957, 408423.424}, {178864.491, 408223.293}, {178537.744, 408101.003}, {178504.209, 407887.598}, {178510.008, 407890.491}, {178542.44, 408098.473}, {178867.788, 408219.534}, {178897.835, 408417.763}, {179170.131, 408181.285}}}, // ccw (with a sharp leg/extension also) + want: map[tms20.TMID][]geom.Polygon{0: { + {{{178976.96, 408487.36}, {178761.92, 408272.32}, {178546.88, 408057.28}, {179192, 408272.32}}}, // ccw + // bunch of points and lines: + {{{179407.04, 408272.32}, {179192, 408272.32}}}, {{{179192, 408272.32}, {179192, 408487.36}}}, {{{179407.04, 408272.32}, {179192, 408272.32}}}, {{{179192, 408272.32}, {178976.96, 408487.36}}}, {{{178976.96, 408487.36}, {178761.92, 408272.32}}}, {{{178761.92, 408272.32}, {178546.88, 408057.28}}}, {{{178546.88, 408057.28}, {178546.88, 407842.24}}}, + }}, + }, + { + name: "one of three split outer rings is cw and turned outer after no matching outer", + tms: loadEmbeddedTileMatrixSet(t, "NetherlandsRDNewQuad"), + tmIDs: []tms20.TMID{0}, + polygon: geom.Polygon{{{88580.011, 439678.996}, {88337.73, 439237.216}, {89273.964, 438026.4}, {89386.079, 438023.335}, {90251.524, 438784.15}, {89852.567, 439284.421}, {89425.263, 439355.284}, {89247.228, 439563.507}, {89089.95, 439692.364}, {88959.832, 439729.531}, {89055.886, 439819.684}, {89466.904, 439382.346}, {89899.488, 439311.969}, {90170.183, 438911.775}, {90329.354, 438821.391}, {90651.094, 438796.963}, {91473.854, 439243.296}, {90632.307, 438747.518}, {90270.708, 438757.632}, {89555.357, 437677.283}, {90499.163, 436096.427}, {91435.651, 435963.019}, {91404.334, 436039.088}, {91254.337, 436091.084}, {90500.745, 436098.362}, {90076.214, 437042.706}, {89870.055, 437307.816}, {89768.94, 437363.42}, {89650.683, 437521.434}, {89640.994, 437568.838}, {89558.222, 437677.647}, {90269.467, 438753.387}, {90632.85, 438744.94}, {91313.174, 439143.369}, {91477.748, 439241.657}, {91475.353, 439245.66}, {91457.592, 439266.852}, {91243.008, 439179.921}, {90710.843, 438897.924}, {90650.175, 438799.288}, {90440.729, 438846.985}, {90395.019, 438846.967}, {90329.938, 438823.822}, {90287.474, 438885.328}, {90172.086, 438913.396}, {90044.257, 439125.421}, {89901.052, 439313.924}, {89885.113, 439321.991}, {89835.824, 439335.083}, {89468.228, 439384.467}, {89173.832, 439758.873}, {89061.413, 439821.909}, {89054.68, 439821.883}, {89023.222, 439801.24}, {88989.659, 439763.597}, {88949.781, 439739.428}, {88958.959, 439726.203}, {89088.39, 439690.41}, {89245.45, 439561.75}, {89388.248, 439376.01}, {89424.081, 439353.075}, {89566.906, 439317.631}, {89851.03, 439282.45}, {90111.766, 438914.525}, {90249.027, 438784.029}, {90211.25, 438760.51}, {90183.492, 438736.293}, {89584.683, 438207.656}, {89384.579, 438025.335}, {89274.819, 438028.749}, {88339.974, 439238.317}, {88419.861, 439377.057}, {88447.454, 439387.602}, {88485.231, 439376.209}, {88505.9, 439379.802}, {88564.366, 439441.722}, {88589.428, 439478.721}, {88598.844, 439504.106}, {88608.517, 439561.563}, {88582.418, 439679.669}, {88565.692, 439724.97}, {88480.367, 439857.335}, {88409.981, 439938.527}, {88412.431, 439940.265}, {88366.171, 440033.682}, {88353.723, 440046.457}, {88356.08, 440054.25}, {88342.856, 440086.861}, {88266.552, 440224.799}, {88252.681, 440243.646}, {88196.44, 440306.135}, {87992.789, 440467.453}, {88250.595, 440274.14}, {88508.083, 439845.775}, {88270.249, 440256.888}, {88194.893, 440335.659}, {88010.485, 440474.349}, {87996.213, 440475.679}, {87990.894, 440469.07}, {88580.011, 439678.996}}}, + want: map[tms20.TMID][]geom.Polygon{}, // want no panicNoMatchingOuterForInnerRing + }, + { + name: "sneaky nested pseudo ring creates more than 1 matching outer ring", + tms: loadEmbeddedTileMatrixSet(t, "NetherlandsRDNewQuad"), + tmIDs: []tms20.TMID{0}, + polygon: geom.Polygon{ + {{198877.1, 506188.635}, {198805.608, 506361.231}, {198633.011, 506432.722}, {198460.415, 506361.23}, {198388.924, 506188.633}, {198460.416, 506016.037}, {198633.013, 505944.546}, {198805.609, 506016.038}}, + {{198429.407, 506188.635}, {198489.229, 506332.782}, {198633.531, 506392.228}, {198777.528, 506332.045}, {198836.612, 506187.594}, {198776.434, 506044.111}, {198632.5, 505985.022}, {198488.864, 506044.832}, {198429.407, 506188.615}, {198551.204, 506045.823}, {198690.244, 506034.324}, {198792.36, 506147.487}, {198748.509, 506305.863}, {198576.128, 506343.056}}, + {{198633.012, 506279.536}, {198766.685, 506188.158}, {198632.396, 506055.195}, {198499.739, 506188.974}}}, + want: map[tms20.TMID][]geom.Polygon{}, // want no panicMoreThanOneMatchingOuterRing + }, + { + name: "nested rings", + tms: newSimpleTileMatrixSet(2, 64), + tmIDs: []tms20.TMID{ + 1, // 32 * 8.0 + }, + polygon: geom.Polygon{ + {{4.0, 124.0}, {4.0, 4.0}, {60.0, 4.0}, {60.0, 124.0}}, // big outer + {{12.0, 52.0}, {12.0, 12.0}, {52.0, 12.0}, {52.0, 52.0}, {30.0, 52.0}, {30.0, 44.0}, {44.0, 44.0}, {44.0, 20.0}, {20.0, 20.0}, {20.0, 44.0}, {27.0, 44.0}, {27.0, 52.0}}, // big letter C that turns into nested rings when snapped + {{12.0, 116.0}, {52.0, 116.0}, {52.0, 76.0}, {30.0, 76.0}, {30.0, 84.0}, {44.0, 84.0}, {44.0, 108.0}, {20.0, 108.0}, {20.0, 84.0}, {27.0, 84.0}, {27.0, 76.0}, {12.0, 76.0}}, // above mirrored vertically + {{30.0, 53.0}, {30.0, 54.0}, {54.0, 54.0}, {54.0, 10.0}, {10.0, 10.0}, {10.0, 54.0}, {27.0, 54.0}, {27.0, 53.0}, {11.0, 53.0}, {11.0, 11.0}, {53.0, 11.0}, {53.0, 53.0}}, // another C around the original C which snaps to a duplicate outer + {{28.0, 28.0}, {36.0, 28.0}, {36.0, 36.0}, {29.0, 36.0}, {29.0, 92.0}, {36.0, 92.0}, {36.0, 100.0}, {28.0, 100.0}}, // dumbbell inside the two C's that also turns into a nested ring when snapped (and some lines) + }, + want: map[tms20.TMID][]geom.Polygon{ + // 1 big outer with 2 holes. and 2 new outers/polygons (inside those holes from the big outer) each with their own hole. + // a duplicate inner/outer pair from the extra C is removed + 1: { + { + {{4.0, 124.0}, {4.0, 4.0}, {60.0, 4.0}, {60.0, 124.0}}, // ccw + {{12.0, 116.0}, {52.0, 116.0}, {52.0, 76.0}, {28.0, 76.0}, {12.0, 76.0}}, // cw + {{28.0, 52.0}, {52.0, 52.0}, {52.0, 12.0}, {12.0, 12.0}, {12.0, 52.0}}, // cw + }, { + {{28.0, 44.0}, {20.0, 44.0}, {20.0, 20.0}, {44.0, 20.0}, {44.0, 44.0}}, // ccw + {{28.0, 36.0}, {36.0, 36.0}, {36.0, 28.0}, {28.0, 28.0}}, // cw + }, { + {{28.0, 84.0}, {44.0, 84.0}, {44.0, 108.0}, {20.0, 108.0}, {20.0, 84.0}}, // ccw + {{28.0, 100.0}, {36.0, 100.0}, {36.0, 92.0}, {28.0, 92.0}}, // cw + }, + // and some lines + {{{28.0, 52.0}, {28.0, 44.0}}}, + {{{28.0, 76.0}, {28.0, 84.0}}}, + {{{28.0, 92.0}, {28.0, 84.0}}}, + {{{28.0, 84.0}, {28.0, 76.0}}}, + {{{28.0, 76.0}, {28.0, 52.0}}}, + {{{28.0, 52.0}, {28.0, 44.0}}}, + {{{28.0, 44.0}, {28.0, 36.0}}}, + }, + }, + }, } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { got := SnapPolygon(tt.polygon, tt.tms, tt.tmIDs) for tmID, wantPoly := range tt.want { if !assert.EqualValues(t, wantPoly, got[tmID]) { - // t.Errorf("snapPolygon(...) = %v, want %v", wkt.MustEncode(got[tmID]), wkt.MustEncode(wantPoly)) - t.Errorf("snapPolygon(...) = %v, want %v", got[tmID], wantPoly) + t.Errorf("snapPolygon(%v, _, %v)\n= %v\nwant: %v", + wkt.MustEncode(tt.polygon), tmID, wktMustEncodePolygonSlice(got[tmID]), wktMustEncodePolygonSlice(wantPoly)) } - log.Printf("snapPolygon(...) = %v", got[tmID]) } }) } @@ -608,3 +688,264 @@ func TestSnap_ringContains(t *testing.T) { }) } } + +func Test_kmpDeduplicate(t *testing.T) { + tests := []struct { + name string + ring [][2]float64 + want [][2]float64 + }{ + { + name: "triangle should stay", + ring: [][2]float64{ + {2, 1}, // A + {1, 1}, // B + {1, 0}, // C + {1, 1}, // B + {0, 1}, // D + {1, 0}, // C + {1, 1}, // B + }, + want: [][2]float64{ + {2, 1}, // A + {1, 1}, // B + {0, 1}, // D + {1, 0}, // C + {1, 1}, // B + }, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + assert.Equalf(t, tt.want, kmpDeduplicate(tt.ring), "kmpDeduplicate(%v)", tt.ring) + }) + } +} + +func Test_dedupeInnersOuters(t *testing.T) { + type args struct { + outers [][][2]float64 + inners [][][2]float64 + } + tests := []struct { + name string + args args + wantOuters [][][2]float64 + wantInners [][][2]float64 + }{ + { + name: "#outer, #inner = 0, 0", + args: args{ + outers: squareRingArray(0, true), + inners: squareRingArray(0, false), + }, + wantOuters: squareRingArray(0, true), + wantInners: squareRingArray(0, false), + }, + { + name: "#outer, #inner = 1, 0", + args: args{ + outers: squareRingArray(1, true), + inners: squareRingArray(0, false), + }, + wantOuters: squareRingArray(1, true), + wantInners: squareRingArray(0, false), + }, + { + name: "#outer, #inner = 1, 1", + args: args{ + outers: squareRingArray(1, true), + inners: squareRingArray(1, false), + }, + wantOuters: squareRingArray(1, true), + wantInners: squareRingArray(1, false), + }, + { + name: "#outer, #inner = 2, 1", + args: args{ + outers: squareRingArray(2, true), + inners: squareRingArray(1, false), + }, + wantOuters: squareRingArray(1, true), + wantInners: squareRingArray(0, false), + }, + { + name: "#outer, #inner = 2, 2", + args: args{ + outers: squareRingArray(2, true), + inners: squareRingArray(2, false), + }, + wantOuters: squareRingArray(1, true), + wantInners: squareRingArray(1, false), + }, + { + name: "#outer, #inner = 0, 1", + args: args{ + outers: squareRingArray(0, true), + inners: squareRingArray(1, false), + }, + wantOuters: squareRingArray(0, true), + wantInners: squareRingArray(1, false), + }, + { + name: "#outer, #inner = 1, 2", + args: args{ + outers: squareRingArray(1, true), + inners: squareRingArray(2, false), + }, + wantOuters: squareRingArray(0, true), + wantInners: squareRingArray(1, false), + }, + { + name: "#outer, #inner = 2, 0", + args: args{ + outers: squareRingArray(2, true), + inners: squareRingArray(0, false), + }, + wantOuters: squareRingArray(2, true), + wantInners: squareRingArray(0, false), + }, + { + name: "#outer, #inner = 0, 2", + args: args{ + outers: squareRingArray(0, true), + inners: squareRingArray(2, false), + }, + wantOuters: squareRingArray(0, true), + wantInners: squareRingArray(2, false), + }, + { + name: "#outer, #inner = 3, 1", + args: args{ + outers: squareRingArray(3, true), + inners: squareRingArray(1, false), + }, + wantOuters: squareRingArray(2, true), + wantInners: squareRingArray(0, false), + }, + { + name: "#outer, #inner = 1, 3", + args: args{ + outers: squareRingArray(1, true), + inners: squareRingArray(3, false), + }, + wantOuters: squareRingArray(0, true), + wantInners: squareRingArray(2, false), + }, + // add different ring, which should be left alone + { + name: "#outer, #inner = 1, 1 + dummy ring", + args: args{ + outers: squareRingArray(1, true), + inners: append(squareRingArray(1, false), [][2]float64{{0, 0}, {1, 0}, {2, 1}}), + }, + wantOuters: squareRingArray(1, true), + wantInners: append(squareRingArray(1, false), [][2]float64{{0, 0}, {1, 0}, {2, 1}}), + }, + { + name: "#outer, #inner = 3, 1 + dummy ring", + args: args{ + outers: squareRingArray(3, true), + inners: append(squareRingArray(1, false), [][2]float64{{0, 0}, {1, 0}, {2, 1}}), + }, + wantOuters: squareRingArray(2, true), + wantInners: append(squareRingArray(0, false), [][2]float64{{0, 0}, {1, 0}, {2, 1}}), + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + got, got1 := dedupeInnersOuters(tt.args.outers, tt.args.inners) + assert.Equalf(t, tt.wantOuters, got, "dedupeInnersOuters(%v, %v)", tt.args.outers, tt.args.inners) + assert.Equalf(t, tt.wantInners, got1, "dedupeInnersOuters(%v, %v)", tt.args.outers, tt.args.inners) + }) + } +} + +// newSimpleTileMatrixSet creates a tms for snap testing purposes +// the effective quadrant amount (one axis) on the deepest level will be 2^maxDepth * 16 (vt internal pixel res) +// the effective quadrant size (one axis) on the deepest level will be cellSize / 16 (vt internal pixel res) +func newSimpleTileMatrixSet(maxDepth Level, cellSize float64) tms20.TileMatrixSet { + zeroZero := tms20.TwoDPoint([2]float64{0.0, 0.0}) + tms := tms20.TileMatrixSet{ + CRS: fakeCRS{}, + OrderedAxes: []string{"X", "Y"}, + TileMatrices: make(map[tms20.TMID]tms20.TileMatrix, maxDepth+1), + } + for tmID := 0; tmID <= int(maxDepth); tmID++ { + // (only values from the root tm are used, for the rest it is assumed to follow quad matrix rules) + tmCellSize := cellSize * float64(pow2(maxDepth-uint(tmID))) + tms.TileMatrices[tmID] = tms20.TileMatrix{ + ID: "0", + ScaleDenominator: tmCellSize / tms20.StandardizedRenderingPixelSize, + CellSize: tmCellSize, + CornerOfOrigin: tms20.BottomLeft, + PointOfOrigin: &zeroZero, + TileWidth: 1, + TileHeight: 1, + MatrixWidth: 1, + MatrixHeight: 1, + } + } + return tms +} + +func wktMustEncodePolygonSlice(geoms []geom.Polygon) string { + s := "" + for i := range geoms { + s += wktMustEncode(geoms[i]) + "\n" + } + return s +} + +func wktMustEncode(g geom.Geometry) (s string) { + p, isPoly := g.(geom.Polygon) + if !isPoly { + return wkt.MustEncode(g) + } + + var lines []geom.LineString + var points []geom.Point + pp := make(geom.Polygon, len(p)) + copy(pp, p) + for r := 0; r < len(pp); r++ { + switch len(pp[r]) { + default: + continue + case 1: + points = append(points, pp[r][0]) + case 2: + lines = append(lines, pp[r]) + } + pp = append(pp[:r], pp[r+1:]...) + r-- + } + + if len(pp) > 0 { + s = wkt.MustEncode(pp) + } + for i := range lines { + s += wkt.MustEncode(lines[i]) + } + for i := range points { + s += wkt.MustEncode(points[i]) + } + return s +} + +func squareRingArray(number int, isOuter bool) [][][2]float64 { + outerSquare := [][2]float64{{0, 0}, {1, 0}, {1, 1}, {0, 1}} // square, counter clockwise + innerSquare := [][2]float64{{0, 0}, {0, 1}, {1, 1}, {1, 0}} // square, clockwise + var squares = [][][2]float64{} + var square [][2]float64 + // outer or inner + if isOuter { + square = outerSquare + } else { + square = innerSquare + } + // add squares + for i := 0; i < number; i++ { + squares = append(squares, square) + } + return squares +}