From 668706ed7123d014a315d7af893bbf3db269d554 Mon Sep 17 00:00:00 2001 From: Michiel Korpel Date: Tue, 14 Nov 2023 09:21:27 +0100 Subject: [PATCH 1/7] Implement winding order check using vector angles --- snap/snap.go | 103 +++++++++++++++-------------------------------- snap/vector2d.go | 29 +++++++++++++ 2 files changed, 61 insertions(+), 71 deletions(-) create mode 100644 snap/vector2d.go diff --git a/snap/snap.go b/snap/snap.go index 4d4cb8b..1f44c0a 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -2,9 +2,9 @@ package snap import ( "fmt" + "log" "math" "slices" - "strings" "github.com/go-spatial/geom" "github.com/pdok/texel/processing" @@ -68,6 +68,7 @@ func addPointsAndSnap(ix *PointIndex, polygon *geom.Polygon) *geom.Polygon { // inner rings (ringIdx != 0) should be clockwise shouldBeClockwise := ringIdx != 0 // winding order is reversed if incorrect + // ensureCorrectWindingOrder(deduplicatedRing, shouldBeClockwise) ensureCorrectWindingOrder(deduplicatedRing, shouldBeClockwise) newPolygon = append(newPolygon, deduplicatedRing) } @@ -78,86 +79,46 @@ func addPointsAndSnap(ix *PointIndex, polygon *geom.Polygon) *geom.Polygon { // validate winding order (CCW for outer rings, CW for inner rings) // if winding order is incorrect, ring is reversed to correct winding order func ensureCorrectWindingOrder(ring [][2]float64, shouldBeClockwise bool) { - steps := []string{} + // 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 + } stepCounts := map[bool]int{true: 0, false: 0} for i := 0; i < len(ring); i++ { - steps = append(steps, directionOfStep(ring[i], ring[(i+1)%len(ring)])) - // need at least a pair of steps to check winding order - if len(steps) <= 1 { - continue - } - // steps back (e.g., 'R,L') or repeated steps (e.g., 'R,R') count as a step for the winding order the ring should have - if steps[len(steps)-2] == oppositeDirectionOf(steps[len(steps)-1]) || steps[len(steps)-2] == steps[len(steps)-1] { - stepCounts[shouldBeClockwise]++ - } else { - stepCounts[isClockwise(steps[len(steps)-2:])]++ - } + // check angle between the vector from the previous point to the current point and the vector from the current point to the next point + points := [3][2]float64{ring[mod(i-1, len(ring))], ring[i], ring[mod(i+1, len(ring))]} + isClockwise := isClockwise(points, shouldBeClockwise) + // TODO: implement penalty for 'bites' (parts of the ring that fall within the ring's extent)? + stepCounts[isClockwise]++ + } + var msg string + if shouldBeClockwise { + msg = "inner ring should be clockwise;" + } else { + msg = "outer ring should be counterclockwise;" } if stepCounts[shouldBeClockwise] < stepCounts[!shouldBeClockwise] { + msg += fmt.Sprintf(" incorrect winding order (correct steps: %d, incorrect steps: %d)", stepCounts[shouldBeClockwise], stepCounts[!shouldBeClockwise]) slices.Reverse(ring) + } else { + msg += fmt.Sprintf(" assuming correct winding order (correct steps: %d, incorrect steps: %d)", stepCounts[shouldBeClockwise], stepCounts[!shouldBeClockwise]) } + log.Printf("%s final ring: %v\n", msg, ring) } -// returns if pair of steps is clockwise -func isClockwise(stepsPair []string) bool { - clockwiseSteps := []string{ - "U,R", // up, then right - "U,RU", // up, then right+up - "U,RD", // up, then right+down - "RU,R", // right+up, then right - "RU,RD", // right+up, then right+down - "RU,D", // right+up, then down - "R,D", // right, then down - "R,RD", // right, then right+down - "R,LD", // right, then left+down - "RD,D", // right+down, then down - "RD,LD", // right+down, then left+down - "RD,L", // right+down, then left - "D,L", // down, then left - "D,LU", // down, then left+up - "D,LD", // down, then left+down - "LD,L", // left+down, then left - "LD,LU", // left+down, then left+up - "LD,U", // left+down, then up - "L,U", // left, then up - "L,LU", // left, then left+up - "L,RU", // left, then right+up - "LU,R", // left+up, then right - "LU,RU", // left+up, then right+up - "LU,U", // left+up, then up - } - return slices.Contains(clockwiseSteps, strings.Join(stepsPair, ",")) -} - -// returns opposite of a direction -func oppositeDirectionOf(direction string) string { - opposites := map[string]string{ - "U": "D", - "D": "U", - "R": "L", - "L": "R", - "RU": "LD", - "RD": "LU", - "LU": "RD", - "LD": "RU", - } - return opposites[direction] -} - -// determines direction of a step from one point to another -func directionOfStep(from, to [2]float64) string { - direction := "" - if to[0] > from[0] { - direction += "R" - } else if to[0] < from[0] { - direction += "L" +// 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, shouldBeClockwise bool) 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) } - if to[1] > from[1] { - direction += "U" - } else if to[1] < from[1] { - direction += "D" + 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 math.Round(relativeAngle) == 0.0 || math.Round(relativeAngle) == 180.0 { + return shouldBeClockwise } - return direction + return math.Round(mod((vector2.angle()-relativeAngle), 360)) != math.Round(vector1.angle()) } // deduplication using an implementation of the Knuth-Morris-Pratt algorithm diff --git a/snap/vector2d.go b/snap/vector2d.go new file mode 100644 index 0000000..f946310 --- /dev/null +++ b/snap/vector2d.go @@ -0,0 +1,29 @@ +package snap + +import "math" + +// vector2d represents the vector between two points in 2D space +type vector2d struct { + x float64 + y float64 +} + +// returns the angle of a vector in relation to the x-axis in degrees (normalised to range 0-360) +func (vec vector2d) angle() float64 { + angleRad := math.Atan2(vec.y, vec.x) + if angleRad < 0 { + angleRad += (2 * math.Pi) + } + angleDeg := angleRad * (180 / math.Pi) + return angleDeg +} + +// dot product of vector with another vector +func (vec vector2d) dot(otherVec vector2d) float64 { + return (vec.x * otherVec.x) + (vec.y * otherVec.y) +} + +// magnitude of vector +func (vec vector2d) magnitude() float64 { + return math.Sqrt(math.Pow(vec.x, 2) + math.Pow(vec.y, 2)) +} From f3d3ca40bc149fa6edebacf3f18635c0957049a5 Mon Sep 17 00:00:00 2001 From: Michiel Korpel Date: Tue, 14 Nov 2023 12:07:10 +0100 Subject: [PATCH 2/7] Determine winding order by examining direction of vectors around the rightmost lowest point --- snap/snap.go | 44 ++++++++++++++++++-------------------------- 1 file changed, 18 insertions(+), 26 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 1f44c0a..3ede1cd 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -2,7 +2,6 @@ package snap import ( "fmt" - "log" "math" "slices" @@ -64,12 +63,11 @@ func addPointsAndSnap(ix *PointIndex, polygon *geom.Polygon) *geom.Polygon { } default: // deduplicate points in the ring, check winding order, then add to the polygon - deduplicatedRing := kmpDeduplicate(newRing) + deduplicatedRing, rightmostLowestPoint := kmpDeduplicate(newRing) // inner rings (ringIdx != 0) should be clockwise shouldBeClockwise := ringIdx != 0 // winding order is reversed if incorrect - // ensureCorrectWindingOrder(deduplicatedRing, shouldBeClockwise) - ensureCorrectWindingOrder(deduplicatedRing, shouldBeClockwise) + ensureCorrectWindingOrder(deduplicatedRing, rightmostLowestPoint, shouldBeClockwise) newPolygon = append(newPolygon, deduplicatedRing) } } @@ -78,32 +76,18 @@ func addPointsAndSnap(ix *PointIndex, polygon *geom.Polygon) *geom.Polygon { // validate winding order (CCW for outer rings, CW for inner rings) // if winding order is incorrect, ring is reversed to correct winding order -func ensureCorrectWindingOrder(ring [][2]float64, shouldBeClockwise bool) { +func ensureCorrectWindingOrder(ring [][2]float64, rightmostLowestPoint [2]float64, shouldBeClockwise 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 } - stepCounts := map[bool]int{true: 0, false: 0} - for i := 0; i < len(ring); i++ { - // check angle between the vector from the previous point to the current point and the vector from the current point to the next point - points := [3][2]float64{ring[mod(i-1, len(ring))], ring[i], ring[mod(i+1, len(ring))]} - isClockwise := isClockwise(points, shouldBeClockwise) - // TODO: implement penalty for 'bites' (parts of the ring that fall within the ring's extent)? - stepCounts[isClockwise]++ - } - var msg string - if shouldBeClockwise { - msg = "inner ring should be clockwise;" - } else { - msg = "outer ring should be counterclockwise;" - } - if stepCounts[shouldBeClockwise] < stepCounts[!shouldBeClockwise] { - msg += fmt.Sprintf(" incorrect winding order (correct steps: %d, incorrect steps: %d)", stepCounts[shouldBeClockwise], stepCounts[!shouldBeClockwise]) + i := slices.Index(ring, rightmostLowestPoint) + // 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))]} + isClockwise := isClockwise(points, shouldBeClockwise) + if isClockwise != shouldBeClockwise { slices.Reverse(ring) - } else { - msg += fmt.Sprintf(" assuming correct winding order (correct steps: %d, incorrect steps: %d)", stepCounts[shouldBeClockwise], stepCounts[!shouldBeClockwise]) } - log.Printf("%s final ring: %v\n", msg, ring) } // determines whether a pair of vectors turns clockwise by examining the relationship of their relative angle to the angles of the vectors @@ -124,7 +108,8 @@ func isClockwise(points [3][2]float64, shouldBeClockwise bool) bool { // deduplication using an implementation of the Knuth-Morris-Pratt algorithm // //nolint:cyclop,funlen -func kmpDeduplicate(newRing [][2]float64) [][2]float64 { +func kmpDeduplicate(newRing [][2]float64) ([][2]float64, [2]float64) { + rightmostLowestPoint := [2]float64{math.MinInt, math.MaxInt} deduplicatedRing := make([][2]float64, len(newRing)) copy(deduplicatedRing, newRing) // map of indices to remove, sorted by starting index of each sequence to remove @@ -135,6 +120,13 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { visitedPoints := [][2]float64{} for i := 0; i < len(newRing); { vertex := newRing[i] + // check if point 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] + } else if vertex[1] == rightmostLowestPoint[1] && vertex[0] > rightmostLowestPoint[0] { + rightmostLowestPoint[0] = vertex[0] + } // not a step back, continue if len(visitedPoints) <= 1 || visitedPoints[len(visitedPoints)-2] != vertex { visitedPoints = append(visitedPoints, vertex) @@ -275,7 +267,7 @@ func kmpDeduplicate(newRing [][2]float64) [][2]float64 { deduplicatedRing = append(deduplicatedRing[:sequence[0]-offset], deduplicatedRing[sequence[1]-offset:]...) offset += sequence[1] - sequence[0] } - return deduplicatedRing + return deduplicatedRing, rightmostLowestPoint } // repeatedly calls kmpSearch, returning all starting indexes of 'find' in 'corpus' From 184fd1f072e9ef88a5904f750439e560c6aff507 Mon Sep 17 00:00:00 2001 From: Michiel Korpel Date: Tue, 14 Nov 2023 13:58:35 +0100 Subject: [PATCH 3/7] Correct winding order of expected polygon --- snap/snap_test.go | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/snap/snap_test.go b/snap/snap_test.go index 444ed38..830f803 100644 --- a/snap/snap_test.go +++ b/snap/snap_test.go @@ -49,13 +49,13 @@ func Test_snapPolygon(t *testing.T) { {110892.93099999999685679, 504407.8219999999855645}, }}, want: &geom.Polygon{{ - {110899.1928125, 504431.1559375}, - {110906.8709375, 504428.8065625}, // horizontal line still here - {110907.6453125, 504428.8065625}, // horizontal line still here - {110909.4565625, 504436.9046875}, - {110920.0353125, 504436.0778125}, - {110929.4459375, 504407.8196875}, {110892.9321875, 504407.8196875}, + {110929.4459375, 504407.8196875}, + {110920.0353125, 504436.0778125}, + {110909.4565625, 504436.9046875}, + {110907.6453125, 504428.8065625}, // horizontal line still here + {110906.8709375, 504428.8065625}, // horizontal line still here + {110899.1928125, 504431.1559375}, }}, }, { From e8e260be150474d92a92cd3dc6e7420f6b7d947f Mon Sep 17 00:00:00 2001 From: Michiel Korpel Date: Tue, 14 Nov 2023 16:10:37 +0100 Subject: [PATCH 4/7] Minor refactor --- snap/snap.go | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/snap/snap.go b/snap/snap.go index 3ede1cd..2d58420 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -84,8 +84,7 @@ func ensureCorrectWindingOrder(ring [][2]float64, rightmostLowestPoint [2]float6 i := slices.Index(ring, rightmostLowestPoint) // 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))]} - isClockwise := isClockwise(points, shouldBeClockwise) - if isClockwise != shouldBeClockwise { + if isClockwise(points, shouldBeClockwise) != shouldBeClockwise { slices.Reverse(ring) } } @@ -120,7 +119,7 @@ func kmpDeduplicate(newRing [][2]float64) ([][2]float64, [2]float64) { visitedPoints := [][2]float64{} for i := 0; i < len(newRing); { vertex := newRing[i] - // check if point is rightmost lowest point (either y is lower, or y is equal and x is higher) + // 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] From b920b40ca167ede24354801eab32aa90b17b90c2 Mon Sep 17 00:00:00 2001 From: Michiel Korpel Date: Tue, 14 Nov 2023 16:19:20 +0100 Subject: [PATCH 5/7] Minor refactor --- snap/snap.go | 1 + 1 file changed, 1 insertion(+) diff --git a/snap/snap.go b/snap/snap.go index 2d58420..650067c 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -105,6 +105,7 @@ func isClockwise(points [3][2]float64, shouldBeClockwise bool) bool { } // deduplication using an implementation of the Knuth-Morris-Pratt algorithm +// function also determines the rightmost lowest point of the ring, which is used to ensure correct winding order // //nolint:cyclop,funlen func kmpDeduplicate(newRing [][2]float64) ([][2]float64, [2]float64) { From 1d81153945dfcecc3ed4987148f1bbc00946d95a Mon Sep 17 00:00:00 2001 From: Michiel Korpel Date: Tue, 14 Nov 2023 16:34:50 +0100 Subject: [PATCH 6/7] Upgrade golangci-lint version --- .github/workflows/lint-go.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/lint-go.yml b/.github/workflows/lint-go.yml index e773b9c..443b1df 100644 --- a/.github/workflows/lint-go.yml +++ b/.github/workflows/lint-go.yml @@ -17,7 +17,7 @@ jobs: - uses: actions/setup-go@v4 with: - go-version: '1.20' + go-version: '1.21' cache: false - uses: actions/checkout@v3 @@ -26,7 +26,7 @@ jobs: uses: golangci/golangci-lint-action@v3 with: # Optional: version of golangci-lint to use in form of v1.2 or v1.2.3 or `latest` to use the latest version - version: v1.52 + version: v1.54 # Optional: working directory, useful for monorepos # working-directory: somedir From dd0b36614ca1bd82cbe9d7bcc88762b241ae5278 Mon Sep 17 00:00:00 2001 From: Michiel Korpel Date: Tue, 14 Nov 2023 16:35:22 +0100 Subject: [PATCH 7/7] Correct go version in Dockerfile --- Dockerfile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Dockerfile b/Dockerfile index 14debc7..21f6aa6 100644 --- a/Dockerfile +++ b/Dockerfile @@ -26,7 +26,7 @@ RUN go test ./... -covermode=atomic RUN go build -v -ldflags='-s -w -linkmode auto' -a -installsuffix cgo -o /texel . # FROM scratch -FROM golang:1.18-bullseye +FROM golang:1.21-bullseye RUN apt-get update && apt-get install -y \ libsqlite3-mod-spatialite \ && rm -rf /var/lib/apt/lists/* @@ -39,4 +39,4 @@ WORKDIR / # Import from builder. COPY --from=build-env /texel / -CMD ["/texel"] \ No newline at end of file +CMD ["/texel"]