diff --git a/intgeom/intgeom.go b/intgeom/intgeom.go index 8285e28..00ac35e 100644 --- a/intgeom/intgeom.go +++ b/intgeom/intgeom.go @@ -56,7 +56,7 @@ func FromGeomOrd(o float64) M { // TODO implement this with integers func SegmentIntersect(l1, l2 Line) ([2]int64, bool) { intersection, intersects := planar.SegmentIntersect(l1.ToGeomLine(), l2.ToGeomLine()) - intIntersection := [2]int64{FromGeomOrd(intersection[0]), FromGeomOrd(intersection[0])} + intIntersection := [2]int64{FromGeomOrd(intersection[0]), FromGeomOrd(intersection[1])} return intIntersection, intersects } diff --git a/main.go b/main.go index fd9f54d..d619b2f 100644 --- a/main.go +++ b/main.go @@ -36,6 +36,8 @@ const ( IGNOREOUTSIDEGRID string = `ignoreoutsidegrid` REVERSEWINDINGORDER string = `reversewindingorder` ENCODETILES string = `encodetiles` + TILEBUFFER string = `tilebuffer` + USELINETRACE string = `uselinetrace` MVTSOURCE string = `mvtSourceGpkg` MVTOUTDIR string = `mvtOutDir` @@ -133,6 +135,22 @@ func main() { Required: false, EnvVars: []string{strcase.ToScreamingSnake(ENCODETILES)}, }, + &cli.UintFlag{ + Name: TILEBUFFER, + Aliases: []string{"buf"}, + Usage: "Buffer (in internal pixels) used to select the tiles a polygon's geometry touches.", + Value: 0, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(TILEBUFFER)}, + }, + &cli.BoolFlag{ + Name: USELINETRACE, + Aliases: []string{"lt"}, + Usage: "Use precise line-trace tile classification instead of a simple buffered bounding box to select tiles.", + Value: false, + Required: false, + EnvVars: []string{strcase.ToScreamingSnake(USELINETRACE)}, + }, }, Action: func(c *cli.Context) error { tileMatrixSet, err := tms20.LoadEmbeddedTileMatrixSet(c.String(TILEMATRIXSET)) @@ -167,6 +185,8 @@ func main() { IgnoreOutsideGrid: c.Bool(IGNOREOUTSIDEGRID), ReverseWindingOrder: c.Bool(REVERSEWINDINGORDER), EncodeTiles: c.Bool(ENCODETILES), + Buffer: c.Uint(TILEBUFFER), + UseLineTrace: c.Bool(USELINETRACE), } for _, tmID := range tileMatrixIDs { gpkgTargets[tmID] = initGPKGTarget(targetPathFmt, tmID, overwrite, pagesize, c.Bool(ENCODETILES)) @@ -273,7 +293,7 @@ func injectSuffixIntoPath(p string) string { func processBySnapping(source processing.Source, targets map[tms20.TMID]processing.Target, tileMatrixSet tms20.TileMatrixSet, snapConfig snap.Config) { processing.ProcessFeatures(source, targets, func(p geom.Polygon, tmIDs []tms20.TMID) map[tms20.TMID]processing.SnapResult { return snap.SnapPolygon(p, tileMatrixSet, tmIDs, snapConfig) - }, snapConfig.EncodeTiles) + }, snapConfig.EncodeTiles, snapConfig.Buffer) } // Initialize resources for creating vecotrtiles and delegate to processing. diff --git a/pointindex/detect.go b/pointindex/detect.go new file mode 100644 index 0000000..026770f --- /dev/null +++ b/pointindex/detect.go @@ -0,0 +1,344 @@ +package pointindex + +import ( + "math" + + "github.com/go-spatial/geom" + "github.com/pdok/texel/intgeom" + "github.com/pdok/texel/mathhelp" + "github.com/pdok/texel/morton" + "github.com/pdok/texel/tile" + "github.com/pdok/texel/tms20" +) + +type SegmentIdx struct { + ringIdx int + pointIdx int +} + +// RegisterFunc marks a tile at (xCoord, yCoord) (in tile coordinates at +// level l) as touched by the segment identified by segmentIdx. +type RegisterFunc func(xCoord, yCoord uint, l Level, segmentIdx SegmentIdx) + +// Logic for searching tile grid for tiles that intersect a line. +// Detect line direction and trace that direction +func (ix *PointIndex) lineTrace(line geom.Line, tileLevel, intPixLevel Level, ringIdx int, pointIdx int, buffer uint, register RegisterFunc) { + intLine := intgeom.FromGeomLine(line) + idx := SegmentIdx{ + ringIdx: ringIdx, + pointIdx: pointIdx, + } + + var dx, dy int + if intLine.Point1().X() < intLine.Point2().X() { + dx = 1 + } else { + dx = -1 + } + if intLine.Point1().Y() < intLine.Point2().Y() { + dy = 1 + } else { + dy = -1 + } + + startTileX, startTileY := ix.findTile(intLine.Point1(), tileLevel) + startX, startY := int(startTileX), int(startTileY) //nolint:gosec // G115 + bufferSize := ix.getResolution(intPixLevel) * intgeom.M(buffer) //nolint:gosec // G115 + + // Register tiles at start + ix.tryRegisterTile(intLine, startX-dx, startY+dy, tileLevel, bufferSize, register, idx) + ix.tryRegisterTile(intLine, startX-dx, startY, tileLevel, bufferSize, register, idx) + ix.tryRegisterTile(intLine, startX-dx, startY-dy, tileLevel, bufferSize, register, idx) + ix.tryRegisterTile(intLine, startX, startY-dy, tileLevel, bufferSize, register, idx) + ix.tryRegisterTile(intLine, startX+dx, startY-dy, tileLevel, bufferSize, register, idx) + + type coord struct{ x, y int } + + // Advance over an anti-diagonal, keeping track of the last two "hits" + frontier := []coord{{startX, startY}} + prevFrontier := make([]coord, 0) + ix.tryRegisterTile(intLine, startX, startY, tileLevel, bufferSize, register, idx) + + for len(prevFrontier)+len(frontier) > 0 { + nextSet := make(map[coord]bool, len(frontier)*2+len(prevFrontier)) + for _, cur := range frontier { + nextSet[coord{cur.x + dx, cur.y}] = true + nextSet[coord{cur.x, cur.y + dy}] = true + } + for _, prev := range prevFrontier { + nextSet[coord{prev.x + dx, prev.y + dy}] = true + } + + next := make([]coord, 0, len(nextSet)) + for c := range nextSet { + if ix.tryRegisterTile(intLine, c.x, c.y, tileLevel, bufferSize, register, idx) { + next = append(next, c) + } + } + prevFrontier = frontier + frontier = next + } +} + +func (ix *PointIndex) getResolution(level Level) intgeom.M { + return ix.intExtent.XSpan() / int64(mathhelp.Pow2(level)) //nolint:gosec // G115 +} + +func (ix *PointIndex) getInternalPixelLevel(deepestTIMID tms20.TMID) Level { + levelDiff := uint(math.Log2(float64(ix.tilePixels))) + uint(math.Log2(float64(ix.internalPixels))) + return uint(deepestTIMID) + levelDiff //nolint:gosec // G115 +} + +// Logic for deciding whether line intersects tile with buffer. +func (ix *PointIndex) tryRegisterTile(line intgeom.Line, x, y int, l Level, bufferSize intgeom.M, register RegisterFunc, idx SegmentIdx) bool { + maxCoord := int(mathhelp.Pow2(l)) - 1 //nolint:gosec // G115 + if x < 0 || y < 0 || x > maxCoord || y > maxCoord { + return false // out of bounds: no tile there, nothing to register + } + ux, uy := uint(x), uint(y) + extent, _ := ix.getQuadrantExtentAndCentroid(l, ux, uy, ix.intExtent) + if tileIntersectsLine(line, extent, bufferSize) { + register(ux, uy, l, idx) + return true + } + return false +} + +func (ix *PointIndex) findTile(p *intgeom.Point, l Level) (tileX, tileY uint) { + levelDiff := ix.deepestLevel - l + //nolint:gosec // G115 + deepestTileX := uint(p.X()-ix.intExtent.MinX()) / uint(ix.deepestRes) + //nolint:gosec // G115 + deepestTileY := uint(p.Y()-ix.intExtent.MinY()) / uint(ix.deepestRes) + tileX = deepestTileX >> levelDiff + tileY = deepestTileY >> levelDiff + return +} + +func tileIntersectsLine(line intgeom.Line, extent intgeom.Extent, buffer intgeom.M) bool { + bufferedExtent := intgeom.Extent{ + extent.MinX() - buffer, + extent.MinY() - buffer, + extent.MaxX() + buffer, + extent.MaxY() + buffer, + } + + return lineIntersects(line, bufferedExtent) +} + +type TileClassification int + +const ( + ClassificationUnknown TileClassification = iota + ClassificationIntersect + ClassificationInside + ClassificationOutside +) + +// Create register logic, then loop over polygon edges +func (ix *PointIndex) registerPolygonEdges(polygon geom.Polygon, tmsID tms20.TMID, buffer uint) (segments map[morton.Z][]SegmentIdx, classification map[Level]map[morton.Z]TileClassification) { + segments = make(map[morton.Z][]SegmentIdx) + tileLevel := Level(tmsID) //nolint:gosec // G115 + intPixLevel := ix.getInternalPixelLevel(tmsID) + classification = make(map[Level]map[morton.Z]TileClassification, tileLevel+1) + for l := range tileLevel + 1 { + classification[l] = make(map[morton.Z]TileClassification) + } + + // Logic for keeping track of intersected tiles + var markIntersected func(l Level, z morton.Z) + markIntersected = func(l Level, z morton.Z) { + if classification[l][z] == ClassificationIntersect { + return + } + classification[l][z] = ClassificationIntersect + if l == 0 { + return + } + markIntersected(l-1, z>>2) + } + + // Logic for keeping track which segment hits which tile + register := func(x, y uint, l Level, segmentIdx SegmentIdx) { + z := morton.MustToZ(x, y) + segments[z] = append(segments[z], segmentIdx) + markIntersected(l, z) + } + + // Loop + for ringIdx, ring := range polygon.LinearRings() { + for pointIdx := range ring { + line := geom.Line{ring[pointIdx], ring[(pointIdx+1)%len(ring)]} + ix.lineTrace(line, tileLevel, intPixLevel, ringIdx, pointIdx, buffer, register) + } + } + return segments, classification +} + +// Given a tile, binary search for tiles left of it that have segments on them +// Note: tileLevel must be the actual tile level and x and y coordinates for that level +func (ix *PointIndex) findIntersectingTilesLeft(x, y, tileLevel Level, classification map[Level]map[morton.Z]TileClassification) []morton.Z { + intersectingCurrentLevel := []morton.Z{0} + var intersectingNextLevel []morton.Z + var leftChild, rightChild morton.Z + for currentLevel := range tileLevel { + intersectingNextLevel = make([]morton.Z, 0) + + xAtNextLevel := x >> (tileLevel - currentLevel - 1) + yAtNextLevel := y >> (tileLevel - currentLevel - 1) + + nextLevelDown := yAtNextLevel%2 == 0 + for _, z := range intersectingCurrentLevel { + if nextLevelDown { + leftChild = z << 2 + rightChild = (z << 2) + 1 + } else { + leftChild = (z << 2) + 2 + rightChild = (z << 2) + 3 + } + + if _, present := classification[currentLevel+1][leftChild]; present { + intersectingNextLevel = append(intersectingNextLevel, leftChild) + } + rightX, _ := morton.FromZ(rightChild) + if _, present := classification[currentLevel+1][rightChild]; present && rightX <= xAtNextLevel { + intersectingNextLevel = append(intersectingNextLevel, rightChild) + } + } + intersectingCurrentLevel = intersectingNextLevel + } + return intersectingCurrentLevel +} + +// Implements raycast to detect whether a non-intersecting tile lies on the inside or outside. +// tileLevel needs to be the tile level and z needs to be a morton coordinate for this level +func (ix *PointIndex) classifyNonIntersectingTile(z morton.Z, tileLevel Level, segments map[morton.Z][]SegmentIdx, classification map[Level]map[morton.Z]TileClassification, polygon geom.Polygon) TileClassification { + x, y := morton.FromZ(z) + intersectingTilesLeft := ix.findIntersectingTilesLeft(x, y, tileLevel, classification) + + tileHeightCoord := ix.getResolution(tileLevel)*intgeom.M(y) + ix.intExtent.MinY() //nolint:gosec // G115 + + seen := make(map[SegmentIdx]bool) + numIntersections := 0 + + // Raycast and loop over intersecting segments + for _, z := range intersectingTilesLeft { + for _, segment := range segments[z] { + if seen[segment] { + continue + } + seen[segment] = true + ring := polygon.LinearRings()[segment.ringIdx] + + y1 := intgeom.FromGeomOrd(ring[segment.pointIdx][1]) + y2 := intgeom.FromGeomOrd(ring[(segment.pointIdx+1)%len(ring)][1]) + + minY := min(y1, y2) + maxY := max(y1, y2) + + // Logic for deciding relevant segments + // Needed for when ray intersects vertex of polygon + switch { + case minY == maxY: + case maxY < tileHeightCoord: + case minY >= tileHeightCoord: + default: + numIntersections++ + } + } + } + + // Result determined by "Jordan Curve Theorem" + if numIntersections%2 == 0 { + return ClassificationOutside + } + return ClassificationInside +} + +func getChildren(z morton.Z) [4]morton.Z { + shift := z << 2 + return [4]morton.Z{shift, shift + 1, shift + 2, shift + 3} +} + +// Classify all tiles on being inside or outside using quadtree. An "outside" or "inside" +// classification of a parent means children have the same class. +func (ix *PointIndex) classifyNonIntersectingTiles(targetLevel, currentLevel Level, currentZ morton.Z, containsAll bool, segments map[morton.Z][]SegmentIdx, classification map[Level]map[morton.Z]TileClassification, polygon geom.Polygon) { + if targetLevel == currentLevel { + return + } + + children := getChildren(currentZ) + intersectingChildren := make([]morton.Z, 0, 4) + nextLevel := currentLevel + 1 + + // Process unknown children + for _, child := range children { + if _, present := classification[nextLevel][child]; present { + intersectingChildren = append(intersectingChildren, child) + continue + } + if containsAll { + classification[nextLevel][child] = ClassificationOutside + } else { + childAtTargetLevel := child << ((targetLevel - nextLevel) * 2) + classification[nextLevel][child] = ix.classifyNonIntersectingTile(childAtTargetLevel, targetLevel, segments, classification, polygon) + } + } + + containsAll = containsAll && len(intersectingChildren) < 2 + + // Recurse for intersecting children + for _, child := range intersectingChildren { + if _, present := classification[nextLevel][child]; present { + ix.classifyNonIntersectingTiles(targetLevel, nextLevel, child, containsAll, segments, classification, polygon) + } + } +} + +func (ix *PointIndex) ClassifyTiles(polygon geom.Polygon, tmsID tms20.TMID, buffer uint) map[Level]map[morton.Z]TileClassification { + targetLevel := Level(tmsID) //nolint:gosec // G115 + segments, classification := ix.registerPolygonEdges(polygon, tmsID, buffer) + ix.classifyNonIntersectingTiles(targetLevel, 0, 0, true, segments, classification, polygon) + return classification +} + +// Plumbing function +// Process output of ClassifyTiles +func (ix *PointIndex) GetLineTraceResult(polygon geom.Polygon, tmsID tms20.TMID, buffer uint) (tiles []tile.Tile) { + classification := ix.ClassifyTiles(polygon, tmsID, buffer) + tileLevel := Level(tmsID) //nolint:gosec // G115 integers < 40 + tiles = make([]tile.Tile, 0) + + for level, classificationAtLevel := range classification { + for z, class := range classificationAtLevel { + switch class { + case ClassificationInside: + size := mathhelp.Pow2((tileLevel - level) * 2) + baseZ := z << (2 * (tileLevel - level)) + for i := range size { + tiles = append(tiles, ix.makeTile(baseZ+i, tileLevel, true)) + } + case ClassificationOutside: + continue + case ClassificationIntersect: + if level == tileLevel { + tiles = append(tiles, ix.makeTile(z, tileLevel, false)) + } + case ClassificationUnknown: + panic("ClassificationUnknown tile for polygon during linetrace") + } + } + } + return tiles +} + +func (ix *PointIndex) makeTile(z morton.Z, l Level, isContained bool) tile.Tile { + x, y := morton.FromZ(z) + extent, _ := ix.getQuadrantExtentAndCentroid(l, x, y, ix.intExtent) + return tile.Tile{ + Extent: extent, + X: x, + Y: mathhelp.Pow2(l) - 1 - y, + IsContained: isContained, + } +} diff --git a/pointindex/detect_test.go b/pointindex/detect_test.go new file mode 100644 index 0000000..9ebed29 --- /dev/null +++ b/pointindex/detect_test.go @@ -0,0 +1,723 @@ +package pointindex + +import ( + "sort" + "testing" + + "github.com/go-spatial/geom" + "github.com/stretchr/testify/assert" + + "github.com/pdok/texel/intgeom" + "github.com/pdok/texel/mathhelp" + "github.com/pdok/texel/morton" + "github.com/pdok/texel/tms20" +) + +// registeredTile records a single call to a RegisterFunc during a test. +type registeredTile struct { + x, y uint + segmentIdx SegmentIdx +} + +// Register function for testing +func recordingRegister(dst *[]registeredTile) RegisterFunc { + return func(tileX, tileY uint, l Level, segmentIdx SegmentIdx) { + maxCoord := mathhelp.Pow2(l) - 1 + if tileX > maxCoord || tileY > maxCoord { + return + } + *dst = append(*dst, registeredTile{tileX, tileY, segmentIdx}) + } +} + +// order and deduplicate slice of registeredTile (for testing) +func uniqueTileCoords(records []registeredTile) [][2]uint { + seen := map[[2]uint]bool{} + var out [][2]uint + for _, r := range records { + key := [2]uint{r.x, r.y} + if !seen[key] { + seen[key] = true + out = append(out, key) + } + } + sort.Slice(out, func(i, j int) bool { + if out[i][0] != out[j][0] { + return out[i][0] < out[j][0] + } + return out[i][1] < out[j][1] + }) + return out +} + +// Pointindex with centre not at (0,0) +func newOffsetPointIndex(deepestLevel Level, cellSize, originX, originY float64) *PointIndex { + deepestSize := mathhelp.Pow2(deepestLevel) + span := cellSize * float64(deepestSize) + intExtent := intgeom.Extent{ + intgeom.FromGeomOrd(originX), intgeom.FromGeomOrd(originY), + intgeom.FromGeomOrd(originX + span), intgeom.FromGeomOrd(originY + span), + } + return &PointIndex{ + Quadrant: Quadrant{intExtent: intExtent}, + deepestLevel: deepestLevel, + deepestSize: deepestSize, + //nolint:gosec // G115 + deepestRes: intExtent.XSpan() / int64(deepestSize), + } +} + +func TestLineTrace_TilesTouched(t *testing.T) { + tests := []struct { + name string + deepestLevel Level + l Level + cellSize float64 + originX float64 + originY float64 + line geom.Line + buffer uint + want [][2]uint + }{ + { + name: "no buffer, +x+y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{14, 10}}, + buffer: 0, + want: [][2]uint{{0, 0}, {1, 0}, {1, 1}, {2, 1}, {2, 2}, {3, 2}}, + }, + { + name: "no buffer, -x-y diagonal (reverse of +x+y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 10}, geom.Point{2, 2}}, + buffer: 0, + want: [][2]uint{{0, 0}, {1, 0}, {1, 1}, {2, 1}, {2, 2}, {3, 2}}, + }, + { + name: "no buffer, +x-y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 10}, geom.Point{14, 2}}, + buffer: 0, + want: [][2]uint{{0, 2}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {3, 0}}, + }, + { + name: "no buffer, -x+y diagonal (reverse of +x-y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 2}, geom.Point{2, 10}}, + buffer: 0, + want: [][2]uint{{0, 2}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {3, 0}}, + }, + { + name: "buffer, +x+y diagonal (reaches beyond index limit)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{14, 9}}, + buffer: 2, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {2, 2}, {3, 1}, {3, 2}}, + }, + { + name: "buffer, -x-y diagonal (ververse of +x+y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 9}, geom.Point{2, 2}}, + buffer: 2, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}, {1, 2}, {2, 0}, {2, 1}, {2, 2}, {3, 1}, {3, 2}}, + }, + { + name: "buffer, +x-y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 10}, geom.Point{14, 2}}, + buffer: 2, + want: [][2]uint{{0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, {1, 3}, {2, 0}, {2, 1}, {2, 2}, {3, 0}, {3, 1}}, + }, + { + name: "buffer, -x+y diagonal: (reverse of +x-y)", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{14, 2}, geom.Point{2, 10}}, + buffer: 2, + want: [][2]uint{{0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, {1, 3}, {2, 0}, {2, 1}, {2, 2}, {3, 0}, {3, 1}}, + }, + { + name: "horizontal line", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 5}, geom.Point{14, 5}}, + buffer: 0, + want: [][2]uint{{0, 1}, {1, 1}, {2, 1}, {3, 1}}, + }, + { + name: "vertical line", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{5, 2}, geom.Point{5, 14}}, + buffer: 0, + want: [][2]uint{{1, 0}, {1, 1}, {1, 2}, {1, 3}}, + }, + { + name: "end in upper right corner", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{3, 3}}, + buffer: 0, + want: [][2]uint{{0, 0}}, + }, + { + name: "alongside upper edge", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{2, 2}, geom.Point{3, 3}}, + buffer: 1, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + { + name: "start in corner, register tiles behind you", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{3, 3}, geom.Point{2, 2}}, + buffer: 1, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + { + name: "no buffer: walk across tile corner", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{6, 6}, geom.Point{10, 10}}, + buffer: 0, + want: [][2]uint{{1, 1}, {1, 2}, {2, 1}, {2, 2}}, + }, + { + name: "buffer: walk across tile corner", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{6, 6}, geom.Point{7, 7}}, + buffer: 1, + want: [][2]uint{{1, 1}, {1, 2}, {2, 1}, {2, 2}}, + }, + { + name: "buffer clips upper right corner of tile", + deepestLevel: 4, l: 2, cellSize: 1.0, + line: geom.Line{geom.Point{12, 14}, geom.Point{14, 12}}, + buffer: 1, + want: [][2]uint{{2, 2}, {2, 3}, {3, 2}, {3, 3}}, // Arguably {2,2} should not be here + }, + { + name: "non-zero origin: same relative result as +x+y diagonal", + deepestLevel: 4, l: 2, cellSize: 1.0, originX: 100, originY: 100, + line: geom.Line{geom.Point{102, 102}, geom.Point{114, 110}}, + buffer: 0, + want: [][2]uint{{0, 0}, {1, 0}, {1, 1}, {2, 1}, {2, 2}, {3, 2}}, + }, + { + name: "non-zero origin with buffer: same relative result as buffered corner case", + deepestLevel: 4, l: 2, cellSize: 1.0, originX: 100, originY: 100, + line: geom.Line{geom.Point{103, 103}, geom.Point{105, 105}}, + buffer: 1, + want: [][2]uint{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + } + + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + ix := newOffsetPointIndex(tt.deepestLevel, tt.cellSize, tt.originX, tt.originY) + + var recorded []registeredTile + ix.lineTrace(tt.line, tt.l, tt.deepestLevel, 0, 0, tt.buffer, recordingRegister(&recorded)) + + got := uniqueTileCoords(recorded) + assert.Equal(t, tt.want, got) + }) + } +} + +// helper for building classification data +func buildClassification(tuples ...[3]uint) map[Level]map[morton.Z]TileClassification { + classification := make(map[Level]map[morton.Z]TileClassification) + for _, tuple := range tuples { + level, x, y := tuple[0], tuple[1], tuple[2] + if classification[level] == nil { + classification[level] = make(map[morton.Z]TileClassification) + } + classification[level][morton.MustToZ(x, y)] = ClassificationIntersect + } + return classification +} + +func TestPointIndex_findIntersectingTilesLeft(t *testing.T) { + tests := []struct { + name string + x, y uint + targetLevel Level + classification map[Level]map[morton.Z]TileClassification + want []morton.Z + }{ + { + name: "trivial example at level 0", + x: 0, + y: 0, + targetLevel: 0, + classification: buildClassification(), + want: []morton.Z{0}, + }, + { + name: "no tiles", + x: 3, + y: 3, + targetLevel: 2, + classification: buildClassification(), + want: []morton.Z{}, + }, + { + name: "single tile registered, check left from this tile.", + x: 0, + y: 0, + targetLevel: 3, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 0}, + [3]uint{2, 0, 0}, + [3]uint{3, 0, 0}, + ), + want: []morton.Z{0}, + }, + { + name: "single tile registered, check this tile, requires rightChild", + x: 1, + y: 0, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 0}, + ), + want: []morton.Z{1}, + }, + { + name: "tile registered right of target", + x: 0, + y: 0, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 0}, + ), + want: []morton.Z{}, + }, + { + name: "use odd y coordinate", + x: 0, + y: 1, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 1}, + ), + want: []morton.Z{2}, + }, + { + name: "use odd y and rightChild", + x: 0, + y: 1, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 1}, + ), + want: []morton.Z{}, + }, + { + name: "ordd y and rightChild, want to find it", + x: 1, + y: 1, + targetLevel: 1, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 1, 1}, + ), + want: []morton.Z{3}, + }, + { + name: "register entire row", + x: 3, + y: 0, + targetLevel: 2, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 0}, + [3]uint{1, 1, 0}, + [3]uint{2, 0, 0}, + [3]uint{2, 1, 0}, + [3]uint{2, 2, 0}, + [3]uint{2, 3, 0}, + ), + want: []morton.Z{0, 1, 4, 5}, + }, + { + name: "corners of grid populated", + x: 7, + y: 7, + targetLevel: 3, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 0}, [3]uint{1, 0, 1}, [3]uint{1, 1, 0}, [3]uint{1, 1, 1}, + [3]uint{2, 0, 0}, [3]uint{2, 0, 3}, [3]uint{2, 3, 0}, [3]uint{2, 3, 3}, + [3]uint{3, 0, 0}, [3]uint{3, 7, 0}, [3]uint{3, 0, 7}, [3]uint{3, 7, 7}, + ), + want: []morton.Z{42, 63}, + }, + { + name: "general test with two disjoint tiles to pick up", + x: 6, + y: 7, + targetLevel: 3, + classification: buildClassification( + [3]uint{0, 0, 0}, + [3]uint{1, 0, 1}, [3]uint{1, 1, 1}, + [3]uint{2, 0, 3}, [3]uint{2, 2, 3}, [3]uint{2, 3, 3}, + [3]uint{3, 0, 7}, [3]uint{2, 6, 1}, [3]uint{3, 5, 7}, [3]uint{3, 7, 7}, + ), + want: []morton.Z{42, 59}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + ix := newSimplePointIndex(tt.targetLevel, 1.0) + got := ix.findIntersectingTilesLeft(tt.x, tt.y, tt.targetLevel, tt.classification) + t.Logf("findIntersectingTilesLeft(%d, %d, %d) = %v", tt.x, tt.y, tt.targetLevel, got) + assert.Equal(t, tt.want, got) + }) + } +} + +func squarePolygon(minX, minY, maxX, maxY float64) geom.Polygon { + return geom.Polygon{ + {{minX, minY}, {maxX, minY}, {maxX, maxY}, {minX, maxY}}, + } +} + +func squareWithHolePolygon(minX, minY, maxX, maxY, holeMinX, holeMinY, holeMaxX, holeMaxY float64) geom.Polygon { + return geom.Polygon{ + {{minX, minY}, {maxX, minY}, {maxX, maxY}, {minX, maxY}}, + {{holeMinX, holeMinY}, {holeMinX, holeMaxY}, {holeMaxX, holeMaxY}, {holeMaxX, holeMinY}}, + } +} + +func trianglePolygon(x1, y1, x2, y2, x3, y3 float64) geom.Polygon { + return geom.Polygon{ + {{x1, y1}, {x2, y2}, {x3, y3}}, + } +} + +func TestPointIndex_classifyNonIntersectingTile(t *testing.T) { + tests := []struct { + name string + ix *PointIndex + polygon geom.Polygon + tmsID tms20.TMID + buffer uint + tileLevel Level + tileX, tileY uint + want TileClassification + }{ + { + name: "tile outside a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 0, + tileY: 0, + want: ClassificationOutside, + }, + { + name: "tile outside a centered square", + ix: newSimplePointIndex(3, 1.5), + polygon: squarePolygon(3, 3, 9, 9), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 0, + tileY: 0, + want: ClassificationOutside, + }, + { + name: "tile inside a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 3, + tileY: 3, + want: ClassificationInside, + }, + { + name: "higher-level tile inside a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(1, 1, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 2, + tileX: 1, + tileY: 1, + want: ClassificationInside, + }, + { + name: "tile clearly outside, opposite corner", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 7, + tileY: 7, + want: ClassificationOutside, + }, + { + name: "tile inside the hole of a donut polygon", + ix: newSimplePointIndex(3, 1.0), + // This should work with hole 3,3,5,5 but lineIntesects is buggy + polygon: squareWithHolePolygon(0, 0, 7, 7, 3, 3, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 4, + tileY: 4, + want: ClassificationOutside, + }, + { + name: "tile in the solid part of a donut polygon", + ix: newSimplePointIndex(3, 1.0), + polygon: squareWithHolePolygon(0, 0, 8, 8, 3, 3, 5, 5), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 1, + tileY: 1, + want: ClassificationInside, + }, + { + name: "tile in the solid part of a donut polygon, far corner", + ix: newSimplePointIndex(3, 1.0), + polygon: squareWithHolePolygon(0, 0, 8, 8, 3, 3, 6, 6), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 7, + tileY: 7, + want: ClassificationOutside, + }, + { + name: "tile below a diagonal triangle (raycast hits intersection)", + ix: newSimplePointIndex(3, 1.0), + polygon: trianglePolygon(0, 0, 0, 7, 7, 7), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 1, + tileY: 0, + want: ClassificationOutside, + }, + { + name: "tile above a diagonal triangle (raycast hits intersection)", + ix: newSimplePointIndex(3, 1.0), + polygon: trianglePolygon(0, 0, 7, 0, 0, 7), + tmsID: 3, + buffer: 0, + tileLevel: 3, + tileX: 1, + tileY: 7, + want: ClassificationOutside, + }, + { + name: "nonzero buffer around a centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 4, + tileLevel: 3, + tileX: 0, + tileY: 0, + want: ClassificationOutside, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + segments, classification := tt.ix.registerPolygonEdges(tt.polygon, tt.tmsID, tt.buffer) + z := morton.MustToZ(tt.tileX, tt.tileY) + got := tt.ix.classifyNonIntersectingTile(z, tt.tileLevel, segments, classification, tt.polygon) + t.Logf("classifyNonIntersectingTile(tile=(%d,%d)) = %v", tt.tileX, tt.tileY, got) + assert.Equal(t, tt.want, got) + }) + } +} + +// convert classification data to visual test data +func classificationGrid(classification map[Level]map[morton.Z]TileClassification, targetLevel Level) map[Level][][]TileClassification { + grid := make(map[Level][][]TileClassification, targetLevel+1) + for l := Level(0); l <= targetLevel; l++ { + levelClassification := classification[l] + size := uint(1) << l + rows := make([][]TileClassification, size) + for y := range size { + row := make([]TileClassification, size) + for x := range size { + z := morton.MustToZ(x, y) + if c, ok := levelClassification[z]; ok { + row[x] = c + } else { + row[x] = ClassificationUnknown + } + } + rows[y] = row + } + grid[l] = rows + } + return grid +} + +func TestPointIndex_classifyNonIntersectingTiles(t *testing.T) { + const ( + u = ClassificationUnknown + x = ClassificationIntersect + i = ClassificationInside + o = ClassificationOutside + ) + tests := []struct { + name string + ix *PointIndex + polygon geom.Polygon + tmsID tms20.TMID + buffer uint + want map[Level][][]TileClassification + }{ + { + name: "centered square", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(2, 2, 6, 6), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {o, x, x, o}, + {x, x, x, x}, + {x, x, x, x}, + {o, x, x, x}, + }, + 3: { + {u, u, o, o, o, o, u, u}, + {u, u, x, x, x, x, u, u}, + {o, x, x, x, x, x, x, o}, + {o, x, x, i, i, x, x, o}, + {o, x, x, i, i, x, x, o}, + {o, x, x, x, x, x, x, o}, + {u, u, x, x, x, x, x, o}, + {u, u, o, o, o, o, o, o}, + }, + }, + }, + { + name: "donut polygon (square with a square hole)", + ix: newSimplePointIndex(3, 1.0), + polygon: squareWithHolePolygon(0, 0, 7, 7, 3, 3, 5, 5), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {x, x, x, x}, + {x, x, x, x}, + {x, x, x, x}, + {x, x, x, x}, + }, + 3: { + {x, x, x, x, x, x, x, x}, + {x, i, i, i, i, i, x, x}, + {x, i, i, x, x, i, x, x}, + {x, i, x, x, x, x, x, x}, + {x, i, x, x, x, x, x, x}, + {x, i, i, x, x, x, x, x}, + {x, x, x, x, x, x, x, x}, + {x, x, x, x, x, x, x, x}, + }, + }, + }, + { + name: "diagonal triangle", + ix: newSimplePointIndex(3, 1.0), + polygon: trianglePolygon(0, 0, 7, 0, 7, 7), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {x, x, x, x}, + {x, x, x, x}, + {o, x, x, x}, + {o, o, x, x}, + }, + 3: { + {x, x, x, x, x, x, x, x}, + {x, x, x, i, i, i, x, x}, + {o, x, x, x, i, i, x, x}, + {o, o, x, x, x, i, x, x}, + {u, u, o, x, x, x, x, x}, + {u, u, o, o, x, x, x, x}, + {u, u, u, u, o, x, x, x}, + {u, u, u, u, o, o, o, x}, + }, + }, + }, + { + name: "polygon covering the whole extent", + ix: newSimplePointIndex(2, 1.0), + polygon: squarePolygon(0, 0, 4, 4), + tmsID: 2, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, x}, {x, x}}, + 2: { + {x, x, x, x}, + {x, i, i, x}, + {x, i, i, x}, + {x, x, x, x}, + }, + }, + }, + { + name: "tiny polygon in a single corner tile", + ix: newSimplePointIndex(3, 1.0), + polygon: squarePolygon(0.25, 0.25, 0.75, 0.75), + tmsID: 3, + buffer: 0, + want: map[Level][][]TileClassification{ + 0: {{x}}, + 1: {{x, o}, {o, o}}, + 2: { + {x, o, u, u}, + {o, o, u, u}, + {u, u, u, u}, + {u, u, u, u}, + }, + 3: { + {x, o, u, u, u, u, u, u}, + {o, o, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + {u, u, u, u, u, u, u, u}, + }, + }, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + targetLevel := Level(tt.tmsID) //nolint:gosec // G115 + segments, classification := tt.ix.registerPolygonEdges(tt.polygon, tt.tmsID, tt.buffer) + tt.ix.classifyNonIntersectingTiles(targetLevel, 0, 0, true, segments, classification, tt.polygon) + + got := classificationGrid(classification, targetLevel) + for l := Level(0); l <= targetLevel; l++ { + t.Logf("level %d: %v", l, got[l]) + } + assert.Equal(t, tt.want, got) + }) + } +} diff --git a/pointindex/pointindex.go b/pointindex/pointindex.go index 4760c95..020a9da 100644 --- a/pointindex/pointindex.go +++ b/pointindex/pointindex.go @@ -13,6 +13,7 @@ import ( "github.com/pdok/texel/mathhelp" "github.com/pdok/texel/morton" + "github.com/pdok/texel/tile" "github.com/pdok/texel/tms20" "github.com/go-spatial/geom" @@ -74,11 +75,13 @@ type PointIndex struct { Quadrant deepestLevel Level // Number of quadrants (in one direction) on the deepest level (= 2 ^ deepestLevel) - deepestSize uint - deepestRes intgeom.M - quadrants map[Level]map[morton.Z]Quadrant - hitOnce map[Level]map[intgeom.Point][]int - hitMultiple map[Level]map[intgeom.Point][]int + deepestSize uint + deepestRes intgeom.M + quadrants map[Level]map[morton.Z]Quadrant + hitOnce map[Level]map[intgeom.Point][]int + hitMultiple map[Level]map[intgeom.Point][]int + tilePixels uint + internalPixels uint } type ( @@ -89,7 +92,8 @@ type ( func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID) (*PointIndex, error) { // assuming IsQuadTree was tested before rootTM := tileMatrixSet.TileMatrices[0] - levelDiff := uint(math.Log2(float64(rootTM.TileWidth))) + uint(math.Log2(float64(VectorTileInternalPixelResolution))) + tilePixels := rootTM.TileWidth + levelDiff := uint(math.Log2(float64(tilePixels))) + uint(math.Log2(float64(VectorTileInternalPixelResolution))) //nolint:gosec // G115 deepestLevel := uint(deepestTMID) + levelDiff bottomLeft, topRight, err := tileMatrixSet.MatrixBoundingBox(0) @@ -105,8 +109,10 @@ func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID intExtent: intExtent, z: 0, }, - deepestLevel: deepestLevel, - deepestSize: deepestSize, + tilePixels: tilePixels, + internalPixels: VectorTileInternalPixelResolution, + deepestLevel: deepestLevel, + deepestSize: deepestSize, //nolint:gosec // G115 deepestRes: intExtent.XSpan() / int64(deepestSize), quadrants: make(map[Level]map[morton.Z]Quadrant, deepestLevel+1), @@ -119,7 +125,7 @@ func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID } // Temporary function: primitively marks quadrants as relevant for tiling -func (ix *PointIndex) GetPrimitiveQBBox(l Level) []Quadrant { +func (ix *PointIndex) GetPrimitiveQBBox(l Level) []tile.Tile { quadrants := ix.quadrants[l] minX := ^uint(0) @@ -135,21 +141,75 @@ func (ix *PointIndex) GetPrimitiveQBBox(l Level) []Quadrant { maxY = max(y, maxY) } - quadrantSlice := make([]Quadrant, (maxY-minY+1)*(maxX-minX+1)) + quadrantSlice := make([]tile.Tile, (maxY-minY+1)*(maxX-minX+1)) for i := range maxX - minX + 1 { for j := range maxY - minY + 1 { - extent, centroid := ix.getQuadrantExtentAndCentroid(l, minX+i, minY+j, ix.intExtent) - newQuadrant := Quadrant{ - z: morton.MustToZ(minX+i, minY+j), - intExtent: extent, - intCentroid: centroid, + extent, _ := ix.getQuadrantExtentAndCentroid(l, minX+i, minY+j, ix.intExtent) + newTile := tile.Tile{ + Extent: extent, + X: minX + i, + Y: minY + i, + IsContained: false, } - quadrantSlice[i*(maxY-minY+1)+j] = newQuadrant + quadrantSlice[i*(maxY-minY+1)+j] = newTile } } return quadrantSlice } +// Loop over all points to find the extent. Return slice of tiles whose +// buffer intersects extent. Buufer size given in deepstlevel (internal pixels). +func (ix *PointIndex) GetQBBoxWithBuffer(l Level, bufferSize uint) []tile.Tile { + quadrants := ix.quadrants[ix.deepestLevel] + + minX := ^uint(0) + maxX := uint(0) + minY := ^uint(0) + maxY := uint(0) + + for z := range quadrants { + x, y := morton.FromZ(z) + minX = min(x, minX) + maxX = max(x, maxX) + minY = min(y, minY) + maxY = max(y, maxY) + } + + var tileMinX, tileMinY uint + if minX < bufferSize { + tileMinX = 0 + } else { + tileMinX = (minX - bufferSize) >> (ix.deepestLevel - l) + } + if minY < bufferSize { + tileMinY = 0 + } else { + tileMinY = (minY - bufferSize) >> (ix.deepestLevel - l) + } + + maxTileCoord := mathhelp.Pow2(l) - 1 + tileMaxX := min((maxX+bufferSize)>>(ix.deepestLevel-l), maxTileCoord) + tileMaxY := min((maxY+bufferSize)>>(ix.deepestLevel-l), maxTileCoord) + + tiles := make([]tile.Tile, 0, (tileMaxX-tileMinX+1)*(tileMaxY-tileMinY+1)) + for i := range tileMaxX - tileMinX + 1 { + for j := range tileMaxY - tileMinY + 1 { + tileX := tileMinX + i + tileY := tileMinY + j + extent, _ := ix.getQuadrantExtentAndCentroid( + l, tileX, tileY, ix.intExtent) + + tiles = append(tiles, tile.Tile{ + Extent: extent, + X: tileX, + Y: tileY, + IsContained: false, + }) + } + } + return tiles +} + // InsertPolygon inserts all points from a Polygon func (ix *PointIndex) InsertPolygon(polygon geom.Polygon) error { // initialize the quadrants map diff --git a/pointindex/pointindex_test.go b/pointindex/pointindex_test.go index 40927ea..bc7f7e6 100644 --- a/pointindex/pointindex_test.go +++ b/pointindex/pointindex_test.go @@ -8,6 +8,7 @@ import ( "github.com/pdok/texel/mapslicehelp" "github.com/pdok/texel/mathhelp" "github.com/pdok/texel/morton" + "github.com/pdok/texel/tile" "github.com/stretchr/testify/require" @@ -132,6 +133,74 @@ func TestPointIndex_getQuadrantExtentAndCentroid(t *testing.T) { } } +func TestPointIndex_GetQBBoxWithBuffer(t *testing.T) { + type tileCoord struct{ x, y uint } + tests := []struct { + name string + deepestLevel Level + points [][2]int // Points to be inserted at deepestlevel + tileLevel Level + bufferSize uint + wantTiles []tileCoord + }{ + { + // deepestLevel 4 -> 16x16 pixel grid, level 2 -> 4x4 tiles of 4x4 pixels each. + name: "points all within one tile, no buffer", + deepestLevel: 4, + points: [][2]int{{5, 5}, {6, 6}}, + tileLevel: 2, + bufferSize: 0, + wantTiles: []tileCoord{{1, 1}}, + }, + { + name: "point in one tile, buffer expands to multiple tiles", + deepestLevel: 4, + points: [][2]int{{0, 0}}, + tileLevel: 2, + bufferSize: 6, + wantTiles: []tileCoord{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + { + name: "non-square rectangle of tiles", + deepestLevel: 4, + points: [][2]int{{4, 4}, {11, 5}}, + tileLevel: 2, + bufferSize: 0, + wantTiles: []tileCoord{{1, 1}, {2, 1}}, + }, + { + name: "points in diagonal tiles, expands to rectangle", + deepestLevel: 4, + points: [][2]int{{3, 3}, {4, 4}}, + tileLevel: 2, + bufferSize: 0, + wantTiles: []tileCoord{{0, 0}, {0, 1}, {1, 0}, {1, 1}}, + }, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + ix := newSimplePointIndex(tt.deepestLevel, 1.0) + for _, p := range tt.points { + require.NoError(t, ix.InsertCoord(p[0], p[1])) + } + + want := make([]tile.Tile, 0, len(tt.wantTiles)) + for _, tc := range tt.wantTiles { + extent, _ := ix.getQuadrantExtentAndCentroid(tt.tileLevel, tc.x, tc.y, ix.intExtent) + want = append(want, tile.Tile{ + Extent: extent, + X: tc.x, + Y: tc.y, + IsContained: false, + }) + } + + got := ix.GetQBBoxWithBuffer(tt.tileLevel, tt.bufferSize) + assert.ElementsMatch(t, want, got) + }) + } +} + func TestPointIndex_InsertPoint(t *testing.T) { tests := []struct { name string @@ -148,8 +217,10 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 1.0, 1.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{0.5, 0.5}), }, - deepestLevel: 0, - deepestSize: mathhelp.Pow2(0), + deepestLevel: 0, + deepestSize: mathhelp.Pow2(0), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(1.0) / intgeom.M(mathhelp.Pow2(0)), quadrants: map[Level]map[morton.Z]Quadrant{0: {0: Quadrant{ @@ -167,8 +238,10 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 1.0, 1.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{0.5, 0.5}), }, - deepestLevel: 1, - deepestSize: mathhelp.Pow2(1), + deepestLevel: 1, + deepestSize: mathhelp.Pow2(1), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(1.0) / intgeom.M(mathhelp.Pow2(1)), quadrants: map[Level]map[morton.Z]Quadrant{ @@ -193,8 +266,10 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 4.0, 4.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{2.0, 2.0}), }, - deepestLevel: 3, - deepestSize: mathhelp.Pow2(3), + deepestLevel: 3, + deepestSize: mathhelp.Pow2(3), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(4.0) / intgeom.M(mathhelp.Pow2(3)), quadrants: map[Level]map[morton.Z]Quadrant{ @@ -230,8 +305,10 @@ func TestPointIndex_InsertPoint(t *testing.T) { intExtent: intgeom.FromGeomExtent(geom.Extent{0.0, 0.0, 16.0, 16.0}), intCentroid: intgeom.FromGeomPoint(geom.Point{8.0, 8.0}), }, - deepestLevel: 5, - deepestSize: mathhelp.Pow2(5), + deepestLevel: 5, + deepestSize: mathhelp.Pow2(5), + tilePixels: 256, + internalPixels: 16, //nolint:gosec // G115 deepestRes: intgeom.FromGeomOrd(16.0) / intgeom.M(mathhelp.Pow2(5)), quadrants: map[Level]map[morton.Z]Quadrant{ @@ -509,6 +586,17 @@ func TestPointIndex_lineIntersects(t *testing.T) { }, want: false, }, + { + name: "diagonally opposed line between non-inclusive points", + // This test is fragile and depends on floating-point rounding + extent: intgeom.Extent{ + 0o0000000, 0o0000000, 10000000, 10000000, + }, + line: intgeom.Line{ + {0o0000000, 10000000}, {10000000, 0o0000000}, + }, + want: false, // TODO This is an undesired outcome + }, } for _, tt := range tests { t.Run(tt.name, func(t *testing.T) { @@ -532,10 +620,12 @@ func newSimplePointIndex(deepestLevel Level, cellSize float64) *PointIndex { deepestLevel: deepestLevel, deepestSize: deepestSize, //nolint:gosec // G115 - deepestRes: intExtent.XSpan() / int64(deepestSize), - quadrants: make(map[Level]map[morton.Z]Quadrant, deepestLevel+1), - hitOnce: make(map[morton.Z]map[intgeom.Point][]int, 0), - hitMultiple: make(map[morton.Z]map[intgeom.Point][]int, 0), + deepestRes: intExtent.XSpan() / int64(deepestSize), + quadrants: make(map[Level]map[morton.Z]Quadrant, deepestLevel+1), + hitOnce: make(map[morton.Z]map[intgeom.Point][]int, 0), + hitMultiple: make(map[morton.Z]map[intgeom.Point][]int, 0), + tilePixels: 256, + internalPixels: 16, } _, ix.intCentroid = ix.getQuadrantExtentAndCentroid(0, 0, 0, ix.intExtent) return &ix diff --git a/processing/interface.go b/processing/interface.go index 7db78a7..4345b29 100644 --- a/processing/interface.go +++ b/processing/interface.go @@ -2,7 +2,6 @@ package processing import ( "github.com/go-spatial/geom" - "github.com/pdok/texel/pointindex" "github.com/pdok/texel/tile" ) @@ -21,7 +20,7 @@ type FeatureForTileMatrix interface { type SnapResult struct { Geometry geom.Geometry // Tiles is nil if encoding is not desired - Tiles []pointindex.Quadrant + Tiles []tile.Tile } type TileCoord struct { diff --git a/processing/processing.go b/processing/processing.go index b107200..49c5886 100644 --- a/processing/processing.go +++ b/processing/processing.go @@ -116,19 +116,19 @@ func (stats *countStats) countGeometry(g geom.Geometry) { } } -func encodeGeometry(s SnapResult) []tile.EncodedGeometry { +func encodeGeometry(s SnapResult, defaultEnc tile.DefaultEncoding) []tile.EncodedGeometry { encGeoms := make([]tile.EncodedGeometry, len(s.Tiles)) orig := s.Geometry for i, q := range s.Tiles { - encGeoms[i] = tile.MvtEncodeGeometry(q, orig) + encGeoms[i] = tile.MvtEncodeGeometry(q, orig, defaultEnc) } return encGeoms } // processFeatures processes the geometries in the features with the given function -func processFeatures(featuresIn <-chan Feature, featuresOut chan<- FeatureForTileMatrix, tmIDs []tms20.TMID, f processPolygonFunc, encodeTiles bool) { +func processFeatures(featuresIn <-chan Feature, featuresOut chan<- FeatureForTileMatrix, tmIDs []tms20.TMID, f processPolygonFunc, encodeTiles bool, defaultEnc tile.DefaultEncoding) { stats := initStats() for { feature, hasMore := <-featuresIn @@ -146,7 +146,7 @@ func processFeatures(featuresIn <-chan Feature, featuresOut chan<- FeatureForTil for tmID, snapResult := range newGeometriesPerTileMatrix { var encGeoms []tile.EncodedGeometry if encodeTiles { - encGeoms = encodeGeometry(snapResult) + encGeoms = encodeGeometry(snapResult, defaultEnc) } featuresOut <- wrapFeatureForTileMatrix(feature, tmID, snapResult.Geometry, encGeoms) } @@ -205,7 +205,7 @@ func writeFeaturesToTargets(featuresForTileMatrices <-chan FeatureForTileMatrix, type processPolygonFunc func(p geom.Polygon, tileMatrixIDs []tms20.TMID) map[tms20.TMID]SnapResult // ProcessFeatures applies the processing function/operation to each Target. -func ProcessFeatures(source Source, targets map[tms20.TMID]Target, f processPolygonFunc, encodeTiles bool) { +func ProcessFeatures(source Source, targets map[tms20.TMID]Target, f processPolygonFunc, encodeTiles bool, buffer uint) { featuresBefore := make(chan Feature) featuresAfter := make(chan FeatureForTileMatrix) tileMatrixIDs := make([]tms20.TMID, 0, len(targets)) @@ -213,13 +213,23 @@ func ProcessFeatures(source Source, targets map[tms20.TMID]Target, f processPoly tileMatrixIDs = append(tileMatrixIDs, tmID) } + // When tile is contained in polygon, use default geometry (tile-filling square). + var defaultEnc tile.DefaultEncoding + if encodeTiles { + var err error + defaultEnc, err = tile.NewDefaultEncoding(buffer) + if err != nil { + panic(err) + } + } + wg := sync.WaitGroup{} wg.Add(1) go func() { defer wg.Done() writeFeaturesToTargets(featuresAfter, targets) }() - go processFeatures(featuresBefore, featuresAfter, tileMatrixIDs, f, encodeTiles) + go processFeatures(featuresBefore, featuresAfter, tileMatrixIDs, f, encodeTiles, defaultEnc) go readFeaturesFromSource(source, featuresBefore) wg.Wait() diff --git a/processing/processing_test.go b/processing/processing_test.go index f866ca9..17ebfa1 100644 --- a/processing/processing_test.go +++ b/processing/processing_test.go @@ -5,7 +5,7 @@ import ( "testing" "github.com/go-spatial/geom" - "github.com/pdok/texel/pointindex" + "github.com/pdok/texel/tile" "github.com/pdok/texel/tms20" ) @@ -14,8 +14,8 @@ import ( func TestProcessGeometry(t *testing.T) { polyA := geom.Polygon{{{0, 0}, {1, 0}, {1, 1}, {0, 0}}} polyB := geom.Polygon{{{2, 2}, {3, 2}, {3, 3}, {2, 2}}} - tileA := pointindex.Quadrant{} - tileB := pointindex.Quadrant{} + tileA := tile.Tile{} + tileB := tile.Tile{} tests := []struct { name string @@ -32,10 +32,10 @@ func TestProcessGeometry(t *testing.T) { geometry: polyA, tmIDs: []tms20.TMID{1}, callResults: []map[tms20.TMID]SnapResult{ - {1: {Geometry: polyA, Tiles: []pointindex.Quadrant{tileA}}}, + {1: {Geometry: polyA, Tiles: []tile.Tile{tileA}}}, }, want: map[tms20.TMID]SnapResult{ - 1: {Geometry: polyA, Tiles: []pointindex.Quadrant{tileA}}, + 1: {Geometry: polyA, Tiles: []tile.Tile{tileA}}, }, }, { @@ -69,10 +69,10 @@ func TestProcessGeometry(t *testing.T) { geometry: geom.MultiPolygon{polyA}, tmIDs: []tms20.TMID{1}, callResults: []map[tms20.TMID]SnapResult{ - {1: {Geometry: polyA, Tiles: []pointindex.Quadrant{tileA}}}, + {1: {Geometry: polyA, Tiles: []tile.Tile{tileA}}}, }, want: map[tms20.TMID]SnapResult{ - 1: {Geometry: polyA, Tiles: []pointindex.Quadrant{tileA}}, + 1: {Geometry: polyA, Tiles: []tile.Tile{tileA}}, }, }, { @@ -80,11 +80,11 @@ func TestProcessGeometry(t *testing.T) { geometry: geom.MultiPolygon{polyA, polyB}, tmIDs: []tms20.TMID{1}, callResults: []map[tms20.TMID]SnapResult{ - {1: {Geometry: polyA, Tiles: []pointindex.Quadrant{tileA}}}, - {1: {Geometry: polyB, Tiles: []pointindex.Quadrant{tileB}}}, + {1: {Geometry: polyA, Tiles: []tile.Tile{tileA}}}, + {1: {Geometry: polyB, Tiles: []tile.Tile{tileB}}}, }, want: map[tms20.TMID]SnapResult{ - 1: {Geometry: geom.MultiPolygon{polyA, polyB}, Tiles: []pointindex.Quadrant{tileA, tileB}}, + 1: {Geometry: geom.MultiPolygon{polyA, polyB}, Tiles: []tile.Tile{tileA, tileB}}, }, }, } diff --git a/snap/snap.go b/snap/snap.go index ae6fd0e..d4576eb 100644 --- a/snap/snap.go +++ b/snap/snap.go @@ -11,6 +11,7 @@ import ( "github.com/pdok/texel/geomhelp" "github.com/pdok/texel/mapslicehelp" "github.com/pdok/texel/pointindex" + "github.com/pdok/texel/tile" "github.com/tobshub/go-sortedmap" "golang.org/x/exp/maps" //nolint:exptostd @@ -35,6 +36,11 @@ type Config struct { IgnoreOutsideGrid bool ReverseWindingOrder bool EncodeTiles bool + // Buffer is the number of internal pixels that tiles get + // inflated for detecting which geometries lie on them + Buffer uint + // Decide whether to use lineTrace or BBox for tile detection + UseLineTrace bool } // SnapPolygon snaps polygons' points to a tile's internal pixel grid @@ -68,9 +74,14 @@ func SnapPolygon(polygon geom.Polygon, tileMatrixSet tms20.TileMatrixSet, tmIDs newPolygonsPerTileMatrixID := make(map[tms20.TMID]processing.SnapResult, len(newPolygonsPerLevel)) for level, newPolygons := range newPolygonsPerLevel { - var tilesbbox []pointindex.Quadrant + var tilesbbox []tile.Tile if config.EncodeTiles { - tilesbbox = ix.GetPrimitiveQBBox(pointindex.Level(tmIDsByLevels[level])) //nolint:gosec // G115 These are numbers < 40 + tmID := tmIDsByLevels[level] + if config.UseLineTrace { + tilesbbox = ix.GetLineTraceResult(polygon, tmID, config.Buffer) + } else { + tilesbbox = ix.GetQBBoxWithBuffer(pointindex.Level(tmID), config.Buffer) //nolint:gosec // G115 These are numbers < 40 + } } newGeometry := geomhelp.PolygonSliceToGeom(newPolygons) newPolygonsPerTileMatrixID[tmIDsByLevels[level]] = processing.SnapResult{Geometry: newGeometry, Tiles: tilesbbox} diff --git a/tile/assemble.go b/tile/assemble.go index da38601..78106a7 100644 --- a/tile/assemble.go +++ b/tile/assemble.go @@ -9,8 +9,8 @@ import ( "github.com/go-spatial/geom/encoding/mvt" vectorTile "github.com/go-spatial/geom/encoding/mvt/vector_tile" oldproto "github.com/golang/protobuf/proto" //nolint:staticcheck // needed: matches the reflection-based marshaling vector_tile.pb.go relies on + "github.com/pdok/texel/intgeom" "github.com/pdok/texel/mapslicehelp" - "github.com/pdok/texel/pointindex" ) const ( @@ -32,30 +32,64 @@ type EncodedFeatureRow struct { Geom EncodedGeometry } -// Transform geometry to tile extent, then encode. We assume the geometry is -// snapped to the proposed grid, in which case makevalid operations should not -// be necessary. -func MvtEncodeGeometry(q pointindex.Quadrant, g geom.Geometry) EncodedGeometry { - ext := q.Extent() - preparedGeo := mvt.PrepareGeo(g, &ext, float64(precision)) +type Tile struct { + Extent intgeom.Extent + X uint + Y uint + IsContained bool +} + +type DefaultEncoding struct { + Encoding []uint32 + GeomType vectorTile.Tile_GeomType +} + +// Default tile-filling geometry with buffer +func NewDefaultEncoding(buffer uint) (DefaultEncoding, error) { + fbuffer := float64(buffer) + defaultPolygon := geom.Polygon{{ + {-fbuffer, -fbuffer}, + {-fbuffer, precision + fbuffer}, + {precision + fbuffer, precision + fbuffer}, + {precision + fbuffer, -fbuffer}, + }} + + defaultExtent := geom.Extent{ + 0, 0, precision, precision, + } + preparedPolygon := mvt.PrepareGeo(defaultPolygon, &defaultExtent, precision) - // This should not be necessary. - // sg, err := convert.ToTegola(preparedGeo) - // tegolaGeo, err := validate.CleanGeometry(context.TODO(), sg, &ext) - // validatedGeo := convert.ToGeom(tegolaGeo) + encoding, geomType, err := EncodeGeometry(preparedPolygon) + return DefaultEncoding{Encoding: encoding, GeomType: geomType}, err +} - encgeom, geomtype, err := EncodeGeometry(preparedGeo) +// Encode geometry. defaultEnc is the precomputed tile-filling geometry +func MvtEncodeGeometry(t Tile, g geom.Geometry, defaultEnc DefaultEncoding) EncodedGeometry { + var encgeom []uint32 + var geomtype vectorTile.Tile_GeomType + var err error + if t.IsContained { + encgeom, geomtype = defaultEnc.Encoding, defaultEnc.GeomType + } else { + ext := t.Extent.ToGeomExtent() + preparedGeo := mvt.PrepareGeo(g, &ext, float64(precision)) + + // This should not be necessary. + // sg, err := convert.ToTegola(preparedGeo) + // tegolaGeo, err := validate.CleanGeometry(context.TODO(), sg, &ext) + // validatedGeo := convert.ToGeom(tegolaGeo) + + encgeom, geomtype, err = EncodeGeometry(preparedGeo) + } if err != nil { panic(err) } - xTile, yTile := q.Coords() - return EncodedGeometry{ Encoding: encgeom, GeometryType: int32(geomtype), - XTile: xTile, - YTile: yTile, + XTile: t.X, + YTile: t.Y, } }