# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1708429929 -3600
# Node ID 148abae1fcd085f38b61ddfcc061d8cf8cd947e0
# Parent  39d91e76c05c0eab99f57464612f588e5573a11a
Quantize xyz in points in X/Y to 1/QuantScale.

diff -r 39d91e76c05c -r 148abae1fcd0 pkg/imports/sr.go
--- a/pkg/imports/sr.go	Tue Feb 20 09:51:20 2024 +0100
+++ b/pkg/imports/sr.go	Tue Feb 20 12:52:09 2024 +0100
@@ -512,9 +512,14 @@
 	if err := xyz.FromWKB(reproj); err != nil {
 		return nil, err
 	}
-
 	feedback.Info("Reprojecting points to EPSG %d took %v.",
 		epsg, time.Since(start))
+
+	start = time.Now()
+	xyz = xyz.QuantizeXY(mesh.QuantScale)
+	feedback.Info("Quantizing points in X/Y to 1/%d took %v.",
+		mesh.QuantScale, time.Since(start))
+
 	feedback.Info("Number of reprojected points: %d", len(xyz))
 	feedback.Info("Triangulate XYZ data.")
 
diff -r 39d91e76c05c -r 148abae1fcd0 pkg/mesh/meshserialize_v2.go
--- a/pkg/mesh/meshserialize_v2.go	Tue Feb 20 09:51:20 2024 +0100
+++ b/pkg/mesh/meshserialize_v2.go	Tue Feb 20 12:52:09 2024 +0100
@@ -24,10 +24,10 @@
 	"gemma.intevation.de/gemma/pkg/common"
 )
 
-const quantScale = 1_000
+const QuantScale = 1_000
 
-func quant(x float64) int64   { return int64(math.Round(x * quantScale)) }
-func unquant(x int64) float64 { return float64(x) / quantScale }
+func quant(x float64) int64   { return int64(math.Round(x * QuantScale)) }
+func unquant(x int64) float64 { return float64(x) / QuantScale }
 
 func (s *STRTree) serializeV2(w io.Writer) error {
 	return serializer(w, s,
diff -r 39d91e76c05c -r 148abae1fcd0 pkg/mesh/vertex.go
--- a/pkg/mesh/vertex.go	Tue Feb 20 09:51:20 2024 +0100
+++ b/pkg/mesh/vertex.go	Tue Feb 20 12:52:09 2024 +0100
@@ -1161,6 +1161,39 @@
 	return out
 }
 
+// QuantizeXY quantize the X/Y values to scale and
+// removes duplicates in place. The z value of
+// duplicates will be averaged pairwise.
+func (mpz MultiPointZ) QuantizeXY(scale float64) MultiPointZ {
+	type qpoint struct {
+		x int64
+		y int64
+	}
+	m := make(map[qpoint]float64, len(mpz))
+	for _, p := range mpz {
+		k := qpoint{
+			x: int64(math.Round(p.X * scale)),
+			y: int64(math.Round(p.Y * scale)),
+		}
+		if z, ok := m[k]; ok {
+			m[k] = (z + p.Z) * 0.5
+		} else {
+			m[k] = p.Z
+		}
+	}
+	i := 0
+	invScale := 1 / scale
+	for k, z := range m {
+		mpz[i] = Vertex{
+			X: float64(k.x) * invScale,
+			Y: float64(k.y) * invScale,
+			Z: z,
+		}
+		i++
+	}
+	return mpz[:i]
+}
+
 // Filter returns a copy removed the vertices which
 // don't pass the filter test.
 func (mpz MultiPointZ) Filter(filter func(Vertex) bool) MultiPointZ {