--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/isoareas/algo.go	Tue Oct 01 11:07:33 2019 +0200
@@ -0,0 +1,419 @@
+// 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"
+	"container/list"
+	"log"
+	"math"
+	"math/bits"
+	"os"
+	"runtime"
+	"sync"
+
+	"gemma.intevation.de/gemma/pkg/octree"
+)
+
+type indexedArc struct {
+	arc   int32
+	index int32
+}
+
+func glue(a, b octree.LineStringZ) octree.LineStringZ {
+
+	a1, a2 := a[0], a[len(a)-1]
+	b1, b2 := b[0], b[len(b)-1]
+
+	if a1.EpsEquals2D(b2) {
+		c := make(octree.LineStringZ, len(a)-1+len(b))
+		copy(c, b)
+		copy(c[len(b):], a[1:])
+		return c
+	}
+
+	if a2.EpsEquals2D(b1) {
+		c := make(octree.LineStringZ, len(a)-1+len(b))
+		copy(c, a)
+		copy(c[len(a):], b[1:])
+		return c
+	}
+
+	if a1.EpsEquals2D(b1) {
+		c := make(octree.LineStringZ, len(a)-1+len(b))
+		copy(c, b)
+		c[:len(b)].Reverse()
+		copy(c[len(b):], a[1:])
+		return c
+	}
+
+	if a2.EpsEquals2D(b2) {
+		c := make(octree.LineStringZ, len(a)-1+len(b))
+		copy(c, a)
+		copy(c[len(a):], b[:len(b)-1])
+		c[len(a):].Reverse()
+		return c
+	}
+
+	return nil
+}
+
+func connectArcs(tri *octree.Triangulation, cuts []indexedArc, arcs *[]octree.LineStringZ) {
+
+	unique := map[int32]struct{}{}
+	for i := range cuts {
+		unique[cuts[i].arc] = struct{}{}
+	}
+	before := len(unique)
+
+	origLen := int32(len(*arcs))
+
+	for i := range cuts {
+		if cuts[i].arc >= origLen {
+			// already has a connected arc assigned.
+			continue
+		}
+
+		line := (*arcs)[cuts[i].arc]
+
+		var modified []int
+		visited := map[int32]struct{}{}
+
+		var recursive func(int32)
+		recursive = func(idx int32) {
+			visited[idx] = struct{}{}
+
+			ns := neighbors(tri, idx)
+			for _, n := range ns {
+				n /= 3
+				if _, already := visited[n]; already {
+					continue
+				}
+
+				arcIdx := findArc(n, cuts)
+				if arcIdx == -1 {
+					continue
+				}
+				arc := cuts[arcIdx].arc
+				if arc >= origLen {
+					// already a new arc.
+					continue
+				}
+				oline := (*arcs)[arc]
+
+				nline := glue(line, oline)
+				if nline == nil {
+					// not joinable
+					continue
+				}
+				line = nline
+				modified = append(modified, arcIdx)
+				recursive(cuts[arcIdx].index)
+			}
+		}
+		recursive(cuts[i].index)
+		if len(modified) > 0 {
+			// alloc a new line an fix the references.
+			nidx := int32(len(*arcs))
+			*arcs = append(*arcs, line)
+			cuts[i].arc = nidx
+			for _, j := range modified {
+				cuts[j].arc = nidx
+			}
+		}
+	}
+
+	unique = map[int32]struct{}{}
+	for i := range cuts {
+		unique[cuts[i].arc] = struct{}{}
+	}
+	log.Printf("unique arcs: before: %d after %d\n",
+		before, len(unique))
+}
+
+func intersectTriangles(
+	tri *octree.Triangulation,
+	heights []float64,
+	min, max octree.Vertex,
+) [][]octree.LineStringZ {
+
+	cuts := make([][]indexedArc, len(heights))
+
+	classes := make([][]int32, len(heights)+1)
+
+	var arcs []octree.LineStringZ
+
+all:
+	for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
+		ti := tri.Triangles[i : i+3]
+		v0 := tri.Points[ti[0]]
+		v1 := tri.Points[ti[1]]
+		v2 := tri.Points[ti[2]]
+		min := math.Min(v0.Z, math.Min(v1.Z, v2.Z))
+
+		if min > heights[len(heights)-1] {
+			classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
+			continue
+		}
+		max := math.Max(v0.Z, math.Max(v1.Z, v2.Z))
+
+		for j, h := range heights {
+
+			var l float64
+			if j > 0 {
+				l = heights[j-1]
+			} else {
+				l = -math.MaxFloat64
+			}
+
+			if l < min && max < h {
+				classes[j] = append(classes[j], i/3)
+				continue all
+			}
+			if min > h || max < h {
+				continue
+			}
+			t := octree.Triangle{v0, v1, v2}
+			cut := t.IntersectHorizontal(h)
+			if len(cut) >= 2 {
+				arc := int32(len(arcs))
+				arcs = append(arcs, cut)
+				cuts[j] = append(cuts[j], indexedArc{arc: arc, index: i / 3})
+			}
+		}
+	}
+
+	// connect the arcs in a cut list to longer arcs.
+
+	for _, c := range cuts {
+		connectArcs(tri, c, &arcs)
+	}
+
+	func() {
+		out, err := os.Create("arcs.svg")
+		if err != nil {
+			log.Printf("err: %v\n", err)
+			return
+		}
+		defer func() {
+			out.Close()
+			log.Println("writing arcs done")
+		}()
+
+		buf := bufio.NewWriter(out)
+
+		dumpArcsSVG(
+			buf,
+			min, max,
+			cuts,
+			arcs)
+		buf.Flush()
+	}()
+
+	result := make([][]octree.LineStringZ, len(classes))
+
+	jobs := make(chan int)
+
+	var wg sync.WaitGroup
+
+	worker := func() {
+		defer wg.Done()
+
+		pb := polygonBuilder{open: list.New()}
+
+		for i := range jobs {
+			usedArcs := map[int32]struct{}{}
+			var dupes int
+			var isolated, inside, found int
+			c := classes[i]
+
+		allInClass:
+			for _, idx := range c {
+				ns := neighbors(tri, idx)
+				mask := where(ns, c)
+				switch bits.OnesCount8(mask) {
+				case 3:
+					// Totally insides do not contribute to the geometries.
+					inside++
+					continue allInClass
+				case 0:
+					// Isolated are areas w/o connections to the rest.
+					isolated++
+					ti := tri.Triangles[idx*3 : idx*3+3]
+					pb.addTriangle(
+						tri.Points[ti[0]],
+						tri.Points[ti[1]],
+						tri.Points[ti[2]])
+					continue allInClass
+				default:
+					ti := tri.Triangles[idx*3 : idx*3+3]
+					for j := 0; j < 3; j++ {
+						if (mask & (1 << j)) == 0 {
+
+							var curr octree.LineStringZ
+
+							for l := i - 1; l <= i; l++ {
+								if l < 0 || l >= len(cuts) {
+									continue
+								}
+								arcIdx := findArc(ns[j]/3, cuts[l])
+								if arcIdx == -1 {
+									continue
+								}
+								found++
+								aIdx := cuts[l][arcIdx].arc
+								if _, already := usedArcs[aIdx]; already {
+									dupes++
+									continue
+								}
+								usedArcs[aIdx] = struct{}{}
+
+								curr = arcs[aIdx]
+								break
+							}
+
+							if curr == nil {
+								k := (j + 1) % 3
+								front := tri.Points[ti[j]]
+								back := tri.Points[ti[k]]
+								curr = octree.LineStringZ{front, back}
+							}
+
+						segment:
+							for e := pb.open.Front(); e != nil; e = e.Next() {
+								line := e.Value.(octree.LineStringZ)
+								if nline := glue(curr, line); nline != nil {
+									curr = nline
+									pb.open.Remove(e)
+									goto segment
+								}
+							} // all open
+
+							// check if closed
+							if curr[0].EpsEquals2D(curr[len(curr)-1]) {
+								if !curr.CCW() {
+									curr.Reverse()
+								}
+								pb.polygons = append(pb.polygons, curr)
+							} else {
+								pb.open.PushFront(curr)
+							}
+						} // for all border parts
+					}
+				}
+			}
+
+			log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n",
+				i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found)
+			/*
+				for e := pb.open.Front(); e != nil; e = e.Next() {
+					line := e.Value.(octree.LineStringZ)
+					if !line.CCW() {
+						line.Reverse()
+					}
+					pb.polygons = append(pb.polygons, line)
+				}
+			*/
+
+			result[i] = pb.polygons
+			pb.reset()
+		}
+	}
+
+	for n := runtime.NumCPU(); n >= 0; n-- {
+		wg.Add(1)
+		go worker()
+	}
+
+	log.Println("inside classes:")
+	for i := range classes {
+		jobs <- i
+	}
+	close(jobs)
+
+	wg.Wait()
+
+	log.Println("cuts:")
+	for i, c := range cuts {
+		log.Printf("\t%.3f: %d\n", heights[i], len(c))
+	}
+
+	return result
+
+	// TODO: sew segments together.
+
+}
+
+type polygonBuilder struct {
+	polygons []octree.LineStringZ
+
+	open *list.List
+}
+
+func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
+	polygon := octree.LineStringZ{v0, v1, v2, v0}
+	pb.polygons = append(pb.polygons, polygon)
+}
+
+func (pb *polygonBuilder) reset() {
+	pb.polygons = pb.polygons[:0]
+	pb.open.Init()
+}
+
+func neighbors(t *octree.Triangulation, idx int32) []int32 {
+	idx *= 3
+	return t.Halfedges[idx : idx+3]
+}
+
+func findArc(needle int32, haystack []indexedArc) int {
+	lo, hi := 0, len(haystack)-1
+	for lo <= hi {
+		mid := (hi-lo)/2 + lo
+		switch v := haystack[mid].index; {
+		case v < needle:
+			lo = mid + 1
+		case v > needle:
+			hi = mid - 1
+		default:
+			return mid
+		}
+	}
+	return -1
+}
+
+func contains(needle int32, haystack []int32) bool {
+	lo, hi := 0, len(haystack)-1
+	for lo <= hi {
+		mid := (hi-lo)/2 + lo
+		switch v := haystack[mid]; {
+		case v < needle:
+			lo = mid + 1
+		case v > needle:
+			hi = mid - 1
+		default:
+			return true
+		}
+	}
+	return false
+}
+
+func where(neighbors, indices []int32) byte {
+	var mask byte
+	for i, n := range neighbors {
+		if n >= 0 && contains(n/3, indices) {
+			mask |= 1 << i
+		}
+	}
+	return mask
+}

