# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1560353420 -7200
# Node ID 810b28f59b8baacd61ea8f73bcfea0bcc7d0de8b
# Parent  10471aa73ad8cc24b9f57bc7b64f34fbbf9c8e41
Generate random points for second mesh.

diff -r 10471aa73ad8 -r 810b28f59b8b pkg/common/random.go
--- a/pkg/common/random.go	Wed Jun 12 12:57:14 2019 +0200
+++ b/pkg/common/random.go	Wed Jun 12 17:30:20 2019 +0200
@@ -19,6 +19,8 @@
 	"io"
 	"log"
 	"math/big"
+	mrand "math/rand"
+	"time"
 )
 
 // GenerateRandomKey generates a cryptographically secure random key
@@ -64,3 +66,12 @@
 	out[0] = special[0]
 	return string(out)
 }
+
+func Random(low, high float64) func() float64 {
+	if low > high {
+		low, high = high, low
+	}
+	rnd := mrand.New(mrand.NewSource(time.Now().Unix()))
+	m := high - low
+	return func() float64 { return rnd.Float64()*m + low }
+}
diff -r 10471aa73ad8 -r 810b28f59b8b pkg/imports/sr.go
--- 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
 	}
diff -r 10471aa73ad8 -r 810b28f59b8b pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Wed Jun 12 12:57:14 2019 +0200
+++ b/pkg/octree/vertex.go	Wed Jun 12 17:30:20 2019 +0200
@@ -432,14 +432,7 @@
 }
 
 func (ls LineStringZ) Area() float64 {
-	var result float64
-	for i := 0; i < len(ls); i += 3 {
-		p0 := ls[i+0]
-		p1 := ls[i+1]
-		p2 := ls[i+2]
-		result += area(p0, p1, p2)
-	}
-	return result / 2
+	return polygonArea(ls)
 }
 
 func (ls LineStringZ) Reverse() {