--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/srsimplify/main.go	Mon Jul 01 11:10:47 2019 +0200
@@ -0,0 +1,91 @@
+// 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) 2019 by via donau
+//   – Österreichische Wasserstraßen-Gesellschaft mbH
+// Software engineering by Intevation GmbH
+//
+// Author(s):
+//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
+
+package main
+
+import (
+	"bufio"
+	"flag"
+	"fmt"
+	"io"
+	"log"
+	"os"
+	"strconv"
+	"strings"
+
+	"gemma.intevation.de/gemma/pkg/octree"
+)
+
+func loadXYZ(r io.Reader) (octree.MultiPointZ, error) {
+
+	scanner := bufio.NewScanner(r)
+
+	points := make(octree.MultiPointZ, 0, 2000000)
+
+	var x, y, z float64
+	var err error
+
+	for scanner.Scan() {
+		line := strings.TrimSpace(scanner.Text())
+		if len(line) == 0 || strings.HasPrefix(line, "#") {
+			continue
+		}
+		parts := strings.SplitN(line, " ", 3)
+		if len(parts) != 3 {
+			continue
+		}
+
+		if x, err = strconv.ParseFloat(parts[0], 64); err != nil {
+			return nil, err
+		}
+		if y, err = strconv.ParseFloat(parts[1], 64); err != nil {
+			return nil, err
+		}
+		if z, err = strconv.ParseFloat(parts[2], 64); err != nil {
+			return nil, err
+		}
+		points = append(points, octree.Vertex{X: x, Y: y, Z: z})
+	}
+
+	return points, nil
+}
+
+func storeXYZ(points octree.MultiPointZ, w io.Writer) error {
+	out := bufio.NewWriter(w)
+	for i := range points {
+		fmt.Fprintf(out, "%.5f,%.5f,%.5f\n",
+			points[i].X, points[i].Y, points[i].Z)
+	}
+	return out.Flush()
+}
+
+func main() {
+
+	var tolerance float64
+
+	flag.Float64Var(&tolerance, "t", 0.1, "accepted tolerance (shorthand)")
+	flag.Float64Var(&tolerance, "tolerance", 0.1, "accepted tolerance")
+
+	flag.Parse()
+
+	points, err := loadXYZ(os.Stdin)
+	if err != nil {
+		log.Fatalf("err: %v\n", err)
+	}
+
+	points = points.Simplify(tolerance)
+
+	if err := storeXYZ(points, os.Stdout); err != nil {
+		log.Fatalf("err: %v\n", err)
+	}
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/simplify.go	Mon Jul 01 11:10:47 2019 +0200
@@ -0,0 +1,189 @@
+// 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) 2019 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"
+)
+
+func (points MultiPointZ) Simplify(tolerance float64) MultiPointZ {
+
+	if len(points) < 2 {
+		return points
+	}
+
+	if tolerance < 0 {
+		tolerance = -tolerance
+	}
+
+	min := Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}
+	max := Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}
+
+	var maxIdx int
+
+	for i, v := range points {
+		min.Minimize(v)
+
+		if v.X < min.X {
+			min.X = v.X
+		}
+		if v.X > max.X {
+			max.X = v.X
+		}
+		if v.Y < min.Y {
+			min.Y = v.Y
+		}
+		if v.Y > max.Y {
+			max.Y = v.Y
+		}
+		if v.Z < min.Z {
+			min.Z = v.Z
+		}
+		if v.Z > max.Z {
+			max.Z = v.Z
+			maxIdx = i
+		}
+
+		max.Maximize(v)
+	}
+
+	/*
+		log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n",
+			min.X, min.Y, min.Z,
+			max.X, max.Y, max.Z)
+	*/
+
+	below := min.Z - 3*tolerance
+	xMin := min.X - tolerance
+	xMax := max.X + tolerance
+	yMin := min.Y - tolerance
+	yMax := max.Y + tolerance
+
+	corners := []Vertex{
+		{xMin, yMin, below},
+		{xMax, yMin, below},
+		{xMax, yMax, below},
+		{xMin, yMax, below},
+	}
+
+	top := points[maxIdx]
+
+	tris := make([]Triangle, len(corners))
+	planes := make([]Plane3D, len(corners))
+
+	for i, v1 := range corners {
+		v2 := corners[(i+1)%len(corners)]
+		tris[i] = Triangle{v1, v2, top}
+		planes[i] = tris[i].Plane3D()
+	}
+
+	parts := make([][]Vertex, len(tris))
+
+	maxDists := make([]float64, len(planes))
+	maxIdxs := make([]int, len(planes))
+
+nextPoint:
+	for i, v := range points {
+		if i == maxIdx {
+			continue
+		}
+
+		for j := range tris {
+			if tris[j].Contains(v.X, v.Y) {
+				if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
+					maxDists[j] = dist
+					maxIdxs[j] = len(parts[j])
+				}
+				parts[j] = append(parts[j], v)
+				continue nextPoint
+			}
+		}
+	}
+
+	result := make([]Vertex, 0, len(points))
+
+	var handleTriangle func(*Triangle, float64, int, []Vertex) bool
+
+	handleTriangle = func(
+		t *Triangle,
+		maxDist float64, maxIdx int,
+		points []Vertex,
+	) bool {
+		if maxDist <= tolerance {
+			return false
+		}
+
+		if len(points) == 1 {
+			result = append(result, points[0])
+		}
+
+		var (
+			tris     [3]Triangle
+			planes   [3]Plane3D
+			maxDists [3]float64
+			maxIdxs  [3]int
+			parts    [3][]Vertex
+		)
+
+		top := points[maxIdx]
+		for i := 0; i < 3; i++ {
+			tris[i] = Triangle{t[i], t[(i+1)%3], top}
+			planes[i] = tris[i].Plane3D()
+		}
+
+	nextPoint:
+		for i, v := range points {
+			if i == maxIdx {
+				continue
+			}
+
+			for j := range tris {
+				if tris[j].Contains(v.X, v.Y) {
+					if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
+						maxDists[j] = dist
+						maxIdxs[j] = len(parts[j])
+					}
+					parts[j] = append(parts[j], v)
+					continue nextPoint
+				}
+			}
+		}
+
+		var found bool
+		for i, part := range parts {
+			if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
+				found = true
+			}
+		}
+
+		if found {
+			result = append(result, top)
+		}
+
+		return found
+	}
+
+	var found bool
+	for i, part := range parts {
+		if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
+			found = true
+		}
+	}
+
+	if found {
+		result = append(result, top)
+	}
+
+	return result
+}

--- a/pkg/octree/vertex.go	Thu Jun 27 16:27:48 2019 +0200
+++ b/pkg/octree/vertex.go	Mon Jul 01 11:10:47 2019 +0200
@@ -110,6 +110,10 @@
 	return -(p.A*x + p.B*y + p.D) / p.C
 }
 
+func (p Plane3D) Eval(v Vertex) float64 {
+	return p.A*v.X + p.B*v.Y + p.C*v.Z + p.D
+}
+
 func (v Vertex) Normalize() Vertex {
 	s := 1 / v.Length()
 	return Vertex{