--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/mesh/triangulation.go	Tue Nov 05 14:31:22 2019 +0100
@@ -0,0 +1,318 @@
+// Copyright (C) 2018 Michael Fogleman
+//
+// Permission is hereby granted, free of charge, to any person obtaining
+// a copy of this software and associated documentation files (the "Software"),
+// to deal in the Software without restriction, including without limitation
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,
+// and/or sell copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+package mesh
+
+import (
+	"fmt"
+	"log"
+	"math"
+
+	"gonum.org/v1/gonum/stat"
+)
+
+type Triangulation struct {
+	Points     []Vertex
+	ConvexHull []Vertex
+	Triangles  []int32
+	Halfedges  []int32
+}
+
+// Triangulate returns a Delaunay triangulation of the provided points.
+func Triangulate(points []Vertex) (*Triangulation, error) {
+	t := newTriangulator(points)
+	err := t.triangulate()
+	return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err
+}
+
+func (t *Triangulation) EstimateTooLong() float64 {
+
+	num := len(t.Triangles) / 3
+
+	lengths := make([]float64, 0, num)
+
+	points := t.Points
+
+tris:
+	for i := 0; i < num; i++ {
+		idx := i * 3
+		var max float64
+		vs := t.Triangles[idx : idx+3]
+		for j, vj := range vs {
+			if t.Halfedges[idx+j] < 0 {
+				continue tris
+			}
+			k := (j + 1) % 3
+			if l := points[vj].Distance2D(points[vs[k]]); l > max {
+				max = l
+			}
+		}
+		lengths = append(lengths, max)
+	}
+
+	std := stat.StdDev(lengths, nil)
+	return 3.5 * std
+}
+
+func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) {
+
+	if tooLong <= 0 {
+		tooLong = t.EstimateTooLong()
+	}
+
+	tooLong *= tooLong
+
+	var candidates []int32
+
+	points := t.Points
+
+	for i, num := 0, len(t.Triangles)/3; i < num; i++ {
+		idx := i * 3
+		var max float64
+		vs := t.Triangles[idx : idx+3]
+		for j, vj := range vs {
+			k := (j + 1) % 3
+			if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max {
+				max = l
+			}
+		}
+		if max > tooLong {
+			candidates = append(candidates, int32(i))
+		}
+	}
+
+	removed := map[int32]struct{}{}
+
+	isBorder := func(n int32) bool {
+		n *= 3
+		for i := int32(0); i < 3; i++ {
+			e := n + i
+			o := t.Halfedges[e]
+			if o < 0 {
+				return true
+			}
+			if _, found := removed[o/3]; found {
+				return true
+			}
+		}
+		return false
+	}
+
+	var newCandidates []int32
+
+	log.Printf("info: candidates: %d\n", len(candidates))
+	for len(candidates) > 0 {
+
+		oldRemoved := len(removed)
+
+		for _, i := range candidates {
+
+			if isBorder(i) {
+				removed[i] = struct{}{}
+			} else {
+				newCandidates = append(newCandidates, i)
+			}
+		}
+
+		if oldRemoved == len(removed) {
+			break
+		}
+
+		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))
+
+	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{
+		points[curr.a],
+		points[curr.b],
+	}
+
+	for curr = curr.next; curr != rings[0]; curr = curr.next {
+		polygon = append(polygon, 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 {
+	sl := make([][]int32, len(t.Triangles)/3)
+	var j int
+	for i := range sl {
+		sl[i] = t.Triangles[j : j+3]
+		j += 3
+	}
+	return sl
+}
+
+func (t *Triangulation) Tin() *Tin {
+
+	min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64}
+	max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64}
+
+	vertices := t.Points
+
+	for _, v := range vertices {
+		min.Minimize(v)
+		max.Maximize(v)
+	}
+
+	return &Tin{
+		Vertices:  vertices,
+		Triangles: t.TriangleSlices(),
+		Min:       min,
+		Max:       max,
+	}
+}
+
+func (t *Triangulation) area() float64 {
+	var result float64
+	points := t.Points
+	ts := t.Triangles
+	for i := 0; i < len(ts); i += 3 {
+		p0 := points[ts[i+0]]
+		p1 := points[ts[i+1]]
+		p2 := points[ts[i+2]]
+		result += area(p0, p1, p2)
+	}
+	return result / 2
+}
+
+// Validate performs several sanity checks on the Triangulation to check for
+// potential errors. Returns nil if no issues were found. You normally
+// shouldn't need to call this function but it can be useful for debugging.
+func (t *Triangulation) Validate() error {
+	// verify halfedges
+	for i1, i2 := range t.Halfedges {
+		if i1 != -1 && t.Halfedges[i1] != i2 {
+			return fmt.Errorf("invalid halfedge connection")
+		}
+		if i2 != -1 && t.Halfedges[i2] != int32(i1) {
+			return fmt.Errorf("invalid halfedge connection")
+		}
+	}
+
+	// verify convex hull area vs sum of triangle areas
+	hull1 := t.ConvexHull
+	hull2 := ConvexHull(t.Points)
+	area1 := polygonArea(hull1)
+	area2 := polygonArea(hull2)
+	area3 := t.area()
+	if math.Abs(area1-area2) > 1e-9 || math.Abs(area1-area3) > 1e-9 {
+		return fmt.Errorf("hull areas disagree: %f, %f, %f", area1, area2, area3)
+	}
+
+	// verify convex hull perimeter
+	perimeter1 := polygonPerimeter(hull1)
+	perimeter2 := polygonPerimeter(hull2)
+	if math.Abs(perimeter1-perimeter2) > 1e-9 {
+		return fmt.Errorf("hull perimeters disagree: %f, %f", perimeter1, perimeter2)
+	}
+
+	return nil
+}