--- a/pkg/octree/vertex.go	Wed Sep 26 17:04:47 2018 +0200
+++ b/pkg/octree/vertex.go	Wed Sep 26 18:36:12 2018 +0200
@@ -402,3 +402,213 @@
 
 	return p.X / p.Z, p.Y / p.Z, true
 }
+
+type VerticalLine struct {
+	x1 float64
+	y1 float64
+	x2 float64
+	y2 float64
+
+	line Plane2D
+}
+
+func NewVerticalLine(x1, y1, x2, y2 float64) *VerticalLine {
+	return &VerticalLine{
+		x1:   x1,
+		y1:   y1,
+		x2:   x2,
+		y2:   y2,
+		line: NewPlane2D(x1, y1, x2, y2),
+	}
+}
+
+func onPlane(x float64) bool { return math.Abs(x) < epsPlane }
+
+func dist(x1, y1, x2, y2 float64) float64 {
+	return math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))
+}
+
+func interpolate(v1, v2 Vertex) func(x, y float64) float64 {
+
+	length := dist(v1.X, v1.Y, v2.X, v2.Y)
+
+	// f(0) = v1.Z
+	// f(length) = v2.Z
+	// v1.Z = 0*m + b
+	// v2.Z = length*m + b
+
+	b := v1.Z
+	m := (v2.Z - v1.Z) / length
+
+	return func(x, y float64) float64 {
+		l := dist(v1.X, v1.Y, x, y)
+		return m*l + b
+	}
+}
+
+func linear(v1, v2 Vertex) func(t float64) Vertex {
+	return func(t float64) Vertex {
+		return Vertex{
+			(v2.X-v1.X)*t + v1.X,
+			(v2.Y-v1.Y)*t + v1.Y,
+			(v2.Z-v1.Z)*t + v1.Z,
+		}
+	}
+}
+
+func inRange(a, b, c float64) bool {
+	if a > b {
+		a, b = b, a
+	}
+	return c >= a && c <= b
+}
+
+func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ {
+
+	var out LineStringZ
+
+	length := dist(vl.x1, vl.y1, vl.x2, vl.y2)
+
+edges:
+	for i := 0; i < 3 && len(out) < 2; i++ {
+		j := (i + 1) % 3
+		edge := NewPlane2D(t[i].X, t[i].Y, t[j].X, t[j].Y)
+
+		s1 := edge.Eval(vl.x1, vl.y1)
+		s2 := edge.Eval(vl.x2, vl.y2)
+
+		o1 := onPlane(s1)
+		o2 := onPlane(s2)
+
+		switch {
+		case o1 && o2:
+			// XXX: Broken!
+			pos := interpolate(t[i], t[j])
+			t1 := pos(vl.x1, vl.y1)
+			t2 := pos(vl.x2, vl.y2)
+
+			r1 := inRange(0, length, t1)
+			r2 := inRange(0, length, t2)
+
+			switch {
+			case r1 && r2:
+				l := linear(t[i], t[j])
+				// XXX: Broken!
+				out = append(out, l(t1/length), l(t2/length))
+				return out
+
+			case !r1 && !r2: // the hole edge
+				out = append(out, t[i], t[j])
+				return out
+			case !r1:
+				// TODO
+			case !r2:
+				// TODO
+			}
+
+		case o1:
+			t1 := interpolate(t[i], t[j])(vl.x1, vl.y1)
+			if !inRange(0, length, t1) {
+				continue edges
+			}
+			out = append(out, linear(t[i], t[j])(t1))
+
+		case o2:
+			t2 := interpolate(t[i], t[j])(vl.x2, vl.y2)
+			if !inRange(0, length, t2) {
+				continue edges
+			}
+			out = append(out, linear(t[i], t[j])(t2))
+
+		default:
+			x, y, intersects := vl.line.Intersection(edge)
+			if !intersects {
+				continue edges
+			}
+			t1 := interpolate(t[i], t[j])(x, y)
+			if !inRange(0, length, t1) {
+				continue edges
+			}
+
+			// XXX: Ugly!
+			z := (t[j].Z-t[i].Z)*(t1/length) + t[1].Z
+			n := Vertex{x, y, z}
+
+			if math.Signbit(s1) != math.Signbit(s2) {
+				// different sides -> vertex on edge.
+				out = append(out)
+			} else { // same side -> inside. Find peer
+				if len(out) > 0 { // we already have our peer
+					out = append(out, n)
+					return out
+				}
+
+				for k := 0; k < 3; k++ {
+					if k == i {
+						continue
+					}
+					l := (k + 1) % 3
+					other := NewPlane2D(t[k].X, t[k].Y, t[l].X, t[l].Y)
+					xo, yo, intersects := vl.line.Intersection(edge)
+					if !intersects {
+						continue
+					}
+					// XXX: Broken!!!
+					zo := interpolate(t[k], t[j])(xo, yo)
+					if !inRange(t[k].Z, t[k].Z, zo) {
+						continue
+					}
+					m := Vertex{xo, yo, zo}
+
+					var xn, yn, xf, yf float64
+
+					// Which is nearer to current edge?
+					if math.Abs(s1) < math.Abs(s2) {
+						xn, yn = vl.x1, vl.y1
+						xf, yf = vl.x2, vl.y2
+					} else {
+						xn, yn = vl.x2, vl.y2
+						xf, yf = vl.x1, vl.y1
+					}
+
+					if onPlane(other.Eval(xn, yn)) {
+						// triangle intersect in only point
+						// treat as no intersection
+						return out
+					}
+
+					// XXX: Broken!!!
+					zint := interpolate(n, m)
+
+					zn := zint(xn, yn)
+					zf := zint(xf, yf)
+
+					if !inRange(n.Z, m.Z, zn) {
+						// if near is out of range far is, too.
+						return out
+					}
+
+					if inRange(n.Z, m.Z, zf) {
+						m.X = xf
+						m.Y = yf
+						m.Z = zf
+					} // else m is clipping
+
+					n.X = xn
+					n.Y = yn
+					n.Z = zn
+
+					out = append(out, n, m)
+					return out
+				}
+			}
+		}
+	}
+
+	// supress single point touches.
+	if len(out) == 1 {
+		out = out[:0]
+	}
+
+	return out
+}