--- a/cmd/octree2contour/vertex.go	Thu Sep 20 00:13:47 2018 +0200
+++ b/cmd/octree2contour/vertex.go	Thu Sep 20 11:32:03 2018 +0200
@@ -38,6 +38,22 @@
 	}
 }
 
+func (v vertex) add(w vertex) vertex {
+	return vertex{
+		v.x + w.x,
+		v.y + w.y,
+		v.z + w.z,
+	}
+}
+
+func (v vertex) scale(s float64) vertex {
+	return vertex{
+		s * v.x,
+		s * v.y,
+		s * v.z,
+	}
+}
+
 func interpolate(v1, v2 vertex) func(vertex) vertex {
 	v2 = v2.sub(v1)
 	return func(s vertex) vertex {
@@ -52,3 +68,67 @@
 func (a vertex) less(b vertex) bool {
 	return a.x < b.x || a.y < b.y || a.z < b.z
 }
+
+type line [2]vertex
+
+func newLine(p1, p2 vertex) line {
+	return line{
+		p2.sub(p1),
+		p1,
+	}
+}
+
+func (l line) eval(t float64) vertex {
+	return l[0].scale(t).add(l[1])
+}
+
+func (l line) intersectH(h float64) vertex {
+	t := (h - l[1].z) / l[0].z
+	return l.eval(t)
+}
+
+type triangle [3]vertex
+
+func side(z, h float64) int {
+	switch {
+	case z < h:
+		return -1
+	case z > h:
+		return +1
+	}
+	return 0
+}
+
+func (t *triangle) intersectH(h float64, points []vertex) []vertex {
+	sides := [3]int{
+		side(t[0].z, h),
+		side(t[1].z, h),
+		side(t[2].z, h),
+	}
+
+	for i := 0; i < 3; i++ {
+		j := (i + 1) % 3
+		si := sides[i]
+		sj := sides[j]
+
+		switch {
+		case si == 0 && sj == 0:
+			// both on plane
+			points = append(points, t[i], t[j])
+		case si == 0 && sj != 0:
+			// first on plane
+			points = append(points, t[i])
+		case si != 0 && sj == 0:
+			// second on plane
+			points = append(points, t[j])
+		case si == sj:
+			// both on same side
+		default:
+			// real intersection
+			v := newLine(t[i], t[j]).intersectH(h)
+			points = append(points, v)
+		}
+	}
+
+	return points
+}