diff --git a/pointindex/pointindex.go b/pointindex/pointindex.go index 52f7f39..ed82742 100644 --- a/pointindex/pointindex.go +++ b/pointindex/pointindex.go @@ -27,9 +27,9 @@ const ( xAx = 0 yAx = 1 // left = 0b00 - right = 0b01 + // right = 0b01 // bottom = 0b00 - top = 0b10 + // top = 0b10 // bottomleft = bottom | left // 0b00 // bottomright = bottom | right // 0b01 // topleft = top | left // 0b10 @@ -71,10 +71,26 @@ type PointIndex struct { quadrants map[Level]map[morton.Z]Quadrant hitOnce map[Level]map[intgeom.Point][]int hitMultiple map[Level]map[intgeom.Point][]int + + intersect map[Level]map[morton.Z]TileData + registeredSegments map[morton.Z][]SegmentIdx } -type Level = uint -type Q = int // quadrant index (0, 1, 2 or 3) +type SegmentIdx struct { + ringIdx int + pointIdx int +} + +type TileData struct { + isOutside bool + isInside bool + isIntersect bool +} + +type ( + Level = uint + Q = int // quadrant index (0, 1, 2 or 3) +) func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID) (*PointIndex, error) { // assuming IsQuadTree was tested before @@ -212,84 +228,6 @@ func (ix *PointIndex) GetHitMultiple(l Level) map[intgeom.Point][]int { return ix.hitMultiple[l] } -// insertCoord adds a point into this pc, assuming the point is inside its extent -func (ix *PointIndex) insertCoord(deepestX int, deepestY int) { - var l Level - for l = 0; l <= ix.deepestLevel; l++ { //nolint:intrange - //nolint:gosec // G115 - x := uint(deepestX) / mathhelp.Pow2(ix.deepestLevel-l) - //nolint:gosec // G115 - y := uint(deepestY) / mathhelp.Pow2(ix.deepestLevel-l) - z := morton.MustToZ(x, y) - if ix.quadrants[l] == nil { // probably already initialized by InsertPolygon - ix.quadrants[l] = make(map[morton.Z]Quadrant) - } - extent, centroid := ix.getQuadrantExtentAndCentroid(l, x, y, ix.intExtent) - ix.quadrants[l][z] = Quadrant{ - z: z, - intExtent: extent, - intCentroid: centroid, - } - } -} - -func (ix *PointIndex) getQuadrantExtentAndCentroid(level Level, x, y uint, intRootExtent intgeom.Extent) (intgeom.Extent, intgeom.Point) { - //nolint:gosec // G115 - intQuadrantSpan := int64(mathhelp.Pow2(ix.deepestLevel-level)) * ix.deepestRes - intMinX := intRootExtent.MinX() - intMinY := intRootExtent.MinY() - intExtent := intgeom.Extent{ - //nolint:gosec // G115 - intMinX + int64(x)*intQuadrantSpan, // minx - //nolint:gosec // G115 - intMinY + int64(y)*intQuadrantSpan, // miny - //nolint:gosec // G115 - intMinX + int64(x+1)*intQuadrantSpan, // maxx - //nolint:gosec // G115 - intMinY + int64(y+1)*intQuadrantSpan, // maxy - } - intCentroid := intgeom.Point{ - //nolint:gosec // G115 - intMinX + (int64(x))*intQuadrantSpan + intQuadrantSpan/2, // <-- here is the plus 0.5 internal pixel size - //nolint:gosec // G115 - intMinY + (int64(y))*intQuadrantSpan + intQuadrantSpan/2, - } - return intExtent, intCentroid -} - -func (ix *PointIndex) snapClosestPoints(intLine intgeom.Line, levelMap map[Level]any) map[Level][]Quadrant { - if len(levelMap) == 0 || !lineIntersects(intLine, ix.intExtent) { - return nil - } - quadrantsIntersectedPerLevel := make(map[Level][]Quadrant, len(levelMap)) - parents := []Quadrant{ix.Quadrant} - if _, includeLevelZero := levelMap[0]; includeLevelZero { - quadrantsIntersectedPerLevel[0] = parents - } - - var level Level - for level = 1; level <= ix.deepestLevel; level++ { - quadrantsIntersected := make([]Quadrant, 0, 10) // TODO good estimate of expected count, based on line length / quadrant span * geom points? - for _, parent := range parents { - quadrantZs := getQuadrantZs(parent.z) - quadrantsWithPoints := make(map[Q]Quadrant, 4) - for q, quadrantZ := range quadrantZs { - if quadrant, exists := ix.quadrants[level][quadrantZ]; exists { - quadrantsWithPoints[q] = quadrant - } - } - for _, q := range findIntersectingQuadrants(intLine, quadrantsWithPoints, parent) { - quadrantsIntersected = append(quadrantsIntersected, quadrantsWithPoints[q]) - } - } - parents = quadrantsIntersected - if _, isLevelIncluded := levelMap[level]; isLevelIncluded { - quadrantsIntersectedPerLevel[level] = quadrantsIntersected - } - } - return quadrantsIntersectedPerLevel -} - //nolint:cyclop func findIntersectingQuadrants(intLine intgeom.Line, quadrants map[Q]Quadrant, parent Quadrant) []Q { pt1InfiniteQuadrantI := getInfiniteQuadrant(intLine[0], parent.intCentroid) @@ -375,18 +313,13 @@ func findIntersectingQuadrants(intLine intgeom.Line, quadrants map[Q]Quadrant, p } func getQuadrantZs(parentZ morton.Z) [4]morton.Z { - parentX, parentY := morton.FromZ(parentZ) - quadrantZs := [4]morton.Z{} - for i := range 4 { - //nolint:gosec // G115 - x := parentX*2 + uint(oneIfRight(i)) - //nolint:gosec // G115 - y := parentY*2 + uint(oneIfTop(i)) - z := morton.MustToZ(x, y) - //nolint:gosec // G602 - quadrantZs[i] = z + baseChildZ := parentZ << 2 + return [4]morton.Z{ + baseChildZ, + baseChildZ + 1, + baseChildZ + 2, + baseChildZ + 3, } - return quadrantZs } // containsPoint checks whether a point is contained in a quadrant's extent. @@ -519,12 +452,371 @@ func lineOverlapsInclusiveEdge(intLine intgeom.Line, edgeI int, intEdge intgeom. return lOrd1 != lOrd2 && (mathhelp.IBetweenInc(lOrd1, eOrd1, eOrd2) && intLine[0] != exclusiveTip || mathhelp.IBetweenInc(lOrd2, eOrd1, eOrd2) && intLine[1] != exclusiveTip) } -func oneIfRight(quadrantI int) int { - return quadrantI & right +// A version of Amanatides-Woo +// Concept: parametrize line with parameter t, with t = 0 being the start and t = 1 the endpoint. Then compute the t for which tile boundaries are intersected. By taking a suitable multiple of t, this is guaranteed to be integral. Use these to raycast along the tiles. +func (ix *PointIndex) AmanatidesWoo(line geom.Line, l Level, ringIdx int, pointIdx int) { + intLine := intgeom.FromGeomLine(line) + segmentIdx := SegmentIdx{ + ringIdx: ringIdx, + pointIdx: pointIdx, + } + + p1, p2 := intLine.Point1(), intLine.Point2() + + // Find start and end square + levelDiff := ix.deepestLevel - l + xStart, yStart := ix.findTileCoords(*p1, levelDiff) + xEnd, yEnd := ix.findTileCoords(*p2, levelDiff) + + tileSize := ix.deepestRes << levelDiff + + var xLine, yLine, xDir, yDir, xCoord, yCoord intgeom.M + + // Current location + xCoord, yCoord = xStart, yStart + + // Determine direction of `line` and determine which edges of the starting tile are intersected first. + if p1.Y() > p2.Y() { + yLine = yStart * tileSize + yDir = -1 + } else { + yLine = (yStart + 1) * tileSize + yDir = 1 + } + + if p1.X() > p2.X() { + xLine = xStart * tileSize + xDir = -1 + } else { + xLine = (xStart + 1) * tileSize + xDir = 1 + } + + xVal := (xLine - p1.X()) * (p2.Y() - p1.Y()) * xDir * yDir + yVal := (yLine - p1.Y()) * (p2.X() - p1.X()) * xDir * yDir + + xDelta := yDir * (p2.Y() - p1.Y()) * tileSize + yDelta := xDir * (p2.X() - p1.X()) * tileSize + + // Register start tile + ix.registerQuadrantAndNeighbours(xCoord, yCoord, l, segmentIdx) + + // Main traversal loop + for xCoord != xEnd || yCoord != yEnd { + // The first two cases are for when you hit an edge. + // The last three cases are when you hit a corner. + switch { + case xVal < yVal: + xCoord += xDir + xVal += xDelta + case yVal < xVal: + yCoord += yDir + yVal += yDelta + case xDir == 1 && yDir == -1: + xCoord += xDir + xVal += xDelta + case xDir == -1 && yDir == 1: + yCoord += yDir + yVal += yDelta + default: + xCoord += xDir + yCoord += yDir + xVal += xDelta + yVal += yDelta + } + ix.registerQuadrantAndNeighbours(xCoord, yCoord, l, segmentIdx) + } +} + +// insertCoord adds a point into this pc, assuming the point is inside its extent +func (ix *PointIndex) insertCoord(deepestX int, deepestY int) { + var l Level + for l = 0; l <= ix.deepestLevel; l++ { //nolint:intrange + //nolint:gosec // G115 + x := uint(deepestX) / mathhelp.Pow2(ix.deepestLevel-l) + //nolint:gosec // G115 + y := uint(deepestY) / mathhelp.Pow2(ix.deepestLevel-l) + z := morton.MustToZ(x, y) + if ix.quadrants[l] == nil { // probably already initialized by InsertPolygon + ix.quadrants[l] = make(map[morton.Z]Quadrant) + } + extent, centroid := ix.getQuadrantExtentAndCentroid(l, x, y, ix.intExtent) + ix.quadrants[l][z] = Quadrant{ + z: z, + intExtent: extent, + intCentroid: centroid, + } + } } -func oneIfTop(quadrantI int) int { - return (quadrantI & top) >> 1 +func (ix *PointIndex) getQuadrantExtentAndCentroid(level Level, x, y uint, intRootExtent intgeom.Extent) (intgeom.Extent, intgeom.Point) { + //nolint:gosec // G115 + intQuadrantSpan := int64(mathhelp.Pow2(ix.deepestLevel-level)) * ix.deepestRes + intMinX := intRootExtent.MinX() + intMinY := intRootExtent.MinY() + intExtent := intgeom.Extent{ + //nolint:gosec // G115 + intMinX + int64(x)*intQuadrantSpan, // minx + //nolint:gosec // G115 + intMinY + int64(y)*intQuadrantSpan, // miny + //nolint:gosec // G115 + intMinX + int64(x+1)*intQuadrantSpan, // maxx + //nolint:gosec // G115 + intMinY + int64(y+1)*intQuadrantSpan, // maxy + } + intCentroid := intgeom.Point{ + //nolint:gosec // G115 + intMinX + (int64(x))*intQuadrantSpan + intQuadrantSpan/2, // <-- here is the plus 0.5 internal pixel size + //nolint:gosec // G115 + intMinY + (int64(y))*intQuadrantSpan + intQuadrantSpan/2, + } + return intExtent, intCentroid +} + +func (ix *PointIndex) snapClosestPoints(intLine intgeom.Line, levelMap map[Level]any) map[Level][]Quadrant { + if len(levelMap) == 0 || !lineIntersects(intLine, ix.intExtent) { + return nil + } + quadrantsIntersectedPerLevel := make(map[Level][]Quadrant, len(levelMap)) + parents := []Quadrant{ix.Quadrant} + if _, includeLevelZero := levelMap[0]; includeLevelZero { + quadrantsIntersectedPerLevel[0] = parents + } + + var level Level + for level = 1; level <= ix.deepestLevel; level++ { + quadrantsIntersected := make([]Quadrant, 0, 10) // TODO good estimate of expected count, based on line length / quadrant span * geom points? + for _, parent := range parents { + quadrantZs := getQuadrantZs(parent.z) + quadrantsWithPoints := make(map[Q]Quadrant, 4) + for q, quadrantZ := range quadrantZs { + if quadrant, exists := ix.quadrants[level][quadrantZ]; exists { + quadrantsWithPoints[q] = quadrant + } + } + for _, q := range findIntersectingQuadrants(intLine, quadrantsWithPoints, parent) { + quadrantsIntersected = append(quadrantsIntersected, quadrantsWithPoints[q]) + } + } + parents = quadrantsIntersected + if _, isLevelIncluded := levelMap[level]; isLevelIncluded { + quadrantsIntersectedPerLevel[level] = quadrantsIntersected + } + } + return quadrantsIntersectedPerLevel +} + +func (ix *PointIndex) findTileCoords(intPoint intgeom.Point, levelDiff Level) (intgeom.M, intgeom.M) { + // I don't check for out of bounds here; this should not occur. + deepestX := ((intPoint.X() - ix.intExtent.MinX()) / ix.deepestRes) >> levelDiff + deepestY := ((intPoint.Y() - ix.intExtent.MinY()) / ix.deepestRes) >> levelDiff + return deepestX, deepestY +} + +func (ix *PointIndex) registerQuadrantAndNeighbours(x intgeom.M, y intgeom.M, l Level, idx SegmentIdx) { + steps := [3]intgeom.M{-1, 0, 1} + levelDiff := ix.deepestLevel - l + levelSize := ix.deepestSize >> levelDiff + + var xNeighbour, yNeighbour intgeom.M + + for _, xDiff := range steps { + for _, yDiff := range steps { + xNeighbour = x + xDiff + yNeighbour = y + yDiff + if xNeighbour < 0 || xNeighbour > intgeom.M(levelSize)-1 || //nolint:gosec + yNeighbour < 0 || yNeighbour > intgeom.M(levelSize)-1 { //nolint:gosec + continue // Out of bounds + } + ix.registerQuadrant(xNeighbour, yNeighbour, l, idx) + } + } +} + +func (ix *PointIndex) registerQuadrant(x intgeom.M, y intgeom.M, l Level, idx SegmentIdx) { + z := morton.MustToZ(uint(x), uint(y)) //nolint:gosec // Overflows checked by MustToZ + ix.registerSegment(z, idx) + + // Create maps if nonexistent + if ix.intersect == nil { + ix.intersect = make(map[Level]map[morton.Z]TileData, l) + for i := uint(0); i <= l; i++ { + ix.intersect[i] = make(map[morton.Z]TileData) + } + + } + // Insert at each level + for i := uint(0); i <= l; i++ { + levelMap := ix.intersect[l-i] + _, b := levelMap[z] + + // Create if does not exist + if !b { + levelMap[z] = TileData{isIntersect: true} + } else { + return + } + + // Update Morton coordinate + z >>= 2 + } +} + +func (ix *PointIndex) registerSegment(z morton.Z, idx SegmentIdx) { + if ix.registeredSegments == nil { + ix.registeredSegments = make(map[morton.Z][]SegmentIdx) + } + + segments, present := ix.registeredSegments[z] + + if !present { + ix.registeredSegments[z] = []SegmentIdx{idx} + return + } + + lastIdx := segments[len(segments)-1] + + if idx != lastIdx { + ix.registeredSegments[z] = append(segments, idx) + } +} + +func (ix *PointIndex) registerPolygonEdges(polygon geom.Polygon, targetLevel Level) { + rings := polygon.LinearRings() + + for ringidx, ring := range rings { + for vertexidx, vertex := range ring { + nextIdx := (vertexidx + 1) % len(ring) + nextVertex := ring[nextIdx] + line := geom.Line{vertex, nextVertex} + ix.AmanatidesWoo(line, targetLevel, ringidx, vertexidx) + } + } +} + +// Going level-by-level in the tree, finding tiles left of origin marked as intersecting +func (ix *PointIndex) findIntersectingTilesLeft(xOrigin uint, yOrigin uint, targetLevel Level) []morton.Z { + todoListAtCurrentLevel := []morton.Z{0} + var todoListAtNextLevel []morton.Z + var leftChild, rightChild morton.Z + for currentLevel := range targetLevel { + todoListAtNextLevel = make([]morton.Z, 0) // currentLen is relatively arbitrary (check if works) + + xOriginAtNextLevel := xOrigin >> (targetLevel - currentLevel - 1) + yOriginAtNextLevel := yOrigin >> (targetLevel - currentLevel - 1) + + nextLevelDown := yOriginAtNextLevel%2 == 0 + for _, currentZ := range todoListAtCurrentLevel { + + if nextLevelDown { + leftChild = currentZ << 2 + rightChild = (currentZ << 2) + 1 + } else { + leftChild = (currentZ << 2) + 2 + rightChild = (currentZ << 2) + 3 + } + _, presentleft := ix.intersect[currentLevel+1][leftChild] + _, presentright := ix.intersect[currentLevel+1][rightChild] + + if presentleft { + todoListAtNextLevel = append(todoListAtNextLevel, leftChild) + } + + rightX, _ := morton.FromZ(rightChild) + + if presentright && (rightX <= xOriginAtNextLevel) { + todoListAtNextLevel = append(todoListAtNextLevel, rightChild) + } + } + todoListAtCurrentLevel = todoListAtNextLevel + } + return todoListAtCurrentLevel +} + +func (ix *PointIndex) determineTileOnOutside(zOrigin morton.Z, targetLevel Level, polygon geom.Polygon) bool { + xOrigin, yOrigin := morton.FromZ(zOrigin) + + // Binary search to find all tiles left of `zOrigin` that are intersected + intersectingTilesLeft := ix.findIntersectingTilesLeft(xOrigin, yOrigin, targetLevel) + + // Check all found tiles for intersecting segments + seen := make(map[SegmentIdx]bool) + intersectionAmount := 0 + + for _, currentZ := range intersectingTilesLeft { + for _, segmentIdx := range ix.registeredSegments[currentZ] { + if seen[segmentIdx] { + continue + } + seen[segmentIdx] = true + ring := polygon.LinearRings()[segmentIdx.ringIdx] + floatPt1 := ring[segmentIdx.pointIdx] + floatPt2 := ring[(segmentIdx.pointIdx+1)%len(ring)] + + y1 := intgeom.FromGeomOrd(floatPt1[1]) + y2 := intgeom.FromGeomOrd(floatPt2[1]) + + minY := min(y1, y2) + maxY := max(y1, y2) + + // Only count a segment if the northernmost vertex is unique and intersects the ray + switch { + case minY == maxY: + case maxY < intgeom.M(yOrigin): //nolint:gosec // Guaranteed by FromZ + case minY >= intgeom.M(yOrigin): //nolint:gosec // Guaranteed by FromZ + continue + default: + intersectionAmount++ + } + } + } + + return intersectionAmount%2 == 0 +} + +// This function kickstarts the recursion +func (ix *PointIndex) classifyNonIntersectingTiles(polygon geom.Polygon, targetLevel Level) { + ix.classifyNonIntersectingTilesHelper(polygon, targetLevel, 0, 0) +} + +// Recursively walk through quadtree, raycasing on unknown tiles +func (ix *PointIndex) classifyNonIntersectingTilesHelper(polygon geom.Polygon, targetLevel Level, currentLevel Level, currentZ morton.Z) { + tileData, present := ix.intersect[currentLevel][currentZ] + + if present && tileData.isIntersect { + if currentLevel == targetLevel { + return + } + nextLevel := currentLevel + 1 + for _, childZ := range getQuadrantZs(currentZ) { + ix.classifyNonIntersectingTilesHelper(polygon, targetLevel, nextLevel, childZ) + } + return + } + + levelDiff := targetLevel - currentLevel + targetZ := currentZ << (2 * levelDiff) + + isOutside := ix.determineTileOnOutside(targetZ, targetLevel, polygon) + + var newTileData TileData + + if isOutside { + newTileData = TileData{isOutside: true} + } else { + newTileData = TileData{isInside: true} + } + + ix.intersect[currentLevel][currentZ] = newTileData +} + +func (ix *PointIndex) classifyTiles(polygon geom.Polygon, targetLevel Level) { + ix.intersect = make(map[Level]map[morton.Z]TileData) + for i := Level(0); i <= targetLevel; i++ { + ix.intersect[i] = make(map[morton.Z]TileData) + } + ix.registerPolygonEdges(polygon, targetLevel) + ix.classifyNonIntersectingTiles(polygon, targetLevel) } //nolint:nestif diff --git a/pointindex/pointindex_test.go b/pointindex/pointindex_test.go index 40927ea..606dbe8 100644 --- a/pointindex/pointindex_test.go +++ b/pointindex/pointindex_test.go @@ -521,6 +521,493 @@ func TestPointIndex_lineIntersects(t *testing.T) { } } +func TestPointIndex_registerQuadrant(t *testing.T) { + btd := TileData{isIntersect: true} + tests := []struct { + name string + ix *PointIndex + setupIx func(ix *PointIndex) + x intgeom.M + y intgeom.M + level Level + + resultIntersect map[Level]map[morton.Z]TileData + resultSegments map[morton.Z][]SegmentIdx + }{ + { + name: "Registreer op hoogste level", + ix: newSimplePointIndex(2, 1), + x: 0, + y: 0, + level: 0, + resultIntersect: map[Level]map[morton.Z]TileData{0: {0: btd}}, + resultSegments: map[morton.Z][]SegmentIdx{0: {{1, 1}}}, + }, + { + name: "Registreer op dieper level", + ix: newSimplePointIndex(4, 1), + x: 1, + y: 1, + level: 2, + resultIntersect: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd}, + 2: {3: btd}, + }, + resultSegments: map[morton.Z][]SegmentIdx{3: {{1, 1}}}, + }, + { + name: "Registreer met aanwezige data", + ix: newSimplePointIndex(4, 1), + setupIx: func(ix *PointIndex) { + ix.registerQuadrant(2, 3, 2, SegmentIdx{1, 2}) + }, + x: 1, + y: 1, + level: 2, + resultIntersect: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd, 3: btd}, + 2: {3: btd, 14: btd}, + }, + resultSegments: map[morton.Z][]SegmentIdx{3: {{1, 1}}, 14: {{1, 2}}}, + }, + { + name: "Registreer met aanwezige data op hetzelfde punt", + ix: newSimplePointIndex(4, 1), + setupIx: func(ix *PointIndex) { + ix.registerQuadrant(1, 1, 2, SegmentIdx{1, 2}) + }, + x: 1, + y: 1, + level: 2, + resultIntersect: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd}, + 2: {3: btd}, + }, + resultSegments: map[morton.Z][]SegmentIdx{3: {{1, 2}, {1, 1}}}, + }, + { + name: "Registreer met aanwezige data op hetzelfde punt", + ix: newSimplePointIndex(4, 1), + setupIx: func(ix *PointIndex) { + ix.registerQuadrant(1, 1, 2, SegmentIdx{1, 1}) + }, + x: 1, + y: 1, + level: 2, + resultIntersect: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd}, + 2: {3: btd}, + }, + resultSegments: map[morton.Z][]SegmentIdx{3: {{1, 1}}}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + if tt.setupIx != nil { + tt.setupIx(tt.ix) + } + tt.ix.registerQuadrant(tt.x, tt.y, tt.level, SegmentIdx{1, 1}) + assert.Equal(t, tt.resultIntersect, tt.ix.intersect) + assert.Equal(t, tt.resultSegments, tt.ix.registeredSegments) + }) + } +} + +func TestPointIndex_registerQuadrantAndNeighbours(t *testing.T) { + btd := TileData{isIntersect: true} + tests := []struct { + name string + ix *PointIndex + x intgeom.M + y intgeom.M + level Level + result map[Level]map[morton.Z]TileData + }{ + { + name: "Insert at highest level", + ix: newSimplePointIndex(2, 1), + x: 0, + y: 0, + level: 0, + result: map[Level]map[morton.Z]TileData{0: {0: btd}}, + }, + { + name: "Insert maximal amount of neighbours.", + ix: newSimplePointIndex(3, 1), + x: 1, + y: 1, + level: 2, + result: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd, 1: btd, 2: btd, 3: btd}, + 2: { + 0: btd, 1: btd, 2: btd, + 3: btd, 4: btd, 6: btd, + 8: btd, 9: btd, 12: btd, + }, + }, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + tt.ix.registerQuadrantAndNeighbours(tt.x, tt.y, tt.level, SegmentIdx{1, 1}) + assert.Equal(t, tt.result, tt.ix.intersect) + }) + } +} + +func TestPointIndex_amanatides_woo(t *testing.T) { + btd := TileData{isIntersect: true} + tests := []struct { + name string + ix *PointIndex + line geom.Line + level Level + result map[Level]map[morton.Z]TileData + }{ + { + name: "Inserting a line", + ix: newSimplePointIndex(3, intgeom.ToGeomOrd(1)), + line: intgeom.Line{{1, 1}, {3, 5}}.ToGeomLine(), + level: 2, + result: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd, 1: btd, 2: btd, 3: btd}, + 2: { + 0: btd, 1: btd, 2: btd, 3: btd, + 4: btd, 6: btd, + 8: btd, 9: btd, 10: btd, 11: btd, + 12: btd, 14: btd, + }, + }, + }, + { + name: "Inserting a line near edge of index", + ix: newSimplePointIndex(3, intgeom.ToGeomOrd(1)), + line: intgeom.Line{{0, 0}, {2, 4}}.ToGeomLine(), + level: 2, + result: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd, 1: btd, 2: btd, 3: btd}, + 2: { + 0: btd, 1: btd, 2: btd, 3: btd, + 6: btd, + 8: btd, 9: btd, 10: btd, 11: btd, + 12: btd, 14: btd, + }, + }, + }, + { + name: "Downward line along right edge of index", + ix: newSimplePointIndex(2, intgeom.ToGeomOrd(1)), + line: intgeom.Line{{4, 4}, {4, 0}}.ToGeomLine(), + level: 2, + result: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {1: btd, 3: btd}, + 2: {5: btd, 7: btd, 13: btd, 15: btd}, + }, + }, + { + name: "Line within one tile", + ix: newSimplePointIndex(3, intgeom.ToGeomOrd(1)), + line: intgeom.Line{{1, 7}, {0, 8}}.ToGeomLine(), + level: 2, + result: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {2: btd}, + 2: {8: btd, 9: btd, 10: btd, 11: btd}, + }, + }, + { + name: "Intersecting diagonally in a down-right direction.", + ix: newSimplePointIndex(3, intgeom.ToGeomOrd(1)), + line: intgeom.Line{{3, 5}, {5, 3}}.ToGeomLine(), + level: 2, + result: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd, 1: btd, 2: btd, 3: btd}, + 2: { + 1: btd, 2: btd, 3: btd, + 4: btd, 5: btd, 6: btd, 7: btd, + 8: btd, 9: btd, 10: btd, 11: btd, + 12: btd, 13: btd, 14: btd, 15: btd, + }, + }, + }, + { + name: "Intersecting diagonally in a up-left direction.", + ix: newSimplePointIndex(3, intgeom.ToGeomOrd(1)), + line: intgeom.Line{{5, 3}, {3, 5}}.ToGeomLine(), + level: 2, + result: map[Level]map[morton.Z]TileData{ + 0: {0: btd}, + 1: {0: btd, 1: btd, 2: btd, 3: btd}, + 2: { + 1: btd, 2: btd, 3: btd, + 4: btd, 5: btd, 6: btd, 7: btd, + 8: btd, 9: btd, 10: btd, 11: btd, + 12: btd, 13: btd, 14: btd, 15: btd, + }, + }, + }, + } + + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + tt.ix.AmanatidesWoo(tt.line, tt.level, 1, 1) + assert.Equal(t, tt.result, tt.ix.intersect) + }) + } +} + +func TestPointIndex_findIntersectingTilesLeft(t *testing.T) { + tests := []struct { + name string + ix *PointIndex + xOrigin uint + yOrigin uint + targetLevel Level + tilesHit []morton.Z + }{ + { + name: "Empty tiles hit", + ix: newSimplePointIndex(3, 1), + xOrigin: 2, + yOrigin: 1, + targetLevel: 2, + tilesHit: []morton.Z{}, + }, + { + name: "Single tile hit", + ix: newSimplePointIndex(3, 1), + xOrigin: 2, + yOrigin: 1, + targetLevel: 2, + tilesHit: []morton.Z{3}, + }, + { + name: "Single tile hit, but missed by algorithm", + ix: newSimplePointIndex(3, 1), + xOrigin: 2, + yOrigin: 1, + targetLevel: 2, + tilesHit: []morton.Z{1, 3}, + }, + { + name: "Single tile hit, but missed by algorithm", + ix: newSimplePointIndex(3, 1), + xOrigin: 2, + yOrigin: 1, + targetLevel: 2, + tilesHit: []morton.Z{1, 3, 7}, + }, + } + for _, tt := range tests { + expected := []morton.Z{} + t.Run(tt.name, func(t *testing.T) { + expected = []morton.Z{} + for _, currentZ := range tt.tilesHit { + currentX, currentY := morton.FromZ(currentZ) + tt.ix.registerQuadrant(intgeom.M(currentX), intgeom.M(currentY), tt.targetLevel, SegmentIdx{0, 0}) //nolint:gosec // Guaranteed by FromZ + if currentY == tt.yOrigin && currentX <= tt.xOrigin { + expected = append(expected, currentZ) + } + } + actual := tt.ix.findIntersectingTilesLeft(tt.xOrigin, tt.yOrigin, tt.targetLevel) + assert.Equal(t, expected, actual) + }) + } +} + +func TestPointIndex_determineTileOnOutside(t *testing.T) { + tests := []struct { + name string + ix *PointIndex + polygon geom.Polygon + targetLevel Level + targetCell morton.Z + expected bool + }{ + { + name: "Empty polygon", + ix: newSimplePointIndex(3, 1), + targetLevel: 2, + polygon: geom.Polygon{}, + targetCell: 10, + expected: true, + }, + { + name: "Inside square", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + targetLevel: 3, + polygon: geom.Polygon{{ + {0, 0}, + {0, intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), 0}, + }}, + targetCell: 15, + expected: false, + }, + { + name: "Edge lies on ray", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + targetLevel: 3, + polygon: geom.Polygon{{ + {0, 0}, + {0, intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(7), intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(7), 0}, + }}, + targetCell: 16, + expected: true, + }, + { + name: "Ray passes through corner, no hit", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + targetLevel: 3, + polygon: geom.Polygon{{ + {0, 0}, + {0, intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), 0}, + {intgeom.ToGeomOrd(1), intgeom.ToGeomOrd(6)}, + }}, + targetCell: 15, + expected: false, + }, + { + name: "Ray passes through corner, hit", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + targetLevel: 3, + polygon: geom.Polygon{{ + {0, 0}, + {intgeom.ToGeomOrd(1), intgeom.ToGeomOrd(6)}, + {0, intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), 0}, + }}, + targetCell: 15, + expected: false, + }, + { + name: "Within a hole", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + targetLevel: 3, + polygon: geom.Polygon{ + { + {0, 0}, + {0, intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), 0}, + }, + { + {intgeom.ToGeomOrd(1), intgeom.ToGeomOrd(1)}, + {intgeom.ToGeomOrd(1), intgeom.ToGeomOrd(14)}, + {intgeom.ToGeomOrd(14), intgeom.ToGeomOrd(14)}, + {intgeom.ToGeomOrd(14), intgeom.ToGeomOrd(1)}, + }, + }, + targetCell: 15, + expected: true, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + tt.ix.registerPolygonEdges(tt.polygon, tt.targetLevel) + actual := tt.ix.determineTileOnOutside(tt.targetCell, tt.targetLevel, tt.polygon) + assert.Equal(t, tt.expected, actual) + }) + } +} + +func TestPointIndex_classifyTiles(t *testing.T) { + iI := TileData{isIntersect: true} + io := TileData{isOutside: true} + ii := TileData{isInside: true} + tests := []struct { + name string + ix *PointIndex + polygon geom.Polygon + targetLevel Level + expected map[Level]map[morton.Z]TileData + }{ + { + name: "Empty polygon", + ix: newSimplePointIndex(2, 1), + polygon: geom.Polygon{}, + targetLevel: 1, + expected: map[Level]map[morton.Z]TileData{ + 0: {0: TileData{isOutside: true}}, + 1: {}, + }, + }, + { + name: "Triangle with interior tile", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + polygon: geom.Polygon{{ + {intgeom.ToGeomOrd(0), intgeom.ToGeomOrd(0)}, + {intgeom.ToGeomOrd(5), intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), intgeom.ToGeomOrd(5)}, + }}, + targetLevel: 3, + expected: map[Level]map[morton.Z]TileData{ + 0: {0: iI}, + 1: {0: iI, 1: iI, 2: iI, 3: iI}, + 2: {0: iI, 1: iI, 2: iI, 3: iI, 4: iI, 5: iI, 6: iI, 7: iI, 8: iI, 9: iI, 10: iI, 11: iI, 12: iI, 13: iI, 14: iI, 15: iI}, + 3: {0: iI, 1: iI, 2: iI, 3: iI, 4: iI, 5: iI, 6: iI, 7: iI, 8: iI, 9: iI, 10: iI, 11: iI, 12: iI, 13: iI, 14: iI, 15: ii, 16: iI, 17: iI, 18: iI, 19: iI, 20: iI, 21: io, 22: iI, 23: iI, 24: iI, 25: iI, 26: iI, 27: iI, 28: iI, 29: iI, 30: iI, 31: iI, 32: iI, 33: iI, 34: iI, 35: iI, 36: iI, 37: iI, 38: iI, 39: iI, 40: iI, 41: iI, 42: io, 43: iI, 44: iI, 45: iI, 46: iI, 47: iI, 48: iI, 49: iI, 50: iI, 51: iI, 52: iI, 53: iI, 54: iI, 55: iI, 56: iI, 57: iI, 58: iI, 59: iI, 60: iI, 61: io, 62: io, 63: io}, + }, + }, + { + name: "Triangle with pinched ring", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + polygon: geom.Polygon{{ + {intgeom.ToGeomOrd(0), intgeom.ToGeomOrd(0)}, + {intgeom.ToGeomOrd(5), intgeom.ToGeomOrd(15)}, + {intgeom.ToGeomOrd(15), intgeom.ToGeomOrd(5)}, + }, { + {intgeom.ToGeomOrd(0), intgeom.ToGeomOrd(0)}, + {intgeom.ToGeomOrd(5), intgeom.ToGeomOrd(14)}, + {intgeom.ToGeomOrd(14), intgeom.ToGeomOrd(5)}, + }}, + targetLevel: 3, + expected: map[Level]map[morton.Z]TileData{ + 0: {0: iI}, + 1: {0: iI, 1: iI, 2: iI, 3: iI}, + 2: {0: iI, 1: iI, 2: iI, 3: iI, 4: iI, 5: iI, 6: iI, 7: iI, 8: iI, 9: iI, 10: iI, 11: iI, 12: iI, 13: iI, 14: iI, 15: iI}, + 3: {0: iI, 1: iI, 2: iI, 3: iI, 4: iI, 5: iI, 6: iI, 7: iI, 8: iI, 9: iI, 10: iI, 11: iI, 12: iI, 13: iI, 14: iI, 15: iI, 16: iI, 17: iI, 18: iI, 19: iI, 20: iI, 21: io, 22: iI, 23: iI, 24: iI, 25: iI, 26: iI, 27: iI, 28: iI, 29: iI, 30: iI, 31: iI, 32: iI, 33: iI, 34: iI, 35: iI, 36: iI, 37: iI, 38: iI, 39: iI, 40: iI, 41: iI, 42: io, 43: iI, 44: iI, 45: iI, 46: iI, 47: iI, 48: iI, 49: iI, 50: iI, 51: iI, 52: iI, 53: iI, 54: iI, 55: iI, 56: iI, 57: iI, 58: iI, 59: iI, 60: iI, 61: io, 62: io, 63: io}, + }, + }, + { + name: "Outside detection at higher level", + ix: newSimplePointIndex(4, intgeom.ToGeomOrd(1)), + polygon: geom.Polygon{{ + {intgeom.ToGeomOrd(1), intgeom.ToGeomOrd(1)}, + {intgeom.ToGeomOrd(1), intgeom.ToGeomOrd(12)}, + {intgeom.ToGeomOrd(2), intgeom.ToGeomOrd(12)}, + {intgeom.ToGeomOrd(2), intgeom.ToGeomOrd(1)}, + }}, + targetLevel: 2, + expected: map[Level]map[morton.Z]TileData{ + 0: {0: iI}, + 1: {0: iI, 1: io, 2: iI, 3: io}, + 2: {0: iI, 1: iI, 2: iI, 3: iI, 8: iI, 9: iI, 10: iI, 11: iI}, + }, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + tt.ix.classifyTiles(tt.polygon, tt.targetLevel) + assert.Equal(t, tt.expected, tt.ix.intersect) + }) + } +} + func newSimplePointIndex(deepestLevel Level, cellSize float64) *PointIndex { deepestSize := mathhelp.Pow2(deepestLevel) span := cellSize * float64(deepestSize)