--- a/pkg/imports/sr.go	Thu Jun 20 19:43:25 2019 +0200
+++ b/pkg/imports/sr.go	Mon Jun 24 11:39:09 2019 +0200
@@ -68,7 +68,6 @@
 )
 
 const (
-	tooLongEdge          = 50.0
 	pointsPerSquareMeter = 2
 )
 
@@ -398,6 +397,7 @@
 
 	if boundary == nil {
 		feedback.Info("No boundary given. Calulate from XYZ data.")
+		tooLongEdge := tri.EstimateTooLong()
 		feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
 
 		var polygon octree.LineStringZ

--- a/pkg/octree/triangulation.go	Thu Jun 20 19:43:25 2019 +0200
+++ b/pkg/octree/triangulation.go	Mon Jun 24 11:39:09 2019 +0200
@@ -23,6 +23,8 @@
 	"fmt"
 	"log"
 	"math"
+
+	"gonum.org/v1/gonum/stat"
 )
 
 type Triangulation struct {
@@ -39,19 +41,54 @@
 	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 := t.Points[vj].SquaredDistance2D(t.Points[vs[k]]); l > max {
+			if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max {
 				max = l
 			}
 		}
@@ -188,12 +225,12 @@
 	curr := rings[0]
 
 	polygon := LineStringZ{
-		t.Points[curr.a],
-		t.Points[curr.b],
+		points[curr.a],
+		points[curr.b],
 	}
 
 	for curr = curr.next; curr != rings[0]; curr = curr.next {
-		polygon = append(polygon, t.Points[curr.b])
+		polygon = append(polygon, points[curr.b])
 	}
 
 	polygon = append(polygon, t.Points[rings[0].a])