--- a/pkg/imports/sr.go	Wed Jun 12 12:57:14 2019 +0200
+++ b/pkg/imports/sr.go	Wed Jun 12 17:30:20 2019 +0200
@@ -24,12 +24,15 @@
 	"errors"
 	"fmt"
 	"io"
+	"log"
 	"math"
 	"os"
 	"path"
 	"path/filepath"
+	"runtime"
 	"strconv"
 	"strings"
+	"sync"
 	"time"
 
 	shp "github.com/jonas-p/go-shp"
@@ -67,6 +70,11 @@
 	contourTolerance = 0.1
 )
 
+const (
+	tooLongEdge          = 50.0
+	pointsPerSquareMeter = 5
+)
+
 // SRJobKind is the unique name of a SoundingResult import job.
 const SRJobKind JobKind = "sr"
 
@@ -360,11 +368,122 @@
 	feedback.Info("Triangulation took %v.", time.Since(start))
 	feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
 
+	var clippingPolygon octree.Polygon
+
 	if boundary == nil {
 		feedback.Info("No boundary given.")
-		feedback.Info("Eliminate triangles with long edges.")
+
+		polygon, removed := tri.ConcaveHull(tooLongEdge)
+
+		polygonArea := polygon.Area()
+		if polygonArea < 0.0 { // counter clockwise
+			polygonArea = -polygonArea
+			polygon.Reverse()
+		}
+
+		firstTin := tri.Tin()
+		builder := octree.NewBuilder(firstTin)
+		builder.Build(removed)
+
+		firstTree := builder.Tree()
+
+		bbox := polygon.BBox()
+		bboxArea := bbox.Area()
+
+		log.Printf("bbox area: %.2f\n", bboxArea)
+		log.Printf("polygon area: %.2f\n", polygonArea)
+
+		numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
+
+		log.Printf("random points for density of ~%d points/m²: %d\n",
+			pointsPerSquareMeter, numPoints)
+
+		start := time.Now()
+
+		var wg sync.WaitGroup
+
+		jobs := make(chan int)
+		out := make(chan []octree.Vertex)
+		done := make(chan struct{})
+
+		cpus := runtime.NumCPU()
+
+		free := make(chan []octree.Vertex, cpus)
+
+		for i, n := 0, cpus; i < n; i++ {
+			wg.Add(1)
+			go func() {
+				defer wg.Done()
+				xRange := common.Random(bbox.X1, bbox.X2)
+				yRange := common.Random(bbox.Y1, bbox.Y2)
 
-		tri.ConcaveHull(50.0)
+				for size := range jobs {
+					var vertices []octree.Vertex
+					select {
+					case vertices = <-free:
+					default:
+						vertices = make([]octree.Vertex, 0, 1000)
+					}
+					for len(vertices) < size {
+						x, y := xRange(), yRange()
+						if z, ok := firstTree.Value(x, y); ok {
+							vertices = append(vertices, octree.Vertex{X: x, Y: y, Z: z})
+						}
+					}
+					out <- vertices
+				}
+			}()
+		}
+
+		generated := make(octree.LineStringZ, 0, numPoints+len(polygon)-1)
+
+		go func() {
+			defer close(done)
+			for vertices := range out {
+				generated = append(generated, vertices...)
+				select {
+				case free <- vertices[:0]:
+				default:
+				}
+			}
+		}()
+
+		for remain := numPoints; remain > 0; {
+			if remain > 1000 {
+				jobs <- 1000
+				remain -= 1000
+			} else {
+				jobs <- remain
+				remain = 0
+			}
+		}
+		close(jobs)
+		wg.Wait()
+		close(out)
+		<-done
+
+		log.Printf("Generating %d points took %v\n", len(generated), time.Since(start))
+
+		// Add the boundary to new point cloud.
+		generated = append(generated, polygon[:len(polygon)-1]...)
+
+		xyz = octree.MultiPointZ(generated)
+
+	} else {
+		var hull []byte
+		if err = tx.QueryRowContext(
+			ctx,
+			reprojectPointsSQL,
+			boundary.asWKB(),
+			m.EPSG,
+			epsg,
+		).Scan(&hull); err != nil {
+			return nil, err
+		}
+		if err := clippingPolygon.FromWKB(hull); err != nil {
+			return nil, err
+		}
+		clippingPolygon.Indexify()
 	}
 
 	// TODO: Implement me!
@@ -419,6 +538,7 @@
 	start = time.Now()
 
 	var clippingPolygon octree.Polygon
+
 	if err := clippingPolygon.FromWKB(hull); err != nil {
 		return nil, err
 	}