add tests for DeviationStats and IsQuadTree

PDOK-16504

main.go

@@ -174,7 +174,7 @@ func validateTileMatrixSet(tms tms20.TileMatrixSet, tileMatrixIDs []tms20.TMID)
 		log.Printf("warning, (largest) deviation is larger than 1 tile pixel (%f units) on the deepest matrix (%d)\n", deviationInUnits, deepestTMID)
 		log.Println(stats)
 	}
-	return nil
+	return pointindex.IsQuadTree(tms)
 }
 
 func initGPKGTarget(targetPathFmt string, tmID int, overwrite bool, pagesize int) *gpkg.TargetGeopackage {

mathhelp/mathhelp.go

@@ -1,6 +1,13 @@
 package mathhelp
 
-func BetweenInc(f, p, q int64) bool {
+func IBetweenInc(f, p, q int64) bool {
+	if p <= q {
+		return p <= f && f <= q
+	}
+	return q <= f && f <= p
+}
+
+func FBetweenInc(f, p, q float64) bool {
 	if p <= q {
 		return p <= f && f <= q
 	}

pointindex/pointindex.go

@@ -3,12 +3,13 @@ package pointindex
 import (
 	"errors"
 	"fmt"
-	"golang.org/x/exp/maps"
 	"io"
 	"math"
 	"slices"
 	"strconv"
 
+	"golang.org/x/exp/maps"
+
 	"github.com/pdok/texel/mathhelp"
 	"github.com/pdok/texel/morton"
 	"github.com/pdok/texel/tms20"
@@ -75,16 +76,14 @@ type PointIndex struct {
 type Level = uint
 type Q = int // quadrant index (0, 1, 2 or 3)
 
-func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID) *PointIndex {
-	if err := isQuadTree(tileMatrixSet); err != nil {
-		panic(err)
-	}
+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)))
 	deepestLevel := uint(deepestTMID) + levelDiff
 	bottomLeft, topRight, err := tileMatrixSet.MatrixBoundingBox(0)
 	if err != nil {
-		panic(fmt.Errorf(`could not make PointIndex from TileMatrixSet %v: %w'`, tileMatrixSet.ID, err))
+		return nil, fmt.Errorf(`could not make PointIndex from TileMatrixSet %v: %w'`, tileMatrixSet.ID, err)
 	}
 	intBottomLeft := intgeom.FromGeomPoint(bottomLeft)
 	intTopRight := intgeom.FromGeomPoint(topRight)
@@ -104,7 +103,7 @@ func FromTileMatrixSet(tileMatrixSet tms20.TileMatrixSet, deepestTMID tms20.TMID
 	}
 	_, ix.intCentroid = ix.getQuadrantExtentAndCentroid(0, 0, 0, intExtent)
 
-	return &ix
+	return &ix, nil
 }
 
 // InsertPolygon inserts all points from a Polygon
@@ -471,7 +470,7 @@ func lineOverlapsInclusiveEdge(intLine intgeom.Line, edgeI int, intEdge intgeom.
 	exclusiveTip := getExclusiveTip(edgeI, intEdge)
 	lOrd1 := intLine[0][varAx]
 	lOrd2 := intLine[1][varAx]
-	return lOrd1 != lOrd2 && (mathhelp.BetweenInc(lOrd1, eOrd1, eOrd2) && intLine[0] != exclusiveTip || mathhelp.BetweenInc(lOrd2, eOrd1, eOrd2) && intLine[1] != exclusiveTip)
+	return lOrd1 != lOrd2 && (mathhelp.IBetweenInc(lOrd1, eOrd1, eOrd2) && intLine[0] != exclusiveTip || mathhelp.IBetweenInc(lOrd2, eOrd1, eOrd2) && intLine[1] != exclusiveTip)
 }
 
 func oneIfRight(quadrantI int) int {
@@ -496,7 +495,8 @@ func (ix *PointIndex) ToWkt(writer io.Writer) {
 	}
 }
 
-func isQuadTree(tms tms20.TileMatrixSet) error {
+//nolint:nestif
+func IsQuadTree(tms tms20.TileMatrixSet) error {
 	var previousTMID int
 	var previousTM *tms20.TileMatrix
 	tmIDs := maps.Keys(tms.TileMatrices)
@@ -529,6 +529,15 @@ func isQuadTree(tms tms20.TileMatrixSet) error {
 			if tm.CornerOfOrigin != previousTM.CornerOfOrigin {
 				return errors.New("tile matrixes should have the same corner of origin: " + tm.ID)
 			}
+			if tm.TileHeight != previousTM.TileHeight {
+				return errors.New("tile matrix tiles should stay the same size: " + tm.ID)
+			}
+			if tm.MatrixHeight != 2*previousTM.MatrixHeight {
+				return errors.New("tile matrix should double in size each level: " + tm.ID)
+			}
+			if !mathhelp.FBetweenInc(previousTM.CellSize/tm.CellSize, 1.99, 2.01) { // between because of fp error
+				return errors.New("cell size should half each level: " + tm.ID)
+			}
 		}
 
 		previousTMID = tmID
@@ -546,7 +555,10 @@ func DeviationStats(tms tms20.TileMatrixSet, deepestTMID tms20.TMID) (stats stri
 	if err != nil {
 		return
 	}
-	ix := FromTileMatrixSet(tms, deepestTMID)
+	ix, err := FromTileMatrixSet(tms, deepestTMID)
+	if err != nil {
+		return
+	}
 	p := uint(intgeom.Precision + 1)
 	ps := strconv.Itoa(int(p))
 	stats += fmt.Sprintf("deepest level: %d\n", ix.deepestLevel)

pointindex/pointindex_test.go

@@ -1,6 +1,7 @@
 package pointindex
 
 import (
+	"fmt"
 	"os"
 	"testing"
 
@@ -19,6 +20,10 @@ import (
 	"github.com/stretchr/testify/assert"
 )
 
+func assertNoErr(t assert.TestingT, err error, _ ...any) bool {
+	return assert.Nil(t, err)
+}
+
 func TestPointIndex_containsPoint(t *testing.T) {
 	tests := []struct {
 		name string
@@ -320,7 +325,8 @@ func TestPointIndex_InsertPoint_Deepest(t *testing.T) {
 	for _, tt := range tests {
 		t.Run(tt.name, func(t *testing.T) {
 			tms := loadEmbeddedTileMatrixSet(t, tt.tmsID)
-			ix := FromTileMatrixSet(tms, tt.tmID)
+			ix, err := FromTileMatrixSet(tms, tt.tmID)
+			require.Nil(t, err)
 
 			ix.InsertPoint(tt.point)
 			assert.Equal(t, 1, len(ix.quadrants[ix.deepestLevel]))
@@ -534,5 +540,113 @@ func loadEmbeddedTileMatrixSet(t *testing.T, tmsID string) tms20.TileMatrixSet {
 }
 
 func newPointIndexFromEmbeddedTileMatrixSet(t *testing.T, tmsID string, deepestTMID tms20.TMID) *PointIndex {
-	return FromTileMatrixSet(loadEmbeddedTileMatrixSet(t, tmsID), deepestTMID)
+	tms, err := FromTileMatrixSet(loadEmbeddedTileMatrixSet(t, tmsID), deepestTMID)
+	require.Nil(t, err)
+	return tms
+}
+
+func TestIsQuadTree(t *testing.T) {
+	tests := []struct {
+		name    string
+		tms     tms20.TileMatrixSet
+		wantErr assert.ErrorAssertionFunc
+	}{
+		{
+			name:    "NetherlandsRDNewQuad",
+			tms:     loadEmbeddedTileMatrixSet(t, "NetherlandsRDNewQuad"),
+			wantErr: assertNoErr,
+		}, {
+			name:    "WebMercatorQuad",
+			tms:     loadEmbeddedTileMatrixSet(t, "WebMercatorQuad"),
+			wantErr: assertNoErr,
+		}, {
+			name:    "EuropeanETRS89_LAEAQuad",
+			tms:     loadEmbeddedTileMatrixSet(t, "EuropeanETRS89_LAEAQuad"),
+			wantErr: assertNoErr,
+		}, {
+			name: "GNOSISGlobalGrid",
+			tms:  loadEmbeddedTileMatrixSet(t, "GNOSISGlobalGrid"),
+			wantErr: func(t assert.TestingT, err error, _ ...any) bool {
+				return assert.ErrorContains(t, err, "tile matrix height should be same as width")
+			},
+		}, {
+			name: "LINZAntarticaMapTilegrid",
+			tms:  loadEmbeddedTileMatrixSet(t, "LINZAntarticaMapTilegrid"),
+			wantErr: func(t assert.TestingT, err error, _ ...any) bool {
+				return assert.ErrorContains(t, err, "tile matrix should double in size each level")
+			},
+		}, {
+			name:    "WorldMercatorWGS84Quad",
+			tms:     loadEmbeddedTileMatrixSet(t, "WorldMercatorWGS84Quad"),
+			wantErr: assertNoErr,
+		},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			tt.wantErr(t, IsQuadTree(tt.tms), fmt.Sprintf("IsQuadTree(%v)", tt.tms))
+		})
+	}
+}
+
+func TestDeviationStats(t *testing.T) {
+	tests := []struct {
+		name                  string
+		tms                   tms20.TileMatrixSet
+		deepestTMID           tms20.TMID
+		wantStats             string
+		wantDeviationInUnits  float64
+		wantDeviationInPixels float64
+		margin                float64
+		wantErr               assert.ErrorAssertionFunc
+	}{
+		{
+			name:                  "NetherlandsRDNewQuad",
+			tms:                   loadEmbeddedTileMatrixSet(t, "NetherlandsRDNewQuad"),
+			deepestTMID:           16,
+			wantDeviationInUnits:  0,
+			wantDeviationInPixels: 0,
+			margin:                1e-6, // micrometers
+			wantErr:               assertNoErr,
+		},
+		{
+			name:                  "WebMercatorQuad",
+			tms:                   loadEmbeddedTileMatrixSet(t, "WebMercatorQuad"),
+			deepestTMID:           18,
+			wantDeviationInUnits:  0,
+			wantDeviationInPixels: 0,
+			margin:                1,
+			wantErr:               assertNoErr,
+		},
+		{
+			name:                  "WebMercatorQuad starting from 19 has more than 1 pixel deviation ... ;(",
+			tms:                   loadEmbeddedTileMatrixSet(t, "WebMercatorQuad"),
+			deepestTMID:           19,
+			wantDeviationInUnits:  1,
+			wantDeviationInPixels: 6,
+			margin:                1,
+			wantErr:               assertNoErr,
+		},
+		{
+			name:                  "EuropeanETRS89_LAEAQuad",
+			tms:                   loadEmbeddedTileMatrixSet(t, "EuropeanETRS89_LAEAQuad"),
+			deepestTMID:           15,
+			wantDeviationInUnits:  0,
+			wantDeviationInPixels: 0,
+			margin:                1,
+			wantErr:               assertNoErr,
+		},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			gotStats, gotDeviationInUnits, gotDeviationInPixels, err := DeviationStats(tt.tms, tt.deepestTMID)
+			if !tt.wantErr(t, err, fmt.Sprintf("DeviationStats(%v, %v)", tt.tms.ID, tt.deepestTMID)) {
+				return
+			}
+			if tt.wantStats != "" {
+				assert.Containsf(t, tt.wantStats, gotStats, "DeviationStats(%v, %v)", tt.tms.ID, tt.deepestTMID)
+			}
+			assert.True(t, mathhelp.FBetweenInc(gotDeviationInUnits, tt.wantDeviationInUnits-tt.margin, tt.wantDeviationInUnits+tt.margin), "DeviationStats(%v, %v)", tt.tms.ID, tt.deepestTMID)
+			assert.True(t, mathhelp.FBetweenInc(gotDeviationInPixels, tt.wantDeviationInPixels-tt.margin, tt.wantDeviationInPixels+tt.margin), "DeviationStats(%v, %v)", tt.tms.ID, tt.deepestTMID)
+		})
+	}
 }

snap/snap.go

@@ -35,7 +35,10 @@ type IsOuter = bool
 //nolint:revive
 func SnapPolygon(polygon geom.Polygon, tileMatrixSet tms20.TileMatrixSet, tmIDs []tms20.TMID, keepPointsAndLines bool) map[tms20.TMID][]geom.Polygon {
 	deepestID := slices.Max(tmIDs)
-	ix := pointindex.FromTileMatrixSet(tileMatrixSet, deepestID)
+	ix, err := pointindex.FromTileMatrixSet(tileMatrixSet, deepestID)
+	if err != nil {
+		panic(err) // TODO let processing.processPolygonFunc return err
+	}
 	tmIDsByLevels := tileMatrixIDsByLevels(tileMatrixSet, tmIDs)
 	levels := make([]pointindex.Level, 0, len(tmIDsByLevels))
 	for level := range tmIDsByLevels {

snap/snap_test.go

@@ -1,6 +1,7 @@
 package snap
 
 import (
+	"strconv"
 	"testing"
 
 	"github.com/pdok/texel/geomhelp"
@@ -928,7 +929,7 @@ func newSimpleTileMatrixSet(deepestTMID pointindex.Level, cellSize float64) tms2
 		// (only values from the root tm are used, for the rest it is assumed to follow quad matrix rules)
 		tmCellSize := cellSize * float64(mathhelp.Pow2(deepestTMID-uint(tmID)))
 		tms.TileMatrices[tmID] = tms20.TileMatrix{
-			ID:               "0",
+			ID:               strconv.Itoa(tmID),
 			ScaleDenominator: tmCellSize / tms20.StandardizedRenderingPixelSize,
 			CellSize:         tmCellSize,
 			CornerOfOrigin:   tms20.BottomLeft,