--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/vertex.go	Sat Sep 22 21:57:30 2018 +0200
@@ -0,0 +1,195 @@
+package octree
+
+import (
+	"bytes"
+	"encoding/binary"
+	"io"
+	"math"
+)
+
+type (
+	Vertex struct {
+		X float64
+		Y float64
+		Z float64
+	}
+
+	Triangle [3]Vertex
+
+	Line [2]Vertex
+
+	LineStringZ      []Vertex
+	MultiLineStringZ []LineStringZ
+)
+
+const (
+	wkbNDR              byte   = 1
+	wkbPointZ           uint32 = 1000 + 1
+	wkbLineStringZ      uint32 = 1000 + 2
+	wkbMultiLineStringZ uint32 = 1000 + 5
+)
+
+func (v *Vertex) Minimize(w Vertex) {
+	if w.X < v.X {
+		v.X = w.X
+	}
+	if w.Y < v.Y {
+		v.Y = w.Y
+	}
+	if w.Z < v.Z {
+		v.Z = w.Z
+	}
+}
+
+func (v *Vertex) Maximize(w Vertex) {
+	if w.X > v.X {
+		v.X = w.X
+	}
+	if w.Y > v.Y {
+		v.Y = w.Y
+	}
+	if w.Z > v.Z {
+		v.Z = w.Z
+	}
+}
+
+func (v Vertex) Sub(w Vertex) Vertex {
+	return Vertex{
+		v.X - w.X,
+		v.Y - w.Y,
+		v.Z - w.Z,
+	}
+}
+
+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 {
+		return Vertex{
+			v2.X*s.X + v1.X,
+			v2.Y*s.Y + v1.Y,
+			v2.Z*s.Z + v1.Z,
+		}
+	}
+}
+
+func (a Vertex) Less(b Vertex) bool {
+	return a.X < b.X || a.Y < b.Y || a.Z < b.Z
+}
+
+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) IntersectHorizontal(h float64) Vertex {
+	t := (h - l[1].Z) / l[0].Z
+	return l.Eval(t)
+}
+
+func side(z, h float64) int {
+	switch {
+	case z < h:
+		return -1
+	case z > h:
+		return +1
+	}
+	return 0
+}
+
+func (t *Triangle) IntersectHorizontal(h float64) LineStringZ {
+	sides := [3]int{
+		side(t[0].Z, h),
+		side(t[1].Z, h),
+		side(t[2].Z, h),
+	}
+
+	var points LineStringZ
+
+	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]).IntersectHorizontal(h)
+			points = append(points, v)
+		}
+	}
+
+	return points
+}
+
+func (v *Vertex) read(r io.Reader) error {
+	var buf [8]byte
+	b := buf[:]
+	if _, err := io.ReadFull(r, b); err != nil {
+		return nil
+	}
+	v.X = math.Float64frombits(binary.LittleEndian.Uint64(b))
+	if _, err := io.ReadFull(r, b); err != nil {
+		return nil
+	}
+	v.Y = math.Float64frombits(binary.LittleEndian.Uint64(b))
+	if _, err := io.ReadFull(r, b); err != nil {
+		return nil
+	}
+	v.Z = math.Float64frombits(binary.LittleEndian.Uint64(b))
+	return nil
+}
+
+func (mls MultiLineStringZ) AsWKB() []byte {
+
+	var buf bytes.Buffer
+
+	binary.Write(&buf, binary.LittleEndian, wkbNDR)
+	binary.Write(&buf, binary.LittleEndian, wkbMultiLineStringZ)
+	binary.Write(&buf, binary.LittleEndian, uint32(len(mls)))
+
+	for _, ml := range mls {
+		binary.Write(&buf, binary.LittleEndian, wkbNDR)
+		binary.Write(&buf, binary.LittleEndian, wkbLineStringZ)
+		binary.Write(&buf, binary.LittleEndian, uint32(len(ml)))
+		for _, p := range ml {
+			binary.Write(&buf, binary.LittleEndian, math.Float64bits(p.X))
+			binary.Write(&buf, binary.LittleEndian, math.Float64bits(p.Y))
+			binary.Write(&buf, binary.LittleEndian, math.Float64bits(p.Z))
+		}
+	}
+
+	return buf.Bytes()
+}