--- a/cmd/srsimplify/main.go	Wed Jun 26 12:06:28 2019 +0200
+++ b/cmd/srsimplify/main.go	Mon Jul 01 11:08:19 2019 +0200
@@ -19,7 +19,6 @@
 	"fmt"
 	"io"
 	"log"
-	"math"
 	"os"
 	"strconv"
 	"strings"
@@ -70,175 +69,6 @@
 	return out.Flush()
 }
 
-func simplify(points octree.MultiPointZ, tolerance float64) octree.MultiPointZ {
-
-	if len(points) < 2 {
-		return points
-	}
-
-	if tolerance < 0 {
-		tolerance = -tolerance
-	}
-
-	min := octree.Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}
-	max := octree.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 := []octree.Vertex{
-		{xMin, yMin, below},
-		{xMax, yMin, below},
-		{xMax, yMax, below},
-		{xMin, yMax, below},
-	}
-
-	top := points[maxIdx]
-
-	tris := make([]octree.Triangle, len(corners))
-	planes := make([]octree.Plane3D, len(corners))
-
-	for i, v1 := range corners {
-		v2 := corners[(i+1)%len(corners)]
-		tris[i] = octree.Triangle{v1, v2, top}
-		planes[i] = tris[i].Plane3D()
-	}
-
-	parts := make([][]octree.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([]octree.Vertex, 0, len(points))
-
-	var handleTriangle func(*octree.Triangle, float64, int, []octree.Vertex) bool
-
-	handleTriangle = func(
-		t *octree.Triangle,
-		maxDist float64, maxIdx int,
-		points []octree.Vertex,
-	) bool {
-		if maxDist <= tolerance {
-			return false
-		}
-
-		if len(points) == 1 {
-			result = append(result, points[0])
-		}
-
-		var (
-			tris     [3]octree.Triangle
-			planes   [3]octree.Plane3D
-			maxDists [3]float64
-			maxIdxs  [3]int
-			parts    [3][]octree.Vertex
-		)
-
-		top := points[maxIdx]
-		for i := 0; i < 3; i++ {
-			tris[i] = octree.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
-}
-
 func main() {
 
 	var tolerance float64
@@ -253,7 +83,7 @@
 		log.Fatalf("err: %v\n", err)
 	}
 
-	points = simplify(points, tolerance)
+	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:08:19 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
+}