--- a/cmd/isoareas/main.go	Mon Sep 30 16:34:55 2019 +0200
+++ b/cmd/isoareas/main.go	Tue Oct 01 11:07:33 2019 +0200
@@ -15,18 +15,14 @@
 
 import (
 	"bufio"
-	"container/list"
 	"fmt"
 	"io"
 	"log"
 	"math"
-	"math/bits"
 	"math/rand"
 	"os"
-	"runtime"
 	"strconv"
 	"strings"
-	"sync"
 	"time"
 
 	svg "github.com/ajstarks/svgo"
@@ -91,398 +87,6 @@
 	}
 }
 
-type indexedArc struct {
-	arc   int32
-	index int32
-}
-
-func glue(a, b octree.LineStringZ) octree.LineStringZ {
-
-	a1, a2 := a[0], a[len(a)-1]
-	b1, b2 := b[0], b[len(b)-1]
-
-	if a1.EpsEquals2D(b2) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, b)
-		copy(c[len(b):], a[1:])
-		return c
-	}
-
-	if a2.EpsEquals2D(b1) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, a)
-		copy(c[len(a):], b[1:])
-		return c
-	}
-
-	if a1.EpsEquals2D(b1) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, b)
-		c[:len(b)].Reverse()
-		copy(c[len(b):], a[1:])
-		return c
-	}
-
-	if a2.EpsEquals2D(b2) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, a)
-		copy(c[len(a):], b[:len(b)-1])
-		c[len(a):].Reverse()
-		return c
-	}
-
-	return nil
-}
-
-func connectArcs(tri *octree.Triangulation, cuts []indexedArc, arcs *[]octree.LineStringZ) {
-
-	unique := map[int32]struct{}{}
-	for i := range cuts {
-		unique[cuts[i].arc] = struct{}{}
-	}
-	before := len(unique)
-
-	origLen := int32(len(*arcs))
-
-	for i := range cuts {
-		if cuts[i].arc >= origLen {
-			// already has a connected arc assigned.
-			continue
-		}
-
-		line := (*arcs)[cuts[i].arc]
-
-		var modified []int
-		visited := map[int32]struct{}{}
-
-		var recursive func(int32)
-		recursive = func(idx int32) {
-			visited[idx] = struct{}{}
-
-			ns := neighbors(tri, idx)
-			for _, n := range ns {
-				n /= 3
-				if _, already := visited[n]; already {
-					continue
-				}
-
-				arcIdx := findArc(n, cuts)
-				if arcIdx == -1 {
-					continue
-				}
-				arc := cuts[arcIdx].arc
-				if arc >= origLen {
-					// already a new arc.
-					continue
-				}
-				oline := (*arcs)[arc]
-
-				nline := glue(line, oline)
-				if nline == nil {
-					// not joinable
-					continue
-				}
-				line = nline
-				modified = append(modified, arcIdx)
-				recursive(cuts[arcIdx].index)
-			}
-		}
-		recursive(cuts[i].index)
-		if len(modified) > 0 {
-			// alloc a new line an fix the references.
-			nidx := int32(len(*arcs))
-			*arcs = append(*arcs, line)
-			cuts[i].arc = nidx
-			for _, j := range modified {
-				cuts[j].arc = nidx
-			}
-		}
-	}
-
-	unique = map[int32]struct{}{}
-	for i := range cuts {
-		unique[cuts[i].arc] = struct{}{}
-	}
-	log.Printf("unique arcs: before: %d after %d\n",
-		before, len(unique))
-}
-
-func intersectTriangles(
-	tri *octree.Triangulation,
-	heights []float64,
-	min, max octree.Vertex,
-) [][]octree.LineStringZ {
-
-	cuts := make([][]indexedArc, len(heights))
-
-	classes := make([][]int32, len(heights)+1)
-
-	var arcs []octree.LineStringZ
-
-all:
-	for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
-		ti := tri.Triangles[i : i+3]
-		v0 := tri.Points[ti[0]]
-		v1 := tri.Points[ti[1]]
-		v2 := tri.Points[ti[2]]
-		min := math.Min(v0.Z, math.Min(v1.Z, v2.Z))
-
-		if min > heights[len(heights)-1] {
-			classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
-			continue
-		}
-		max := math.Max(v0.Z, math.Max(v1.Z, v2.Z))
-
-		for j, h := range heights {
-
-			var l float64
-			if j > 0 {
-				l = heights[j-1]
-			} else {
-				l = -math.MaxFloat64
-			}
-
-			if l < min && max < h {
-				classes[j] = append(classes[j], i/3)
-				continue all
-			}
-			if min > h || max < h {
-				continue
-			}
-			t := octree.Triangle{v0, v1, v2}
-			cut := t.IntersectHorizontal(h)
-			if len(cut) >= 2 {
-				arc := int32(len(arcs))
-				arcs = append(arcs, cut)
-				cuts[j] = append(cuts[j], indexedArc{arc: arc, index: i / 3})
-			}
-		}
-	}
-
-	// connect the arcs in a cut list to longer arcs.
-
-	for _, c := range cuts {
-		connectArcs(tri, c, &arcs)
-	}
-
-	func() {
-		out, err := os.Create("arcs.svg")
-		if err != nil {
-			log.Printf("err: %v\n", err)
-			return
-		}
-		defer func() {
-			out.Close()
-			log.Println("writing arcs done")
-		}()
-
-		buf := bufio.NewWriter(out)
-
-		dumpArcsSVG(
-			buf,
-			min, max,
-			cuts,
-			arcs)
-		buf.Flush()
-	}()
-
-	result := make([][]octree.LineStringZ, len(classes))
-
-	jobs := make(chan int)
-
-	var wg sync.WaitGroup
-
-	worker := func() {
-		defer wg.Done()
-
-		pb := polygonBuilder{open: list.New()}
-
-		for i := range jobs {
-			usedArcs := map[int32]struct{}{}
-			var dupes int
-			var isolated, inside, found int
-			c := classes[i]
-
-		allInClass:
-			for _, idx := range c {
-				ns := neighbors(tri, idx)
-				mask := where(ns, c)
-				switch bits.OnesCount8(mask) {
-				case 3:
-					// Totally insides do not contribute to the geometries.
-					inside++
-					continue allInClass
-				case 0:
-					// Isolated are areas w/o connections to the rest.
-					isolated++
-					ti := tri.Triangles[idx*3 : idx*3+3]
-					pb.addTriangle(
-						tri.Points[ti[0]],
-						tri.Points[ti[1]],
-						tri.Points[ti[2]])
-					continue allInClass
-				default:
-					ti := tri.Triangles[idx*3 : idx*3+3]
-					for j := 0; j < 3; j++ {
-						if (mask & (1 << j)) == 0 {
-
-							var curr octree.LineStringZ
-
-							for l := i - 1; l <= i; l++ {
-								if l < 0 || l >= len(cuts) {
-									continue
-								}
-								arcIdx := findArc(ns[j]/3, cuts[l])
-								if arcIdx == -1 {
-									continue
-								}
-								found++
-								aIdx := cuts[l][arcIdx].arc
-								if _, already := usedArcs[aIdx]; already {
-									dupes++
-									continue
-								}
-								usedArcs[aIdx] = struct{}{}
-
-								curr = arcs[aIdx]
-								break
-							}
-
-							if curr == nil {
-								k := (j + 1) % 3
-								front := tri.Points[ti[j]]
-								back := tri.Points[ti[k]]
-								curr = octree.LineStringZ{front, back}
-							}
-
-						segment:
-							for e := pb.open.Front(); e != nil; e = e.Next() {
-								line := e.Value.(octree.LineStringZ)
-								if nline := glue(curr, line); nline != nil {
-									curr = nline
-									pb.open.Remove(e)
-									goto segment
-								}
-							} // all open
-
-							// check if closed
-							if curr[0].EpsEquals2D(curr[len(curr)-1]) {
-								if !curr.CCW() {
-									curr.Reverse()
-								}
-								pb.polygons = append(pb.polygons, curr)
-							} else {
-								pb.open.PushFront(curr)
-							}
-						} // for all border parts
-					}
-				}
-			}
-
-			log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n",
-				i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found)
-			/*
-				for e := pb.open.Front(); e != nil; e = e.Next() {
-					line := e.Value.(octree.LineStringZ)
-					if !line.CCW() {
-						line.Reverse()
-					}
-					pb.polygons = append(pb.polygons, line)
-				}
-			*/
-
-			result[i] = pb.polygons
-			pb.reset()
-		}
-	}
-
-	for n := runtime.NumCPU(); n >= 0; n-- {
-		wg.Add(1)
-		go worker()
-	}
-
-	log.Println("inside classes:")
-	for i := range classes {
-		jobs <- i
-	}
-	close(jobs)
-
-	wg.Wait()
-
-	log.Println("cuts:")
-	for i, c := range cuts {
-		log.Printf("\t%.3f: %d\n", heights[i], len(c))
-	}
-
-	return result
-
-	// TODO: sew segments together.
-
-}
-
-type polygonBuilder struct {
-	polygons []octree.LineStringZ
-
-	open *list.List
-}
-
-func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
-	polygon := octree.LineStringZ{v0, v1, v2, v0}
-	pb.polygons = append(pb.polygons, polygon)
-}
-
-func (pb *polygonBuilder) reset() {
-	pb.polygons = pb.polygons[:0]
-	pb.open.Init()
-}
-
-func neighbors(t *octree.Triangulation, idx int32) []int32 {
-	idx *= 3
-	return t.Halfedges[idx : idx+3]
-}
-
-func findArc(needle int32, haystack []indexedArc) int {
-	lo, hi := 0, len(haystack)-1
-	for lo <= hi {
-		mid := (hi-lo)/2 + lo
-		switch v := haystack[mid].index; {
-		case v < needle:
-			lo = mid + 1
-		case v > needle:
-			hi = mid - 1
-		default:
-			return mid
-		}
-	}
-	return -1
-}
-
-func contains(needle int32, haystack []int32) bool {
-	lo, hi := 0, len(haystack)-1
-	for lo <= hi {
-		mid := (hi-lo)/2 + lo
-		switch v := haystack[mid]; {
-		case v < needle:
-			lo = mid + 1
-		case v > needle:
-			hi = mid - 1
-		default:
-			return true
-		}
-	}
-	return false
-}
-
-func where(neighbors, indices []int32) byte {
-	var mask byte
-	for i, n := range neighbors {
-		if n >= 0 && contains(n/3, indices) {
-			mask |= 1 << i
-		}
-	}
-	return mask
-}
-
 func linear(x1, y1, x2, y2 float64) func(float64) float64 {
 	// f(x1) = y1
 	// f(x2) = y2