# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1560337034 -7200
# Node ID 10471aa73ad8cc24b9f57bc7b64f34fbbf9c8e41
# Parent  2a4216c81e7bc6fa42e432133f6c974d74ecc3cd# Parent  89a39783c20a0798811baccf5940358fac1a667c
Merged default into single-beam branch.

diff -r 89a39783c20a -r 10471aa73ad8 3rdpartylibs.sh
--- a/3rdpartylibs.sh	Tue Jun 11 20:32:25 2019 +0200
+++ b/3rdpartylibs.sh	Wed Jun 12 12:57:14 2019 +0200
@@ -44,6 +44,10 @@
 go get -u -v gonum.org/v1/gonum/stat
 # BSD-3-Clause
 
+# Only needed when generating SVG graphics for debugging.
+# go get -u -v github.com/ajstarks/svgo
+# Attribution 3.0 United States (CC BY 3.0 US)
+
 ## list of additional licenses that get fetched and installed as dependencies
 # github.com/fsnotify/fsnotify/ BSD-3-Clause
 # github.com/hashicorp/hcl/ MPL-2.0
diff -r 89a39783c20a -r 10471aa73ad8 pkg/octree/triangulation.go
--- a/pkg/octree/triangulation.go	Tue Jun 11 20:32:25 2019 +0200
+++ b/pkg/octree/triangulation.go	Wed Jun 12 12:57:14 2019 +0200
@@ -39,11 +39,11 @@
 	return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err
 }
 
-func (t *Triangulation) ConcaveHull(tooLong float64) []int32 {
+func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) {
 
 	tooLong *= tooLong
 
-	var oldCandidates []int32
+	var candidates []int32
 
 	for i, num := 0, len(t.Triangles)/3; i < num; i++ {
 		idx := i * 3
@@ -56,17 +56,21 @@
 			}
 		}
 		if max > tooLong {
-			oldCandidates = append(oldCandidates, int32(i))
+			candidates = append(candidates, int32(i))
 		}
 	}
 
-	removed := map[int32]bool{}
+	removed := map[int32]struct{}{}
 
 	isBorder := func(n int32) bool {
 		n *= 3
 		for i := int32(0); i < 3; i++ {
 			e := n + i
-			if t.Halfedges[e] == -1 || removed[e] {
+			o := t.Halfedges[e]
+			if o < 0 {
+				return true
+			}
+			if _, found := removed[o/3]; found {
 				return true
 			}
 		}
@@ -75,15 +79,15 @@
 
 	var newCandidates []int32
 
-	for len(oldCandidates) > 0 {
-		log.Printf("candidates: %d\n", len(oldCandidates))
+	log.Printf("info: candidates: %d\n", len(candidates))
+	for len(candidates) > 0 {
 
 		oldRemoved := len(removed)
 
-		for _, i := range oldCandidates {
+		for _, i := range candidates {
 
 			if isBorder(i) {
-				removed[i] = true
+				removed[i] = struct{}{}
 			} else {
 				newCandidates = append(newCandidates, i)
 			}
@@ -93,13 +97,110 @@
 			break
 		}
 
-		oldCandidates = newCandidates
+		candidates = newCandidates
 		newCandidates = newCandidates[:0]
 	}
 
+	log.Printf("info: candidates left: %d\n", len(candidates))
+	log.Printf("info: triangles: %d\n", len(t.Triangles)/3)
 	log.Printf("info: triangles to remove: %d\n", len(removed))
 
-	return nil
+	type edge struct {
+		a, b       int32
+		prev, next *edge
+	}
+
+	isClosed := func(e *edge) bool {
+		for curr := e.next; curr != nil; curr = curr.next {
+			if curr == e {
+				return true
+			}
+		}
+		return false
+	}
+
+	open := map[int32]*edge{}
+	var rings []*edge
+
+	for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ {
+		if _, found := removed[i]; found {
+			continue
+		}
+		n := i * 3
+		for j := int32(0); j < 3; j++ {
+			e := n + j
+			f := t.Halfedges[e]
+			if f >= 0 {
+				if _, found := removed[f/3]; !found {
+					continue
+				}
+			}
+			a := t.Triangles[e]
+			b := t.Triangles[n+(j+1)%3]
+
+			curr := &edge{a: a, b: b}
+
+			if old := open[a]; old != nil {
+				delete(open, a)
+				if old.a == a {
+					old.prev = curr
+					curr.next = old
+				} else {
+					old.next = curr
+					curr.prev = old
+				}
+
+				if isClosed(old) {
+					rings = append(rings, old)
+				}
+			} else {
+				open[a] = curr
+			}
+
+			if old := open[b]; old != nil {
+				delete(open, b)
+				if old.b == b {
+					old.next = curr
+					curr.prev = old
+				} else {
+					old.prev = curr
+					curr.next = old
+				}
+
+				if isClosed(old) {
+					rings = append(rings, old)
+				}
+			} else {
+				open[b] = curr
+			}
+		}
+	}
+
+	if len(open) > 0 {
+		log.Printf("warn: open vertices left: %d\n", len(open))
+	}
+
+	if len(rings) == 0 {
+		log.Println("warn: no ring found")
+		return nil, removed
+	}
+
+	curr := rings[0]
+
+	polygon := LineStringZ{
+		t.Points[curr.a],
+		t.Points[curr.b],
+	}
+
+	for curr = curr.next; curr != rings[0]; curr = curr.next {
+		polygon = append(polygon, t.Points[curr.b])
+	}
+
+	polygon = append(polygon, t.Points[rings[0].a])
+
+	log.Printf("length of boundary: %d\n", len(polygon))
+
+	return polygon, removed
 }
 
 func (t *Triangulation) TriangleSlices() [][]int32 {
diff -r 89a39783c20a -r 10471aa73ad8 pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Tue Jun 11 20:32:25 2019 +0200
+++ b/pkg/octree/vertex.go	Wed Jun 12 12:57:14 2019 +0200
@@ -413,6 +413,41 @@
 	}
 }
 
+func (ls LineStringZ) BBox() Box2D {
+
+	min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64}
+	max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64}
+
+	for _, v := range ls {
+		min.Minimize(v)
+		max.Maximize(v)
+	}
+
+	return Box2D{
+		X1: min.X,
+		Y1: min.Y,
+		X2: max.X,
+		Y2: max.Y,
+	}
+}
+
+func (ls LineStringZ) Area() float64 {
+	var result float64
+	for i := 0; i < len(ls); i += 3 {
+		p0 := ls[i+0]
+		p1 := ls[i+1]
+		p2 := ls[i+2]
+		result += area(p0, p1, p2)
+	}
+	return result / 2
+}
+
+func (ls LineStringZ) Reverse() {
+	for i, j := 0, len(ls)-1; i < j; i, j = i+1, j-1 {
+		ls[i], ls[j] = ls[j], ls[i]
+	}
+}
+
 func (ls LineStringZ) order(position func(Vertex) float64) {
 	type posVertex struct {
 		pos float64
@@ -738,6 +773,10 @@
 	}
 }
 
+func (a Box2D) Area() float64 {
+	return (a.X2 - a.X1) * (a.Y2 - a.Y1)
+}
+
 // NewPlane2D creates a new Plane2D from two given points.
 func NewPlane2D(x1, y1, x2, y2 float64) Plane2D {
 	b := x2 - x1