--- a/pkg/octree/strtree.go	Mon Oct 14 13:18:15 2019 +0200
+++ b/pkg/octree/strtree.go	Mon Oct 14 14:58:04 2019 +0200
@@ -32,12 +32,72 @@
 	bboxes  []Box2D
 }
 
+// Vertical does a vertical cross cut from (x1, y1) to (x2, y2).
+func (s *STRTree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {
+	box := Box2D{
+		X1: math.Min(x1, x2),
+		Y1: math.Min(y1, y2),
+		X2: math.Max(x1, x2),
+		Y2: math.Max(y1, y2),
+	}
+
+	// out of bounding box
+	if box.X2 < s.tin.Min.X || s.tin.Max.X < box.X1 ||
+		box.Y2 < s.tin.Min.Y || s.tin.Max.Y < box.Y1 {
+		return
+	}
+
+	line := NewPlane2D(x1, y1, x2, y2)
+
+	vertices := s.tin.Vertices
+
+	stack := []int32{s.index[0]}
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		if top > 0 { // node
+			if b := s.bbox(top); b.Intersects(box) && b.IntersectsPlane(line) {
+				n := s.index[top+1]
+				stack = append(stack, s.index[top+2:top+2+n]...)
+			}
+		} else { // leaf
+			top = -top - 1
+			if b := s.bbox(top); !b.Intersects(box) || !b.IntersectsPlane(line) {
+				continue
+			}
+			n := s.index[top+1]
+			for _, idx := range s.index[top+2 : top+2+n] {
+				ti := s.tin.Triangles[idx]
+				t := Triangle{
+					vertices[ti[0]],
+					vertices[ti[1]],
+					vertices[ti[2]],
+				}
+				v0 := line.Eval(t[0].X, t[0].Y)
+				v1 := line.Eval(t[1].X, t[1].Y)
+				v2 := line.Eval(t[2].X, t[2].Y)
+
+				if onPlane(v0) || onPlane(v1) || onPlane(v2) ||
+					sides(sides(sides(0, v0), v1), v2) == 3 {
+					fn(&t)
+				}
+			}
+		}
+	}
+}
+
 func (s *STRTree) Min() Vertex  { return s.tin.Min }
 func (s *STRTree) Max() Vertex  { return s.tin.Max }
 func (s *STRTree) EPSG() uint32 { return s.tin.EPSG }
 
 func (s *STRTree) Value(x, y float64) (float64, bool) {
 
+	if len(s.index) == 0 {
+		return 0, false
+	}
+
 	stack := []int32{s.index[0]}
 
 	vertices := s.tin.Vertices
@@ -56,8 +116,8 @@
 			if !s.bbox(top).Contains(x, y) {
 				continue
 			}
-			for i, n := int32(0), s.index[top+1]; i < n; i++ {
-				idx := s.index[top+2+i]
+			n := s.index[top+1]
+			for _, idx := range s.index[top+2 : top+2+n] {
 				ti := s.tin.Triangles[idx]
 				t := Triangle{
 					vertices[ti[0]],
@@ -412,3 +472,51 @@
 	}
 	return int32(-(pos + 1))
 }
+
+func (s *STRTree) Diff(other *STRTree) PointMap {
+
+	firstVs, secondVs := s.tin.Vertices, other.tin.Vertices
+
+	result := make(PointMap, len(firstVs)+len(secondVs))
+
+	sliceWork(
+		firstVs,
+		result,
+		func(slice []Vertex, turn func([]Vertex) []Vertex) {
+			p := turn(nil)
+			for i := range slice {
+				v := &slice[i]
+				if z, found := other.Value(v.X, v.Y); found {
+					p = append(p, Vertex{v.X, v.Y, v.Z - z})
+					if len(p) == cap(p) {
+						p = turn(p)
+					}
+				}
+			}
+			if len(p) > 0 {
+				turn(p)
+			}
+		})
+
+	sliceWork(
+		secondVs,
+		result,
+		func(
+			slice []Vertex, turn func([]Vertex) []Vertex) {
+			p := turn(nil)
+			for i := range slice {
+				v := &slice[i]
+				if z, found := s.Value(v.X, v.Y); found {
+					p = append(p, Vertex{v.X, v.Y, z - v.Z})
+					if len(p) == cap(p) {
+						p = turn(p)
+					}
+				}
+			}
+			if len(p) > 0 {
+				turn(p)
+			}
+		})
+
+	return result
+}