--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/strtree.go	Tue Mar 05 13:05:50 2019 +0100
@@ -0,0 +1,194 @@
+// This is Free Software under GNU Affero General Public License v >= 3.0
+// without warranty, see README.md and license for details.
+//
+// SPDX-License-Identifier: AGPL-3.0-or-later
+// License-Filename: LICENSES/AGPL-3.0.txt
+//
+// Copyright (C) 2018 by via donau
+//   – Österreichische Wasserstraßen-Gesellschaft mbH
+// Software engineering by Intevation GmbH
+//
+// Author(s):
+//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
+
+package octree
+
+import (
+	"math"
+	"sort"
+)
+
+const numEntries = 8
+
+type STRTree struct {
+	tri       *Triangulation
+	index     []int32
+	triangles [][]int32
+	bboxes    []Box2D
+}
+
+func (s *STRTree) Build(t *Triangulation) {
+
+	s.tri = t
+
+	s.triangles = t.TriangleSlices()
+
+	all := make([]int32, len(s.triangles))
+
+	for i := range all {
+		all[i] = int32(i)
+	}
+
+	s.index = append(s.index, 0)
+
+	root := s.build(all)
+
+	s.index[0] = root
+}
+
+func (s *STRTree) build(items []int32) int32 {
+	sort.Slice(items, func(i, j int) bool {
+		ti := s.triangles[items[i]]
+		tj := s.triangles[items[j]]
+		return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X
+	})
+
+	P := int(math.Ceil(float64(len(items)) / float64(numEntries)))
+	S := int(math.Ceil(math.Sqrt(float64(P))))
+
+	sm := S * numEntries
+
+	slices := make([][]int32, S)
+	for i := range slices {
+		end := (i + 1) * sm
+		if end > len(items) {
+			end = len(items)
+		}
+		slices[i] = items[i*sm : end]
+	}
+
+	leaves := make([]int32, 0, S*S)
+
+	for _, slice := range slices {
+		sort.Slice(slice, func(i, j int) bool {
+			ti := s.triangles[slice[i]]
+			tj := s.triangles[slice[j]]
+			return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X
+		})
+
+		for len(slice) > 0 {
+			n := numEntries
+			if l := len(slice); l < numEntries {
+				n = l
+			}
+			leaves = append(leaves, s.allocLeaf(slice[:n]))
+			slice = slice[n:]
+		}
+	}
+
+	return s.buildNodes(leaves)
+}
+
+func bboxIndex(idx int32) int32 {
+	if idx < 0 {
+		return -idx - 1
+	}
+	return idx
+}
+
+func (s *STRTree) buildNodes(items []int32) int32 {
+
+	if len(items) <= numEntries {
+		return s.allocNode(items)
+	}
+
+	sort.Slice(items, func(i, j int) bool {
+		bi := s.bboxes[s.index[bboxIndex(items[i])]]
+		bj := s.bboxes[s.index[bboxIndex(items[j])]]
+		return bi.X1 < bj.X1
+	})
+
+	P := int(math.Ceil(float64(len(items)) / float64(numEntries)))
+	S := int(math.Ceil(math.Sqrt(float64(P))))
+
+	sm := S * numEntries
+
+	slices := make([][]int32, S)
+	for i := range slices {
+		end := (i + 1) * sm
+		if end > len(items) {
+			end = len(items)
+		}
+		slices[i] = items[i*sm : end]
+	}
+
+	nodes := make([]int32, 0, S*S)
+
+	for _, slice := range slices {
+		sort.Slice(slice, func(i, j int) bool {
+			bi := s.bboxes[s.index[bboxIndex(slice[i])]]
+			bj := s.bboxes[s.index[bboxIndex(slice[j])]]
+			return bi.X1 < bj.X1
+		})
+
+		for len(slice) > 0 {
+			n := numEntries
+			if l := len(slice); l < numEntries {
+				n = l
+			}
+			nodes = append(nodes, s.allocNode(slice[:n]))
+			slice = slice[n:]
+		}
+	}
+
+	return s.buildNodes(nodes)
+}
+
+func (s *STRTree) allocNode(items []int32) int32 {
+	pos := len(s.index)
+	s.index = append(s.index, 0)
+	s.index = append(s.index, items...)
+	if len(items) > 0 {
+		box := s.bboxes[s.index[bboxIndex(items[0])]]
+		for i := 1; i < len(items); i++ {
+			it := items[i]
+			b := s.bboxes[s.index[bboxIndex(it)]]
+			box = box.Union(b)
+		}
+		s.index[pos] = int32(s.allocBBox(box))
+	}
+	return int32(pos)
+}
+
+func (s *STRTree) allocBBox(box Box2D) int {
+	pos := len(s.bboxes)
+	s.bboxes = append(s.bboxes, box)
+	return pos
+}
+
+func (s *STRTree) allocLeaf(items []int32) int32 {
+	pos := len(s.index)
+	s.index = append(s.index, 0)
+	s.index = append(s.index, items...)
+	if len(items) > 0 {
+		ti := s.triangles[items[0]]
+		t := Triangle{
+			s.tri.Points[ti[0]],
+			s.tri.Points[ti[1]],
+			s.tri.Points[ti[2]],
+		}
+		box := t.BBox()
+		for i := 1; i < len(items); i++ {
+			it := items[i]
+			ti := s.triangles[it]
+			t := Triangle{
+				s.tri.Points[ti[0]],
+				s.tri.Points[ti[1]],
+				s.tri.Points[ti[2]],
+			}
+			box = box.Union(t.BBox())
+		}
+		s.index[pos] = int32(s.allocBBox(box))
+	}
+	return -int32(pos) - 1
+}

--- a/pkg/octree/vertex.go	Mon Mar 04 18:16:13 2019 +0100
+++ b/pkg/octree/vertex.go	Tue Mar 05 13:05:50 2019 +0100
@@ -723,6 +723,15 @@
 	return a.Y2
 }
 
+func (a Box2D) Union(b Box2D) Box2D {
+	return Box2D{
+		X1: math.Min(a.X1, b.X1),
+		Y1: math.Min(a.Y1, b.Y1),
+		X2: math.Max(a.X2, b.X2),
+		Y2: math.Max(a.Y2, b.Y2),
+	}
+}
+
 // NewPlane2D creates a new Plane2D from two given points.
 func NewPlane2D(x1, y1, x2, y2 float64) Plane2D {
 	b := x2 - x1