# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1560336978 -7200
# Node ID 2a4216c81e7bc6fa42e432133f6c974d74ecc3cd
# Parent  db7136854a51b117fdf6e07eccba79fee4d00b28
Extract the removed triangles from first triangulation, too. Useful to build a artifical DEM for second pass.

diff -r db7136854a51 -r 2a4216c81e7b pkg/octree/triangulation.go
--- a/pkg/octree/triangulation.go	Tue Jun 11 18:20:41 2019 +0200
+++ b/pkg/octree/triangulation.go	Wed Jun 12 12:56:18 2019 +0200
@@ -39,7 +39,7 @@
 	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
 
@@ -60,13 +60,17 @@
 		}
 	}
 
-	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 o := t.Halfedges[e]; o < 0 || removed[o/3] {
+			o := t.Halfedges[e]
+			if o < 0 {
+				return true
+			}
+			if _, found := removed[o/3]; found {
 				return true
 			}
 		}
@@ -83,7 +87,7 @@
 		for _, i := range candidates {
 
 			if isBorder(i) {
-				removed[i] = true
+				removed[i] = struct{}{}
 			} else {
 				newCandidates = append(newCandidates, i)
 			}
@@ -119,51 +123,55 @@
 	var rings []*edge
 
 	for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ {
-		if removed[i] {
+		if _, found := removed[i]; found {
 			continue
 		}
 		n := i * 3
 		for j := int32(0); j < 3; j++ {
 			e := n + j
-			if f := t.Halfedges[e]; f < 0 || removed[f/3] {
-				a := t.Triangles[e]
-				b := t.Triangles[n+(j+1)%3]
-
-				curr := &edge{a: a, b: b}
+			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]
 
-				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
-					}
+			curr := &edge{a: a, b: b}
 
-					if isClosed(old) {
-						rings = append(rings, old)
-					}
+			if old := open[a]; old != nil {
+				delete(open, a)
+				if old.a == a {
+					old.prev = curr
+					curr.next = old
 				} else {
-					open[a] = curr
+					old.next = curr
+					curr.prev = old
 				}
 
-				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[a] = curr
+			}
 
-					if isClosed(old) {
-						rings = append(rings, old)
-					}
+			if old := open[b]; old != nil {
+				delete(open, b)
+				if old.b == b {
+					old.next = curr
+					curr.prev = old
 				} else {
-					open[b] = curr
+					old.prev = curr
+					curr.next = old
 				}
+
+				if isClosed(old) {
+					rings = append(rings, old)
+				}
+			} else {
+				open[b] = curr
 			}
 		}
 	}
@@ -174,19 +182,25 @@
 
 	if len(rings) == 0 {
 		log.Println("warn: no ring found")
-		return nil
+		return nil, removed
 	}
 
 	curr := rings[0]
-	vertices := []int32{curr.a, curr.b}
+
+	polygon := LineStringZ{
+		t.Points[curr.a],
+		t.Points[curr.b],
+	}
 
 	for curr = curr.next; curr != rings[0]; curr = curr.next {
-		vertices = append(vertices, curr.b)
+		polygon = append(polygon, t.Points[curr.b])
 	}
 
-	log.Printf("length of boundary: %d\n", len(vertices))
+	polygon = append(polygon, t.Points[rings[0].a])
 
-	return vertices
+	log.Printf("length of boundary: %d\n", len(polygon))
+
+	return polygon, removed
 }
 
 func (t *Triangulation) TriangleSlices() [][]int32 {
diff -r db7136854a51 -r 2a4216c81e7b pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Tue Jun 11 18:20:41 2019 +0200
+++ b/pkg/octree/vertex.go	Wed Jun 12 12:56:18 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