--- a/pkg/mesh/vertex.go	Wed Feb 14 16:56:28 2024 +0100
+++ b/pkg/mesh/vertex.go	Tue Feb 20 21:28:56 2024 +0100
@@ -88,7 +88,24 @@
 
 	v0 := t[1].Sub(t[0])
 	v1 := t[2].Sub(t[0])
-	n := v0.Cross(v1).Normalize()
+	n := v0.Cross(v1)
+
+	// If the length of normal is near zero assume we have
+	// a plane constant in z with an average z value of
+	// the three vertices.
+	// This should protect us from the effects of collinear
+	// geometries.
+	l := n.Length()
+	if l < 1e-7 {
+		sum := t[0].Z + t[1].Z + t[2].Z
+		return Plane3D{
+			A: 0,
+			B: 0,
+			C: 1,
+			D: -sum / 3,
+		}
+	}
+	n = n.Scale(1 / l)
 
 	// x*nx+ y*ny+ z*nz + d = 0
 	// d = - (x*nx+ y*ny + z*nz)
@@ -1144,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 {