--- a/pkg/octree/triangulation.go	Wed Jun 05 09:38:59 2019 +0200
+++ b/pkg/octree/triangulation.go	Wed Jun 05 11:39:25 2019 +0200
@@ -21,7 +21,11 @@
 
 import (
 	"fmt"
+	"log"
 	"math"
+	"sort"
+
+	"gonum.org/v1/gonum/stat"
 )
 
 type Triangulation struct {
@@ -38,6 +42,82 @@
 	return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err
 }
 
+func (t *Triangulation) ConcaveHull(quantile float64) []int32 {
+
+	numTris := len(t.Triangles) / 3
+
+	// Calculate the max edge lengths per triangle.
+	lengths := make([]float64, numTris)
+	for i := range lengths {
+		idx := i * 3
+		var max float64
+		vs := t.Triangles[idx : idx+3]
+		for j, vj := range vs {
+			k := (j + 1) % 3
+			if l := t.Points[vj].SquaredDistance2D(t.Points[vs[k]]); l > max {
+				max = l
+			}
+		}
+		lengths[i] = max
+	}
+
+	qs := make([]float64, len(lengths))
+	copy(qs, lengths)
+	sort.Float64s(qs)
+	q := stat.Quantile(quantile, stat.LinInterp, qs, nil)
+	log.Printf("info: %.2f%% quantile: %.2f\n", quantile, math.Sqrt(q))
+
+	removed := map[int32]bool{}
+
+	isBorder := func(n int32) bool {
+		n *= 3
+		for i := int32(0); i < 3; i++ {
+			e := n + i
+			if t.Halfedges[e] == -1 || removed[e] {
+				return true
+			}
+		}
+		return false
+	}
+
+	oldCandidates := make([]int32, 0, int(math.Ceil(float64(numTris)*(1.0-quantile))))
+	log.Printf("info: cap candidates: %d\n", cap(oldCandidates))
+
+	for i, length := range lengths {
+		if length >= q {
+			oldCandidates = append(oldCandidates, int32(i))
+		}
+	}
+
+	var newCandidates []int32
+
+	for len(oldCandidates) > 0 {
+		log.Printf("candidates: %d\n", len(oldCandidates))
+
+		oldRemoved := len(removed)
+
+		for _, i := range oldCandidates {
+
+			if isBorder(i) {
+				removed[i] = true
+			} else {
+				newCandidates = append(newCandidates, i)
+			}
+		}
+
+		if oldRemoved == len(removed) {
+			break
+		}
+
+		oldCandidates = newCandidates
+		newCandidates = newCandidates[:0]
+	}
+
+	log.Printf("info: triangles to remove: %d\n", len(removed))
+
+	return nil
+}
+
 func (t *Triangulation) TriangleSlices() [][]int32 {
 	sl := make([][]int32, len(t.Triangles)/3)
 	var j int