--- a/pkg/imports/sr.go	Wed Jun 05 09:38:59 2019 +0200
+++ b/pkg/imports/sr.go	Wed Jun 05 11:39:25 2019 +0200
@@ -28,13 +28,11 @@
 	"os"
 	"path"
 	"path/filepath"
-	"sort"
 	"strconv"
 	"strings"
 	"time"
 
 	shp "github.com/jonas-p/go-shp"
-	"gonum.org/v1/gonum/stat"
 
 	"gemma.intevation.de/gemma/pkg/common"
 	"gemma.intevation.de/gemma/pkg/misc"
@@ -366,36 +364,7 @@
 		feedback.Info("No boundary given.")
 		feedback.Info("Eliminate triangles with long edges.")
 
-		lens := make([]float64, len(tri.Triangles)/3)
-		for i := 0; i < len(tri.Triangles); i += 3 {
-			vs := tri.Triangles[i : i+3]
-			var max float64
-			for j := 0; j < 3; j++ {
-				k := (j + 1) % 3
-				if l := tri.Points[vs[j]].SquaredDistance2D(tri.Points[vs[k]]); l > max {
-					max = l
-				}
-			}
-			lens[i/3] = max
-		}
-		qs := make([]float64, len(lens))
-		copy(qs, lens)
-		sort.Float64s(qs)
-		q := stat.Quantile(0.98, stat.LinInterp, qs, nil)
-		feedback.Info("98%% quantile of the edge lengths: %.2f\n", math.Sqrt(q))
-
-		ts := tri.Triangles
-
-		for i := len(lens) - 1; i >= 0; i-- {
-			if lens[i] >= q {
-				if i != len(lens)-1 {
-					copy(ts[i*3:], ts[(i+1)*3:])
-				}
-				ts = ts[:len(ts)-3]
-			}
-		}
-		feedback.Info("Remaining triangles: %d (%.2f%%)\n",
-			len(ts)/3, float64(len(ts))/float64(len(tri.Triangles))*100.0)
+		tri.ConcaveHull(0.98)
 	}
 
 	// TODO: Implement me!

--- 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