--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/vertex.go	Tue Sep 18 00:11:32 2018 +0200
@@ -0,0 +1,107 @@
+package main
+
+import "math"
+
+type vertex struct {
+	x float64
+	y float64
+	z float64
+}
+
+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
+	}
+}
+
+type plane struct {
+	a float64
+	b float64
+	c float64
+	d float64
+}
+
+func (a vertex) cross(b vertex) vertex {
+	return vertex{
+		a.y*b.z - a.z*b.y,
+		a.z*b.x - a.x*b.z,
+		a.x*b.y - a.y*b.x,
+	}
+}
+
+func (a vertex) sub(b vertex) vertex {
+	return vertex{
+		a.x - b.x,
+		a.y - b.y,
+		a.z - b.z,
+	}
+}
+
+func (a vertex) scale(s float64) vertex {
+	return vertex{
+		s * a.x,
+		s * a.y,
+		s * a.z,
+	}
+}
+
+func (a vertex) dot(b vertex) float64 {
+	return a.x*b.x + a.y*b.y + a.z*b.z
+}
+
+func (a vertex) length() float64 {
+	return math.Sqrt(a.x*a.x + a.y*a.y + a.z*a.z)
+}
+
+func newPlane(v1, v2, v3 vertex) plane {
+	w1 := v2.sub(v1)
+	w2 := v3.sub(v1)
+	p := w1.cross(w2)
+	p = p.scale(1 / p.length())
+
+	// x*p.x + y*p.y + z*p.z + d = 0
+	// d = -(x*p.x + y*p.y + z*p.z)
+	d := -p.dot(v1)
+	return plane{
+		a: p.x,
+		b: p.y,
+		c: p.z,
+		d: d,
+	}
+}
+
+func (p plane) eval(v vertex) float64 {
+	return p.a*v.x + p.b*v.y + p.c*v.z + p.d
+}
+
+func interpolate(v1, v2 vertex) func(vertex) vertex {
+	return func(s vertex) vertex {
+		return vertex{
+			(v2.x-v1.x)*s.x + v1.x,
+			(v2.y-v1.y)*s.y + v1.y,
+			(v2.z-v1.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
+}