From 1d74236551830bcd4b0c3143083023bf5b568402 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Mon, 29 Jan 2024 16:03:00 +0100 Subject: [PATCH 01/20] round floats in isClockWise more precise --- snap/snap.go | 13 ++++++++++--- snap/snap_test.go | 16 +++++++++++++++- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 4b39b0e..fdee94a 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -19,8 +19,9 @@ import ( ) const ( - keepPointsAndLines = true // TODO do something with polys that collapsed into points and lines - internalPixelResolution = 16 + keepPointsAndLines = true // TODO do something with polys that collapsed into points and lines + internalPixelResolution = 16 + roundFloatAgainstFPErrorPrecision = 5 ) // SnapPolygon snaps polygons' points to a tile's internal pixel grid @@ -443,7 +444,8 @@ func isClockwise(points [3][2]float64) bool { if vector1.angle() == 0.0 { return relativeAngle > 180 } - return math.Round(mod((vector2.angle()-relativeAngle), 360)) != math.Round(vector1.angle()) + return roundFloat(mod(vector2.angle()-relativeAngle, 360), roundFloatAgainstFPErrorPrecision) != + roundFloat(vector1.angle(), roundFloatAgainstFPErrorPrecision) } // deduplication using an implementation of the Knuth-Morris-Pratt algorithm @@ -781,3 +783,8 @@ func panicMoreThanOneMatchingOuterRing(polygons [][][][2]float64, innerRing [][2 func truncatedWkt(geom geom.Geometry, width uint) string { return truncate.StringWithTail(wkt.MustEncode(geom), width, "...") } + +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..182cb9f 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -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,6 +564,20 @@ 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 + }, } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { From 4e7e6645ec2b46d2ce33a3239f36c2dd6d656d06 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Tue, 30 Jan 2024 16:56:26 +0100 Subject: [PATCH 02/20] swap outer and inner rings after splitting if incorrect winding order and use geom.winding Co-authored-by: Michiel Korpel --- snap/snap.go | 64 +++++++++++++++-------------------------------- snap/snap_test.go | 11 ++++++++ 2 files changed, 31 insertions(+), 44 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index fdee94a..41c2a2d 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -6,6 +6,8 @@ import ( "slices" "sort" + "github.com/go-spatial/geom/winding" + "github.com/go-spatial/geom/encoding/wkt" "github.com/muesli/reflow/truncate" @@ -19,9 +21,8 @@ import ( ) const ( - keepPointsAndLines = true // TODO do something with polys that collapsed into points and lines - internalPixelResolution = 16 - roundFloatAgainstFPErrorPrecision = 5 + keepPointsAndLines = true // TODO do something with polys that collapsed into points and lines + internalPixelResolution = 16 ) // SnapPolygon snaps polygons' points to a tile's internal pixel grid @@ -327,14 +328,8 @@ func ensureCorrectWindingOrder(ring [][2]float64, shouldBeClockwise bool) { // 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 @@ -429,23 +424,21 @@ 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) + 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) + innerRings = append(innerRings, outerRing) + } + outerRings = make([][][2]float64, 0) } - return roundFloat(mod(vector2.angle()-relativeAngle, 360), roundFloatAgainstFPErrorPrecision) != - roundFloat(vector1.angle(), roundFloatAgainstFPErrorPrecision) + return outerRings, innerRings, pointsAndLines } // deduplication using an implementation of the Knuth-Morris-Pratt algorithm @@ -605,23 +598,6 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { return deduplicatedRing } -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 - } - } - return j, rightmostLowestPoint -} - // repeatedly calls kmpSearch, returning all starting indexes of 'find' in 'corpus' func kmpSearchAll(corpus, find [][2]float64) []int { matches := []int{} diff --git a/snap/snap_test.go b/snap/snap_test.go index 182cb9f..884d2bf 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -578,6 +578,17 @@ func TestSnap_snapPolygon(t *testing.T) { {{{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}}}, + }}, + }, } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { From 6436299067ef72e36c75e4bc7550284c1bee0177 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Tue, 30 Jan 2024 18:05:41 +0100 Subject: [PATCH 03/20] no need to check for points on line in dedupe algo and a little readability improvement PDOK-16135 Co-authored-by: Michiel Korpel --- go.mod | 4 +- go.sum | 4 +- snap/snap.go | 149 ++++++++++++++++++---------------------------- snap/snap_test.go | 33 ++++++++++ 4 files changed, 96 insertions(+), 94 deletions(-) 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/snap.go b/snap/snap.go index 41c2a2d..73af429 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -14,7 +14,7 @@ import ( "github.com/go-spatial/geom" "github.com/pdok/texel/intgeom" "github.com/pdok/texel/tms20" - "github.com/umpc/go-sortedmap" + "github.com/tobshub/go-sortedmap" orderedmap "github.com/wk8/go-ordered-map/v2" "golang.org/x/exp/constraints" "golang.org/x/exp/maps" @@ -219,58 +219,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 @@ -336,10 +336,10 @@ func windingOrderIsCorrect(ring [][2]float64, shouldBeClockwise bool) bool { 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 @@ -444,17 +444,16 @@ func splitRing(ring [][2]float64, isOuter bool, hitMultiple map[intgeom.Point][] // 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) @@ -465,7 +464,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 @@ -481,7 +480,7 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { 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 @@ -491,8 +490,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 { @@ -500,7 +499,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 @@ -512,11 +511,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{} @@ -525,11 +520,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{} @@ -552,50 +543,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 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] + } + newRing = append(newRing, ring[keepFrom:]...) + return newRing } // repeatedly calls kmpSearch, returning all starting indexes of 'find' in 'corpus' diff --git a/snap/snap_test.go b/snap/snap_test.go index 884d2bf..c19f0ba 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -633,3 +633,36 @@ 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) + }) + } +} From 8c8878edfa81a09edf53e406bca71b4298a95e90 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Mon, 5 Feb 2024 08:57:40 +0100 Subject: [PATCH 04/20] re-reverse inner-turned-outer if no matching outer is found PDOK-16135 --- snap/snap.go | 84 ++++++++++++++++++++++++++++++++++++++++------- snap/snap_test.go | 7 ++++ 2 files changed, 79 insertions(+), 12 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 73af429..bbb2514 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -2,6 +2,7 @@ package snap import ( "fmt" + "log" "math" "slices" "sort" @@ -131,11 +132,9 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ for _, outerRing := range outerRings { // (there are only outer rings if isOuter) newPolygons[level] = append(newPolygons[level], [][][2]float64{outerRing}) } - if len(newPolygons[level]) == 0 && len(innerRings) > 0 { // should never happen - panicInnerRingsButNoOuterRings(level, polygon, innerRings) - continue + if len(innerRings) > 0 { + newPolygons[level] = matchInnersToPolygon(newPolygons[level], innerRings, len(polygon) > 1) } - newPolygons[level] = matchInnersToPolygon(newPolygons[level], innerRings) if keepPointsAndLines { newPointsAndLines[level] = append(newPointsAndLines[level], pointsAndLines...) } @@ -150,14 +149,14 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ return floatPolygonsToGeomPolygonsForAllLevels(newPolygons) } -func matchInnersToPolygon(polygons [][][][2]float64, innerRings [][][2]float64) [][][][2]float64 { +func matchInnersToPolygon(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 } + numDeduped, polygons := dedupePolygonsByOuters(polygons) + + var innersTurnedOuters [][][2]float64 matchInners: for _, innerRing := range innerRings { containsPerPolyI := make(map[int]uint, lenPolygons) @@ -179,17 +178,60 @@ matchInners: continue matchInners } } - // no (single) matching outer ring was found (should never happen) if len(containsPerPolyI) == 0 { - panicNoMatchingOuterForInnerRing(polygons, innerRing) + // 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) + 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) + // multiple matching outer rings were found, possibly because there are duplicates + panicMoreThanOneMatchingOuterRing(polygons, innerRing, numDeduped) continue } + for i := range innersTurnedOuters { + polygons = append(polygons, [][][2]float64{innersTurnedOuters[i]}) + } return polygons } +// helper for matchInnersToPolygon to delete duplicate polygons +// comparing them only by their outers and asserting that a deleted polygon didn't have inner rings appended yet +// yes it's implemented as ~O(n^2), +// but it's expected that the (outer) rings are usually different even from the first point, +// making it still more efficient than using a hashmap of the entire rings +func dedupePolygonsByOuters(polygons [][][][2]float64) (int, [][][][2]float64) { + numPolygons := len(polygons) + numDeleted := 0 + for i := 0; i < numPolygons; i++ { + ring := polygons[i][0] // only check the outer + compareToOther: + for j := i + 1; j < numPolygons; j++ { + ringLen := len(ring) + other := polygons[j][0] + otherLen := len(other) + if ringLen != otherLen { + continue + } + for k := 0; k < min(ringLen, otherLen); k++ { + if ring[k] != other[k] { + continue compareToOther + } + } + // delete + if len(polygons[j]) > 1 { + panicDeletedPolygonHasInnerRing(polygons[j]) + } + polygons = append(polygons[:j], polygons[j+1:]...) + j-- + numPolygons-- + numDeleted++ + } + } + return numDeleted, polygons +} + // 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. @@ -725,10 +767,28 @@ 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 c19f0ba..85a0ff3 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -589,6 +589,13 @@ func TestSnap_snapPolygon(t *testing.T) { {{{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 + }, } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { From 56ff19333c3acf8b364024f59e966bedd5a7958e Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Mon, 5 Feb 2024 08:58:50 +0100 Subject: [PATCH 05/20] test (dupe) inner rings, one turns outer, multi-match PDOK-16135 --- snap/snap_test.go | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/snap/snap_test.go b/snap/snap_test.go index 85a0ff3..7080b01 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -596,6 +596,16 @@ func TestSnap_snapPolygon(t *testing.T) { 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}, {198877.1, 506188.635}}, + {{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}, {198429.407, 506188.635}}, + {{198633.012, 506279.536}, {198766.685, 506188.158}, {198632.396, 506055.195}, {198499.739, 506188.974}, {198633.012, 506279.536}}}, + want: map[tms20.TMID][]geom.Polygon{}, // want no panicMoreThanOneMatchingOuterRing + }, } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { From a2b8386cb14406688e44c3bc73da77b56ef5b138 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Sun, 11 Feb 2024 17:29:20 +0100 Subject: [PATCH 06/20] shoelace Co-authored-by: WouterVisscher --- snap/geom.go | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 snap/geom.go 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) +} From 440cd141878971dc91c138e77a0fc47d6b7fb42b Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Mon, 12 Feb 2024 11:08:37 +0100 Subject: [PATCH 07/20] match inner rings to smallest outer rings also dont reverse rings in place PDOK-16135 --- snap/pointindex_test.go | 25 -------- snap/snap.go | 92 +++++++++++++++++++++--------- snap/snap_test.go | 122 +++++++++++++++++++++++++++++++++++++--- 3 files changed, 179 insertions(+), 60 deletions(-) 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 bbb2514..04e85e9 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -102,7 +102,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 { @@ -154,41 +154,42 @@ func matchInnersToPolygon(polygons [][][][2]float64, innerRings [][][2]float64, if len(innerRings) == 0 { return polygons } - numDeduped, polygons := dedupePolygonsByOuters(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 } } - if len(containsPerPolyI) == 0 { + 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 } - // multiple matching outer rings were found, possibly because there are duplicates - panicMoreThanOneMatchingOuterRing(polygons, innerRing, numDeduped) - 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]}) @@ -196,6 +197,20 @@ matchInners: 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() +} + // helper for matchInnersToPolygon to delete duplicate polygons // comparing them only by their outers and asserting that a deleted polygon didn't have inner rings appended yet // yes it's implemented as ~O(n^2), @@ -362,10 +377,11 @@ 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) @@ -469,13 +485,13 @@ func splitRing(ring [][2]float64, isOuter bool, hitMultiple map[intgeom.Point][] // 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) + 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) + slices.Reverse(outerRing) // in place, not used elsewhere innerRings = append(innerRings, outerRing) } outerRings = make([][][2]float64, 0) @@ -516,7 +532,7 @@ func kmpDeduplicate(ring [][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) @@ -689,23 +705,43 @@ 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 floatPolygonsToGeomPolygonsForAllLevels(floatersPerLevel map[Level][][][][2]float64) map[Level][]geom.Polygon { geomsPerLevel := make(map[Level][]geom.Polygon, len(floatersPerLevel)) for l := range floatersPerLevel { diff --git a/snap/snap_test.go b/snap/snap_test.go index 7080b01..3da9190 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -1,9 +1,10 @@ package snap import ( + "fmt" "testing" - "log" + "github.com/go-spatial/geom/encoding/wkt" "github.com/pdok/texel/tms20" @@ -601,21 +602,57 @@ func TestSnap_snapPolygon(t *testing.T) { 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}, {198877.1, 506188.635}}, - {{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}, {198429.407, 506188.635}}, - {{198633.012, 506279.536}, {198766.685, 506188.158}, {198632.396, 506055.195}, {198499.739, 506188.974}, {198633.012, 506279.536}}}, + {{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{ // TODO think of something with duplicate outer rings + {{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 + {{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. + 1: { + { + {{4.0, 124.0}, {4.0, 4.0}, {60.0, 4.0}, {60.0, 124.0}}, // ccw + {{28.0, 52.0}, {52.0, 52.0}, {52.0, 12.0}, {12.0, 12.0}, {12.0, 52.0}}, // cw + {{12.0, 116.0}, {52.0, 116.0}, {52.0, 76.0}, {28.0, 76.0}, {12.0, 76.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]) } }) } @@ -683,3 +720,74 @@ func Test_kmpDeduplicate(t *testing.T) { }) } } + +// 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 +} From 5c54f85b4fe84bb157edfa35412f875bc672ad2f Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Wed, 14 Feb 2024 12:12:41 +0100 Subject: [PATCH 08/20] WIP handle overlapping inner rings PDOK-16135 --- snap/snap.go | 21 +++++++++++++++------ snap/snap_test.go | 8 +++++++- 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 04e85e9..b9f1a9a 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -93,6 +93,7 @@ func tileMatrixIDsByLevels(tms tms20.TileMatrixSet, tmIDs []tms20.TMID) map[Leve func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[Level][]geom.Polygon { levelMap := asKeys(levels) newPolygons := 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. @@ -129,17 +130,23 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ delete(levelMap, level) // outer ring has become too small continue } - for _, outerRing := range outerRings { // (there are only outer rings if isOuter) + for _, outerRing := range outerRings { newPolygons[level] = append(newPolygons[level], [][][2]float64{outerRing}) } - if len(innerRings) > 0 { - newPolygons[level] = matchInnersToPolygon(newPolygons[level], innerRings, len(polygon) > 1) - } + newInners[level] = append(newInners[level], innerRings...) if keepPointsAndLines { newPointsAndLines[level] = append(newPointsAndLines[level], pointsAndLines...) } } } + + for l := range levelMap { + newPolygonsForLevel := matchInnersToPolygons(newPolygons[l], newInners[l], len(polygon) > 1) + if len(newPolygonsForLevel) > 1 { + newPolygons[l] = newPolygonsForLevel + } + } + // points and lines at the end, as outer rings for level, pointsAndLines := range newPointsAndLines { for _, pointOrLine := range pointsAndLines { @@ -149,11 +156,13 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ return floatPolygonsToGeomPolygonsForAllLevels(newPolygons) } -func matchInnersToPolygon(polygons [][][][2]float64, innerRings [][][2]float64, hasInners bool) [][][][2]float64 { +func matchInnersToPolygons(polygons [][][][2]float64, innerRings [][][2]float64, hasInners bool) [][][][2]float64 { lenPolygons := len(polygons) if len(innerRings) == 0 { return polygons } + _, polygons = dedupePolygonsByOuters(polygons) + var polyISortedByOuterAreaDesc []int var innersTurnedOuters [][][2]float64 matchInners: @@ -211,7 +220,7 @@ func sortPolyIdxsByOuterAreaDesc(polygons [][][][2]float64) []int { return areas.Keys() } -// helper for matchInnersToPolygon to delete duplicate polygons +// helper for matchInnersToPolygons to delete duplicate polygons // comparing them only by their outers and asserting that a deleted polygon didn't have inner rings appended yet // yes it's implemented as ~O(n^2), // but it's expected that the (outer) rings are usually different even from the first point, diff --git a/snap/snap_test.go b/snap/snap_test.go index 3da9190..441f01c 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -613,10 +613,11 @@ func TestSnap_snapPolygon(t *testing.T) { tmIDs: []tms20.TMID{ 1, // 32 * 8.0 }, - polygon: geom.Polygon{ // TODO think of something with duplicate outer rings + 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{ @@ -626,6 +627,11 @@ func TestSnap_snapPolygon(t *testing.T) { {{4.0, 124.0}, {4.0, 4.0}, {60.0, 4.0}, {60.0, 124.0}}, // ccw {{28.0, 52.0}, {52.0, 52.0}, {52.0, 12.0}, {12.0, 12.0}, {12.0, 52.0}}, // cw {{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}, {12.0, 52.0}, {12.0, 12.0}, {52.0, 12.0}, {52.0, 52.0}, {28.0, 52.0}}, // ccw + {{28.0, 52.0}, {52.0, 52.0}, {52.0, 12.0}, {12.0, 12.0}, {12.0, 52.0}, {28.0, 52.0}}, // cw + // TODO order of separate outer rings OK. does this overlap the deepest nested square? + // TODO deduplicate this, remove totally? }, { {{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 From 3114976000710e0aefa2c6e4c9106d80e65b0251 Mon Sep 17 00:00:00 2001 From: kad-dirc Date: Thu, 15 Feb 2024 09:19:24 +0100 Subject: [PATCH 09/20] dedup overlapping inner rings --- snap/snap.go | 92 +++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 84 insertions(+), 8 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index b9f1a9a..39d417b 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -2,6 +2,7 @@ package snap import ( "fmt" + "github.com/tobshub/go-sortedmap" "log" "math" "slices" @@ -15,7 +16,6 @@ import ( "github.com/go-spatial/geom" "github.com/pdok/texel/intgeom" "github.com/pdok/texel/tms20" - "github.com/tobshub/go-sortedmap" orderedmap "github.com/wk8/go-ordered-map/v2" "golang.org/x/exp/constraints" "golang.org/x/exp/maps" @@ -92,7 +92,7 @@ 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)) @@ -126,12 +126,13 @@ 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 { - newPolygons[level] = append(newPolygons[level], [][][2]float64{outerRing}) + newOuters[level] = append(newOuters[level], outerRing) } newInners[level] = append(newInners[level], innerRings...) if keepPointsAndLines { @@ -140,8 +141,10 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ } } + newPolygons := make(map[Level][][][][2]float64, len(levels)) for l := range levelMap { - newPolygonsForLevel := matchInnersToPolygons(newPolygons[l], newInners[l], len(polygon) > 1) + newOuters[l], newInners[l] = dedupeAndSortBySizeInnersOuters(newOuters[l], newInners[l]) + newPolygonsForLevel := matchInnersToPolygons(newPolygonsForLevel, newInners[l], len(polygon) > 1) if len(newPolygonsForLevel) > 1 { newPolygons[l] = newPolygonsForLevel } @@ -150,10 +153,84 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ // points and lines at the end, as outer rings for level, pointsAndLines := range newPointsAndLines { for _, pointOrLine := range pointsAndLines { - newPolygons[level] = append(newPolygons[level], [][][2]float64{pointOrLine}) + newOuters[level] = append(newOuters[level], [][][2]float64{pointOrLine}) + } + } + return floatPolygonsToGeomPolygonsForAllLevels(newOuters) +} + +// lang=python +func dedupeAndSortBySizeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]float64, [][][2]float64) { + // ToDo: optimize by deleting rings from allRings on the fly + allRings := append(outers, inners...) + lenOuters := len(outers) + var indexesToDelete []int + for i := 0; i < len(allRings); i++ { + iIsOuter := i < lenOuters + var equalOutersIndexes []int + var equalInnersIndexes []int + if iIsOuter { + equalOutersIndexes = append(equalOutersIndexes, i) + } else { + equalInnersIndexes = append(equalInnersIndexes, i) + } + compareTwoRings: + for j := i + 1; j < len(allRings); j++ { + jIsOuter := j < lenOuters + // check length + iLen := len(allRings[i]) + jLen := len(allRings[j]) + if iLen != jLen { + continue + } + idx := slices.Index(allRings[j], allRings[i][0]) + if idx < 0 { + continue + } + differentWindingOrder := iIsOuter && !jIsOuter + // Check if rings are equal + for k := 0; k < iLen; k++ { + if !differentWindingOrder && allRings[i][k] != allRings[j][idx+k%iLen] { + continue compareTwoRings + } + if differentWindingOrder && allRings[i][k] != allRings[j][idx-k%iLen] { + continue compareTwoRings + } + } + // they are the same! + if jIsOuter { + equalOutersIndexes = append(equalOutersIndexes, j) + } else { + equalInnersIndexes = append(equalInnersIndexes, j) + } + } + // hier toevoegen aan indexesToDelete. mits meer dan 1 herhaling + difference := int(math.Abs(float64(len(equalOutersIndexes)) - float64(len(equalInnersIndexes)))) + if difference == 0 { + indexesToDelete = append(indexesToDelete, equalOutersIndexes[1:]...) + indexesToDelete = append(indexesToDelete, equalInnersIndexes[1:]...) + } + if difference > 0 { + numToDelete := min(len(equalOutersIndexes), len(equalInnersIndexes)) + indexesToDelete = append(indexesToDelete, equalOutersIndexes[0:numToDelete-1]...) + indexesToDelete = append(indexesToDelete, equalInnersIndexes[0:numToDelete-1]...) + } + } + newOuters := make([][][2]float64, 0, lenOuters) + newInners := make([][][2]float64, 0, len(inners)) + for i, outer := range outers { + if slices.Contains(indexesToDelete, i) { + continue + } + newOuters = append(newOuters, outer) + } + for i, inner := range inners { + if slices.Contains(indexesToDelete, i+lenOuters) { + continue } + newInners = append(newInners, inner) } - return floatPolygonsToGeomPolygonsForAllLevels(newPolygons) + return newOuters, newInners } func matchInnersToPolygons(polygons [][][][2]float64, innerRings [][][2]float64, hasInners bool) [][][][2]float64 { @@ -161,7 +238,6 @@ func matchInnersToPolygons(polygons [][][][2]float64, innerRings [][][2]float64, if len(innerRings) == 0 { return polygons } - _, polygons = dedupePolygonsByOuters(polygons) var polyISortedByOuterAreaDesc []int var innersTurnedOuters [][][2]float64 From 5a5d6b7cde3cc204ea053fbd35600d7934969790 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Thu, 15 Feb 2024 10:08:17 +0100 Subject: [PATCH 10/20] fix plumbing around deduped overlapping rings PDOK-16135 --- snap/snap.go | 54 ++++++++++++----------------------------------- snap/snap_test.go | 7 +----- 2 files changed, 14 insertions(+), 47 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 39d417b..502aa49 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -143,8 +143,8 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ newPolygons := make(map[Level][][][][2]float64, len(levels)) for l := range levelMap { - newOuters[l], newInners[l] = dedupeAndSortBySizeInnersOuters(newOuters[l], newInners[l]) - newPolygonsForLevel := matchInnersToPolygons(newPolygonsForLevel, newInners[l], len(polygon) > 1) + newOuters[l], newInners[l] = dedupeInnersOuters(newOuters[l], newInners[l]) + newPolygonsForLevel := matchInnersToPolygons(outersToPolygons(newOuters[l]), newInners[l], len(polygon) > 1) if len(newPolygonsForLevel) > 1 { newPolygons[l] = newPolygonsForLevel } @@ -153,14 +153,22 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ // points and lines at the end, as outer rings for level, pointsAndLines := range newPointsAndLines { for _, pointOrLine := range pointsAndLines { - newOuters[level] = append(newOuters[level], [][][2]float64{pointOrLine}) + newPolygons[level] = append(newPolygons[level], [][][2]float64{pointOrLine}) } } - return floatPolygonsToGeomPolygonsForAllLevels(newOuters) + return floatPolygonsToGeomPolygonsForAllLevels(newPolygons) +} + +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 } // lang=python -func dedupeAndSortBySizeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]float64, [][][2]float64) { +func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]float64, [][][2]float64) { // ToDo: optimize by deleting rings from allRings on the fly allRings := append(outers, inners...) lenOuters := len(outers) @@ -296,42 +304,6 @@ func sortPolyIdxsByOuterAreaDesc(polygons [][][][2]float64) []int { return areas.Keys() } -// helper for matchInnersToPolygons to delete duplicate polygons -// comparing them only by their outers and asserting that a deleted polygon didn't have inner rings appended yet -// yes it's implemented as ~O(n^2), -// but it's expected that the (outer) rings are usually different even from the first point, -// making it still more efficient than using a hashmap of the entire rings -func dedupePolygonsByOuters(polygons [][][][2]float64) (int, [][][][2]float64) { - numPolygons := len(polygons) - numDeleted := 0 - for i := 0; i < numPolygons; i++ { - ring := polygons[i][0] // only check the outer - compareToOther: - for j := i + 1; j < numPolygons; j++ { - ringLen := len(ring) - other := polygons[j][0] - otherLen := len(other) - if ringLen != otherLen { - continue - } - for k := 0; k < min(ringLen, otherLen); k++ { - if ring[k] != other[k] { - continue compareToOther - } - } - // delete - if len(polygons[j]) > 1 { - panicDeletedPolygonHasInnerRing(polygons[j]) - } - polygons = append(polygons[:j], polygons[j+1:]...) - j-- - numPolygons-- - numDeleted++ - } - } - return numDeleted, polygons -} - // 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. diff --git a/snap/snap_test.go b/snap/snap_test.go index 441f01c..da36288 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -1,7 +1,6 @@ package snap import ( - "fmt" "testing" "github.com/go-spatial/geom/encoding/wkt" @@ -622,16 +621,12 @@ func TestSnap_snapPolygon(t *testing.T) { }, 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 {{28.0, 52.0}, {52.0, 52.0}, {52.0, 12.0}, {12.0, 12.0}, {12.0, 52.0}}, // cw {{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}, {12.0, 52.0}, {12.0, 12.0}, {52.0, 12.0}, {52.0, 52.0}, {28.0, 52.0}}, // ccw - {{28.0, 52.0}, {52.0, 52.0}, {52.0, 12.0}, {12.0, 12.0}, {12.0, 52.0}, {28.0, 52.0}}, // cw - // TODO order of separate outer rings OK. does this overlap the deepest nested square? - // TODO deduplicate this, remove totally? }, { {{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 From 26e13d3c0ea8b0c104d365a1182eb2e93b46afb7 Mon Sep 17 00:00:00 2001 From: kad-dirc Date: Thu, 15 Feb 2024 10:12:02 +0100 Subject: [PATCH 11/20] wip: unittest dedup inner outer rings --- snap/snap.go | 1 - snap/snap_test.go | 1 - 2 files changed, 2 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 39d417b..1da1b45 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -159,7 +159,6 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ return floatPolygonsToGeomPolygonsForAllLevels(newOuters) } -// lang=python func dedupeAndSortBySizeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]float64, [][][2]float64) { // ToDo: optimize by deleting rings from allRings on the fly allRings := append(outers, inners...) diff --git a/snap/snap_test.go b/snap/snap_test.go index 441f01c..2d71c68 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -1,7 +1,6 @@ package snap import ( - "fmt" "testing" "github.com/go-spatial/geom/encoding/wkt" From 5099490450fcdbb91867aef684ef096ababbf23b Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Thu, 15 Feb 2024 12:46:31 +0100 Subject: [PATCH 12/20] fix slice high is exclusive PDOK-16135 --- snap/snap.go | 4 ++-- snap/snap_test.go | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 502aa49..da32c13 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -220,8 +220,8 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl } if difference > 0 { numToDelete := min(len(equalOutersIndexes), len(equalInnersIndexes)) - indexesToDelete = append(indexesToDelete, equalOutersIndexes[0:numToDelete-1]...) - indexesToDelete = append(indexesToDelete, equalInnersIndexes[0:numToDelete-1]...) + indexesToDelete = append(indexesToDelete, equalOutersIndexes[0:numToDelete]...) + indexesToDelete = append(indexesToDelete, equalInnersIndexes[0:numToDelete]...) } } newOuters := make([][][2]float64, 0, lenOuters) diff --git a/snap/snap_test.go b/snap/snap_test.go index da36288..c81bf48 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -625,8 +625,8 @@ func TestSnap_snapPolygon(t *testing.T) { 1: { { {{4.0, 124.0}, {4.0, 4.0}, {60.0, 4.0}, {60.0, 124.0}}, // ccw - {{28.0, 52.0}, {52.0, 52.0}, {52.0, 12.0}, {12.0, 12.0}, {12.0, 52.0}}, // cw {{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 From d25422984d185564257ea7240ebd975f8f5d9d09 Mon Sep 17 00:00:00 2001 From: kad-dirc Date: Thu, 15 Feb 2024 12:49:39 +0100 Subject: [PATCH 13/20] wip: unittest dedup inner outer rings --- snap/snap_test.go | 57 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/snap/snap_test.go b/snap/snap_test.go index da36288..680ca77 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -722,6 +722,63 @@ func Test_kmpDeduplicate(t *testing.T) { } } +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: [][][2]float64{{{}}}, + inners: [][][2]float64{{{}}}, + }, + wantOuters: [][][2]float64{{{}}}, + wantInners: [][][2]float64{}, + }, + { + name: "#outer, #inner = 1, 0", + args: args{ + outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + inners: [][][2]float64{{{}}}, + }, + wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + wantInners: [][][2]float64{}, + }, + { + name: "#outer, #inner = 1, 1", + args: args{ + outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + inners: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + }, + wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + wantInners: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + }, + { + name: "#outer, #inner = 2, 1", + args: args{ + outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}, {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + inners: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + }, + wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + wantInners: [][][2]float64{{{}}}, + }, + } + 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) From 0f42c96e35459299025ddfbe8d144a6bab0d5b9b Mon Sep 17 00:00:00 2001 From: kad-dirc Date: Thu, 15 Feb 2024 13:22:23 +0100 Subject: [PATCH 14/20] fix dedupInnerOuters: use only positive modulus values. + first unittest case --- snap/snap.go | 2 +- snap/snap_test.go | 32 +++++++------------------------- 2 files changed, 8 insertions(+), 26 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 4887111..02c6e9d 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -200,7 +200,7 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl if !differentWindingOrder && allRings[i][k] != allRings[j][idx+k%iLen] { continue compareTwoRings } - if differentWindingOrder && allRings[i][k] != allRings[j][idx-k%iLen] { + if differentWindingOrder && allRings[i][k] != allRings[j][(idx+iLen-k)%iLen] { // "+iLen" to ensure the modulus returns a positive number continue compareTwoRings } } diff --git a/snap/snap_test.go b/snap/snap_test.go index 1b2eae1..22d4c11 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -733,41 +733,23 @@ func Test_dedupeInnersOuters(t *testing.T) { wantOuters [][][2]float64 wantInners [][][2]float64 }{ - { - name: "#outer, #inner = 0, 0", - args: args{ - outers: [][][2]float64{{{}}}, - inners: [][][2]float64{{{}}}, - }, - wantOuters: [][][2]float64{{{}}}, - wantInners: [][][2]float64{}, - }, { name: "#outer, #inner = 1, 0", args: args{ - outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, // square, counter clockwise inners: [][][2]float64{{{}}}, }, - wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - wantInners: [][][2]float64{}, + wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + wantInners: [][][2]float64{{{}}}, }, { name: "#outer, #inner = 1, 1", args: args{ - outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - inners: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, // square, counter clockwise + inners: [][][2]float64{{{0, 0}, {0, 1}, {1, 1}, {1, 0}}}, // square, clockwise }, - wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - wantInners: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - }, - { - name: "#outer, #inner = 2, 1", - args: args{ - outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}, {{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - inners: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - }, - wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - wantInners: [][][2]float64{{{}}}, + wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, + wantInners: [][][2]float64{{{0, 0}, {0, 1}, {1, 1}, {1, 0}}}, }, } for _, tt := range tests { From 8c0b7eabfd2226636830308c25aa5aa80bc77df5 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Thu, 15 Feb 2024 14:22:30 +0100 Subject: [PATCH 15/20] skip already deleted rings PDOK-16135 Co-authored-by: kad-dirc --- snap/snap.go | 143 +++++++++++++++++++++++++++++++-------------------- 1 file changed, 88 insertions(+), 55 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index da32c13..0ae4409 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -167,78 +167,90 @@ func outersToPolygons(outers [][][2]float64) [][][][2]float64 { return polygons } -// lang=python func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]float64, [][][2]float64) { // ToDo: optimize by deleting rings from allRings on the fly - allRings := append(outers, inners...) lenOuters := len(outers) - var indexesToDelete []int - for i := 0; i < len(allRings); i++ { - iIsOuter := i < lenOuters - var equalOutersIndexes []int - var equalInnersIndexes []int - if iIsOuter { - equalOutersIndexes = append(equalOutersIndexes, i) - } else { - equalInnersIndexes = append(equalInnersIndexes, i) + lenInners := len(inners) + lenAll := lenOuters + lenInners + indexesToDelete := make(map[int]bool) // true means outer + for i := 0; i < lenAll; i++ { + if _, deleted := indexesToDelete[i]; deleted { + continue } - compareTwoRings: - for j := i + 1; j < len(allRings); j++ { - jIsOuter := j < lenOuters - // check length - iLen := len(allRings[i]) - jLen := len(allRings[j]) - if iLen != jLen { - continue - } - idx := slices.Index(allRings[j], allRings[i][0]) - if idx < 0 { + iIsOuter := i < lenOuters + equalIndexes := make(map[int]bool) // true means outer + equalIndexes[i] = iIsOuter + for j := i + 1; j < lenAll; j++ { + if _, deleted := indexesToDelete[j]; deleted { continue } - differentWindingOrder := iIsOuter && !jIsOuter - // Check if rings are equal - for k := 0; k < iLen; k++ { - if !differentWindingOrder && allRings[i][k] != allRings[j][idx+k%iLen] { - continue compareTwoRings - } - if differentWindingOrder && allRings[i][k] != allRings[j][idx-k%iLen] { - continue compareTwoRings - } + jIsOuter := j < lenOuters + var ringI [][2]float64 + if iIsOuter { + ringI = outers[i] + } else { + ringI = inners[i-lenOuters] } - // they are the same! + var ringJ [][2]float64 if jIsOuter { - equalOutersIndexes = append(equalOutersIndexes, j) + ringJ = outers[j] } else { - equalInnersIndexes = append(equalInnersIndexes, j) + ringJ = inners[j-lenOuters] + } + if !ringsAreEqual(ringI, ringJ, iIsOuter, jIsOuter) { + continue } + equalIndexes[j] = jIsOuter } - // hier toevoegen aan indexesToDelete. mits meer dan 1 herhaling - difference := int(math.Abs(float64(len(equalOutersIndexes)) - float64(len(equalInnersIndexes)))) + lenEqualOuters := countVals(equalIndexes, true) + lenEqualInners := countVals(equalIndexes, false) + difference := int(math.Abs(float64(lenEqualOuters) - float64(lenEqualInners))) + var numOutersToDelete, numInnersToDelete int if difference == 0 { - indexesToDelete = append(indexesToDelete, equalOutersIndexes[1:]...) - indexesToDelete = append(indexesToDelete, equalInnersIndexes[1:]...) - } - if difference > 0 { - numToDelete := min(len(equalOutersIndexes), len(equalInnersIndexes)) - indexesToDelete = append(indexesToDelete, equalOutersIndexes[0:numToDelete]...) - indexesToDelete = append(indexesToDelete, equalInnersIndexes[0:numToDelete]...) + // delete all but one + numOutersToDelete = lenEqualOuters - 1 + numInnersToDelete = lenEqualInners - 1 + } else { // difference > 0 + // delete the surplus + numOutersToDelete = min(lenEqualOuters, lenEqualInners) + numInnersToDelete = numOutersToDelete + } + for equalI, isOuter := range equalIndexes { + if isOuter && numOutersToDelete > 0 { + indexesToDelete[equalI] = isOuter + numOutersToDelete-- + } else if !isOuter && numInnersToDelete > 0 { + indexesToDelete[equalI] = isOuter + numInnersToDelete-- + } } } - newOuters := make([][][2]float64, 0, lenOuters) - newInners := make([][][2]float64, 0, len(inners)) - for i, outer := range outers { - if slices.Contains(indexesToDelete, i) { - continue - } - newOuters = append(newOuters, outer) + 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 } - for i, inner := range inners { - if slices.Contains(indexesToDelete, i+lenOuters) { - continue + 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-k%ringLen] { + return false } - newInners = append(newInners, inner) } - return newOuters, newInners + return true } func matchInnersToPolygons(polygons [][][][2]float64, innerRings [][][2]float64, hasInners bool) [][][][2]float64 { @@ -799,6 +811,27 @@ func lastMatch[T comparable](haystack, needle []T) 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 map[K]V, v V) int { + n := 0 + for k := range m { + if m[k] == 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 { From ee965f2e203cdf9b1ffec3d21c3b0f90b07d94a3 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Thu, 15 Feb 2024 14:43:07 +0100 Subject: [PATCH 16/20] fix plumbing around deduped overlapping rings PDOK-16135 --- snap/snap.go | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/snap/snap.go b/snap/snap.go index c68add2..b1aa4e1 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -145,7 +145,7 @@ func addPointsAndSnap(ix *PointIndex, polygon geom.Polygon, levels []Level) map[ 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) > 1 { + if len(newPolygonsForLevel) > 0 { newPolygons[l] = newPolygonsForLevel } } From b6dc854a27352e6ebc7302646ff86beca00f1840 Mon Sep 17 00:00:00 2001 From: kad-dirc Date: Thu, 15 Feb 2024 14:54:15 +0100 Subject: [PATCH 17/20] dedupInnerOuters add unittest cases --- snap/snap.go | 6 +++ snap/snap_test.go | 134 +++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 132 insertions(+), 8 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 02c6e9d..1781207 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -173,6 +173,9 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl lenOuters := len(outers) var indexesToDelete []int for i := 0; i < len(allRings); i++ { + if slices.Contains(indexesToDelete, i) { + continue + } iIsOuter := i < lenOuters var equalOutersIndexes []int var equalInnersIndexes []int @@ -183,6 +186,9 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl } compareTwoRings: for j := i + 1; j < len(allRings); j++ { + if slices.Contains(indexesToDelete, j) { + continue + } jIsOuter := j < lenOuters // check length iLen := len(allRings[i]) diff --git a/snap/snap_test.go b/snap/snap_test.go index 22d4c11..faa0994 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -733,23 +733,123 @@ func Test_dedupeInnersOuters(t *testing.T) { 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: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, // square, counter clockwise - inners: [][][2]float64{{{}}}, + outers: squareRingArray(1, true), + inners: squareRingArray(0, false), }, - wantOuters: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - wantInners: [][][2]float64{{{}}}, + wantOuters: squareRingArray(1, true), + wantInners: squareRingArray(0, false), }, { name: "#outer, #inner = 1, 1", args: args{ - outers: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, // square, counter clockwise - inners: [][][2]float64{{{0, 0}, {0, 1}, {1, 1}, {1, 0}}}, // square, clockwise + 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: [][][2]float64{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, - wantInners: [][][2]float64{{{0, 0}, {0, 1}, {1, 1}, {1, 0}}}, + 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 { @@ -831,3 +931,21 @@ func wktMustEncode(g geom.Geometry) (s string) { } 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 +} From 3a772a42ac9015f580fa1c5d5e02e8873e4ac335 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Thu, 15 Feb 2024 14:58:47 +0100 Subject: [PATCH 18/20] fix parentheses PDOK-16135 --- snap/snap.go | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/snap/snap.go b/snap/snap.go index b1aa4e1..35a0a26 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -243,7 +243,7 @@ func ringsAreEqual(ringI, ringJ [][2]float64, iIsOuter, jIsOuter bool) bool { differentWindingOrder := iIsOuter && !jIsOuter // Check if rings are equal for k := 0; k < ringLen; k++ { - if !differentWindingOrder && ringI[k] != ringJ[idx+k%ringLen] { + 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 From 2e31eef48ffc9483146ca2cdf4d3230b8f7e0ce0 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Thu, 15 Feb 2024 15:01:42 +0100 Subject: [PATCH 19/20] skip already processed rings PDOK-16135 --- snap/snap.go | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 35a0a26..60efa5f 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -172,16 +172,17 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl lenOuters := len(outers) lenInners := len(inners) lenAll := lenOuters + lenInners - indexesToDelete := make(map[int]bool) // true means outer + processedIndexes := make(map[int]bool) // true means outer + indexesToDelete := make(map[int]bool) // true means outer for i := 0; i < lenAll; i++ { - if _, deleted := indexesToDelete[i]; deleted { + if _, processed := processedIndexes[i]; processed { continue } iIsOuter := i < lenOuters equalIndexes := make(map[int]bool) // true means outer equalIndexes[i] = iIsOuter for j := i + 1; j < lenAll; j++ { - if _, deleted := indexesToDelete[j]; deleted { + if _, processed := processedIndexes[j]; processed { continue } jIsOuter := j < lenOuters @@ -202,6 +203,11 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl } equalIndexes[j] = jIsOuter } + if len(equalIndexes) <= 1 { + continue + } + maps.Copy(processedIndexes, equalIndexes) + lenEqualOuters := countVals(equalIndexes, true) lenEqualInners := countVals(equalIndexes, false) difference := int(math.Abs(float64(lenEqualOuters) - float64(lenEqualInners))) @@ -225,6 +231,10 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl } } } + + if len(indexesToDelete) == 0 { + return outers, inners + } newOuters := deleteFromSliceByIndex(outers, indexesToDelete, 0) newInners := deleteFromSliceByIndex(inners, indexesToDelete, lenOuters) return newOuters, newInners From 1993bbb9fa1f9aadc1b73a211db45e9cceec4553 Mon Sep 17 00:00:00 2001 From: Roel Arents Date: Thu, 15 Feb 2024 23:06:47 +0100 Subject: [PATCH 20/20] fix map ranging in order in dedupeInnersOuters --- snap/snap.go | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 60efa5f..a9c7a47 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -26,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. // @@ -172,15 +174,15 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl lenOuters := len(outers) lenInners := len(inners) lenAll := lenOuters + lenInners - processedIndexes := make(map[int]bool) // true means outer - indexesToDelete := make(map[int]bool) // true means outer + 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 := make(map[int]bool) // true means outer - equalIndexes[i] = iIsOuter + 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 @@ -201,12 +203,11 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl if !ringsAreEqual(ringI, ringJ, iIsOuter, jIsOuter) { continue } - equalIndexes[j] = jIsOuter + equalIndexes.Set(j, jIsOuter) } - if len(equalIndexes) <= 1 { + if equalIndexes.Len() <= 1 { continue } - maps.Copy(processedIndexes, equalIndexes) lenEqualOuters := countVals(equalIndexes, true) lenEqualInners := countVals(equalIndexes, false) @@ -221,7 +222,9 @@ func dedupeInnersOuters(outers [][][2]float64, inners [][][2]float64) ([][][2]fl numOutersToDelete = min(lenEqualOuters, lenEqualInners) numInnersToDelete = numOutersToDelete } - for equalI, isOuter := range equalIndexes { + 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-- @@ -832,10 +835,10 @@ func deleteFromSliceByIndex[V any, X any](s []V, indexesToDelete map[int]X, inde return r } -func countVals[K, V comparable](m map[K]V, v V) int { +func countVals[K, V comparable](m *orderedmap.OrderedMap[K, V], v V) int { n := 0 - for k := range m { - if m[k] == v { + for p := m.Oldest(); p != nil; p = p.Next() { + if p.Value == v { n++ } }