--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/octree2contour/octree.go	Fri Sep 21 11:02:24 2018 +0200
@@ -0,0 +1,81 @@
+package main
+
+import (
+	"math"
+)
+
+type octree struct {
+	epsg uint32
+
+	vertices  []vertex
+	triangles [][]int32
+	index     []int32
+
+	min vertex
+	max vertex
+}
+
+func (ot *octree) horizontal(h float64, fn func(*triangle)) {
+
+	type frame struct {
+		pos int32
+		min float64
+		max float64
+	}
+
+	if h < ot.min.z || ot.max.z < h {
+		return
+	}
+
+	stack := []frame{{1, ot.min.z, ot.max.z}}
+
+	dupes := map[int32]struct{}{}
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		pos := top.pos
+		if pos == 0 {
+			continue
+		}
+		min, max := top.min, top.max
+
+		if pos > 0 { // node
+			if mid := (max-min)*0.5 + min; h >= mid {
+				pos += 4 // nodes with z-bit set
+				min = mid
+			} else {
+				max = mid
+			}
+			stack = append(stack,
+				frame{ot.index[pos+0], min, max},
+				frame{ot.index[pos+1], min, max},
+				frame{ot.index[pos+2], min, max},
+				frame{ot.index[pos+3], min, max})
+		} else { // leaf
+			pos = -pos - 1
+			n := ot.index[pos]
+			//log.Printf("%d %d %d\n", pos, n, len(ot.index))
+			indices := ot.index[pos+1 : pos+1+n]
+
+			for _, idx := range indices {
+				if _, found := dupes[idx]; found {
+					continue
+				}
+				tri := ot.triangles[idx]
+				t := triangle{
+					ot.vertices[tri[0]],
+					ot.vertices[tri[1]],
+					ot.vertices[tri[2]],
+				}
+
+				if !(math.Min(t[0].z, math.Min(t[1].z, t[2].z)) > h ||
+					math.Max(t[0].z, math.Max(t[1].z, t[2].z)) < h) {
+					dupes[idx] = struct{}{}
+					fn(&t)
+				}
+			}
+		}
+	}
+}