--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/isoareas/main.go	Tue Sep 24 17:19:52 2019 +0200
@@ -0,0 +1,111 @@
+// 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"
+	"io"
+	"log"
+	"os"
+	"strconv"
+	"strings"
+
+	"gemma.intevation.de/gemma/pkg/octree"
+)
+
+const classBreaks = `1:#ff00dd,1.5,1.7,1.9,2.1,2.3,2.5:#f25f20,2.7,2.9,3.1:#f7e40e,3.3,3.5,4:#8ad51a,4.5,5,5.5,6,6.5,7:#1414ff`
+
+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 minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) {
+	if len(points) == 0 {
+		return octree.Vertex{}, octree.Vertex{}
+	}
+
+	min, max := points[0], points[0]
+
+	for _, v := range points {
+		min.Minimize(v)
+		max.Maximize(v)
+	}
+
+	return min, max
+}
+
+func check(err error) {
+	if err != nil {
+		log.Fatalf("error: %v\n", err)
+	}
+}
+
+func main() {
+
+	heights, err := octree.ParseClassBreaks(classBreaks)
+	check(err)
+
+	xyz, err := loadXYZ(os.Stdin)
+	check(err)
+
+	log.Printf("num vertices: %d\n", len(xyz))
+
+	min, max := minMax(xyz)
+
+	log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
+		min.X, min.Y, min.Z,
+		max.X, max.Y, max.Z)
+
+	heights = octree.ExtrapolateClassBreaks(heights, min.Z, max.Z)
+
+	log.Println("classes:")
+	for i, v := range heights {
+		log.Printf("\t%d: %.3f\n", i+1, v)
+	}
+
+	tri, err := octree.Triangulate(xyz)
+	check(err)
+
+	_ = tri
+}