# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1560506627 -7200
# Node ID b30c0b02e8fbcfdeeb8dfff067649bf838df1496
# Parent  7126fdc5779a02801fef9e8494b02af4667df6d5# Parent  66f2cb7899058e9dad43addda79f8f498292ba33
Merged single-beam branch into default.

diff -r 7126fdc5779a -r b30c0b02e8fb 3rdpartylibs.sh
--- a/3rdpartylibs.sh	Fri Jun 14 11:16:09 2019 +0200
+++ b/3rdpartylibs.sh	Fri Jun 14 12:03:47 2019 +0200
@@ -44,6 +44,10 @@
 go get -u -v gonum.org/v1/gonum/stat
 # BSD-3-Clause
 
+# Only needed when generating SVG graphics for debugging.
+# go get -u -v github.com/ajstarks/svgo
+# Attribution 3.0 United States (CC BY 3.0 US)
+
 ## list of additional licenses that get fetched and installed as dependencies
 # github.com/fsnotify/fsnotify/ BSD-3-Clause
 # github.com/hashicorp/hcl/ MPL-2.0
diff -r 7126fdc5779a -r b30c0b02e8fb pkg/common/random.go
--- a/pkg/common/random.go	Fri Jun 14 11:16:09 2019 +0200
+++ b/pkg/common/random.go	Fri Jun 14 12:03:47 2019 +0200
@@ -18,7 +18,10 @@
 	"crypto/rand"
 	"io"
 	"log"
+	"math"
 	"math/big"
+	mrand "math/rand"
+	"time"
 )
 
 // GenerateRandomKey generates a cryptographically secure random key
@@ -64,3 +67,20 @@
 	out[0] = special[0]
 	return string(out)
 }
+
+func Random(low, high float64) func() float64 {
+	if low > high {
+		low, high = high, low
+	}
+
+	var seed int64
+	if seedInt, err := rand.Int(rand.Reader, big.NewInt(math.MaxInt64)); err != nil {
+		log.Printf("warn: Generating good random seed failed: %v\n", err)
+		seed = time.Now().Unix()
+	} else {
+		seed = seedInt.Int64()
+	}
+	rnd := mrand.New(mrand.NewSource(seed))
+	m := high - low
+	return func() float64 { return rnd.Float64()*m + low }
+}
diff -r 7126fdc5779a -r b30c0b02e8fb pkg/imports/sr.go
--- a/pkg/imports/sr.go	Fri Jun 14 11:16:09 2019 +0200
+++ b/pkg/imports/sr.go	Fri Jun 14 12:03:47 2019 +0200
@@ -67,6 +67,11 @@
 	contourTolerance = 0.1
 )
 
+const (
+	tooLongEdge          = 50.0
+	pointsPerSquareMeter = 2
+)
+
 // SRJobKind is the unique name of a SoundingResult import job.
 const SRJobKind JobKind = "sr"
 
@@ -122,7 +127,7 @@
     CASE WHEN $5::bytea IS NULL THEN
       ST_Transform(ST_ConcaveHull(ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography
     ELSE
-      ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography
+      ST_MakeValid(ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326))::geography
     END)
 FROM waterway.bottlenecks
 WHERE objnam = $1 AND validity @> CAST($2 AS timestamptz)
@@ -135,13 +140,24 @@
 `
 
 	reprojectPointsSQL = `
-SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))
-`
+SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))`
+
+	reprojectPointsBufferedSQL = `
+SELECT
+  ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer)),
+  ST_AsBinary(ST_Buffer(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer), 0.1)),
+  ST_Area(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))`
+
 	insertOctreeSQL = `
 UPDATE waterway.sounding_results SET
   octree_checksum = $2, octree_index = $3
 WHERE id = $1`
 
+	repairBoundarySQL = `
+SELECT
+  ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.0)),
+  ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.1))`
+
 	insertContourSQL = `
 INSERT INTO waterway.sounding_results_contour_lines (
   sounding_result_id,
@@ -338,6 +354,7 @@
 ) (interface{}, error) {
 
 	feedback.Info("Processing as single beam scan.")
+	feedback.Info("Reproject XYZ data.")
 
 	start := time.Now()
 
@@ -361,6 +378,7 @@
 	feedback.Info("Reprojecting points to EPSG %d took %v.",
 		epsg, time.Since(start))
 	feedback.Info("Number of reprojected points: %d", len(xyz))
+	feedback.Info("Triangulate XYZ data.")
 
 	start = time.Now()
 	tri, err := octree.Triangulate(xyz)
@@ -370,15 +388,229 @@
 	feedback.Info("Triangulation took %v.", time.Since(start))
 	feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
 
+	var (
+		clippingPolygon         octree.Polygon
+		clippingPolygonBuffered octree.Polygon
+		removed                 map[int32]struct{}
+		polygonArea             float64
+		clippingPolygonWKB      []byte
+	)
+
 	if boundary == nil {
-		feedback.Info("No boundary given.")
-		feedback.Info("Eliminate triangles with long edges.")
+		feedback.Info("No boundary given. Calulate from XYZ data.")
+		feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
+
+		var polygon octree.LineStringZ
+		start = time.Now()
+		polygon, removed = tri.ConcaveHull(tooLongEdge)
+
+		polygonArea = polygon.Area()
+		if polygonArea < 0.0 { // counter clockwise
+			polygonArea = -polygonArea
+			polygon.Reverse()
+		}
+
+		clippingPolygon.FromLineStringZ(polygon)
+
+		var buffered []byte
+		if err := tx.QueryRowContext(
+			ctx,
+			repairBoundarySQL,
+			clippingPolygon.AsWKB(),
+			epsg,
+		).Scan(&clippingPolygonWKB, &buffered); err != nil {
+			return nil, err
+		}
+
+		if err := clippingPolygon.FromWKB(clippingPolygonWKB); err != nil {
+			return nil, err
+		}
+		if err := clippingPolygonBuffered.FromWKB(buffered); err != nil {
+			return nil, err
+		}
+
+		feedback.Info("Calculating took %v.", time.Since(start))
 
-		tri.ConcaveHull(50.0)
+	} else { // has Boundary
+		feedback.Info("Using uploaded boundary polygon.")
+		var buffered []byte
+		if err = tx.QueryRowContext(
+			ctx,
+			reprojectPointsBufferedSQL,
+			boundary.asWKB(),
+			m.EPSG,
+			epsg,
+		).Scan(&clippingPolygonWKB, &buffered, &polygonArea); err != nil {
+			return nil, err
+		}
+		if err := clippingPolygon.FromWKB(clippingPolygonWKB); err != nil {
+			return nil, err
+		}
+		if err := clippingPolygonBuffered.FromWKB(buffered); err != nil {
+			return nil, err
+		}
+
+		tin := tri.Tin()
+		tin.EPSG = uint32(epsg)
+
+		var str octree.STRTree
+		str.Build(tin)
+
+		removed = str.Clip(&clippingPolygon)
+	}
+
+	// Build the first mesh to generate random points on.
+
+	feedback.Info("Build virtual DEM based on original XYZ data.")
+
+	start = time.Now()
+
+	var tree *octree.Tree
+	{
+		builder := octree.NewBuilder(tri.Tin())
+		builder.Build(removed)
+		tree = builder.Tree()
 	}
 
-	// TODO: Implement me!
-	return nil, errors.New("Not implemented, yet!")
+	feedback.Info("Building took %v", time.Since(start))
+
+	feedback.Info("Boundary area: %.2fm²", polygonArea)
+
+	numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
+
+	feedback.Info("Generate %d random points for an average density of ~%d points/m².",
+		numPoints, pointsPerSquareMeter)
+
+	start = time.Now()
+
+	generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0))
+
+	tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
+		generated = append(generated, vertices...)
+	})
+
+	feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
+
+	// Add the boundary to new point cloud.
+	dupes := map[[2]float64]struct{}{}
+	clippingPolygon.Vertices(0, func(x, y float64) {
+		key := [2]float64{x, y}
+		if _, found := dupes[key]; found {
+			return
+		}
+		dupes[key] = struct{}{}
+		if z, ok := tree.Value(x, y); ok {
+			generated = append(generated, octree.Vertex{X: x, Y: y, Z: z})
+		}
+	})
+
+	feedback.Info("Triangulate new point cloud.")
+
+	start = time.Now()
+
+	xyz = octree.MultiPointZ(generated)
+
+	tri, err = octree.Triangulate(xyz)
+	if err != nil {
+		return nil, err
+	}
+	feedback.Info("Second triangulation took %v.", time.Since(start))
+	feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
+	feedback.Info("Clipping triangles from new mesh.")
+
+	start = time.Now()
+	tin := tri.Tin()
+	tin.EPSG = uint32(epsg)
+
+	var str octree.STRTree
+	str.Build(tin)
+	feedback.Info("Building STR tree took %v", time.Since(start))
+
+	start = time.Now()
+
+	clippingPolygonBuffered.Indexify()
+	removed = str.Clip(&clippingPolygonBuffered)
+	feedback.Info("Clipping STR tree took %v.", time.Since(start))
+	feedback.Info("Number of triangles to clip %d.", len(removed))
+
+	start = time.Now()
+
+	feedback.Info("Build final octree index")
+
+	builder := octree.NewBuilder(tin)
+	builder.Build(removed)
+	octreeIndex, err := builder.Bytes()
+	if err != nil {
+		return nil, err
+	}
+	feedback.Info("Building octree took %v.", time.Since(start))
+	feedback.Info("Store octree.")
+
+	start = time.Now()
+
+	var (
+		id       int64
+		dummy    uint32
+		lat, lon float64
+	)
+
+	var hull []byte
+
+	if err := tx.QueryRowContext(
+		ctx,
+		insertHullSQL,
+		m.Bottleneck,
+		m.Date.Time,
+		m.DepthReference,
+		nil,
+		clippingPolygonWKB,
+		epsg,
+	).Scan(
+		&id,
+		&lat,
+		&lon,
+		&dummy,
+		&hull,
+	); err != nil {
+		return nil, err
+	}
+
+	h := sha1.New()
+	h.Write(octreeIndex)
+	checksum := hex.EncodeToString(h.Sum(nil))
+	_, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
+	if err != nil {
+		return nil, err
+	}
+	feedback.Info("Storing octree index took %s.", time.Since(start))
+	feedback.Info("Generate contour lines")
+
+	start = time.Now()
+	err = generateContours(ctx, tx, builder.Tree(), id)
+	if err != nil {
+		return nil, err
+	}
+	feedback.Info("Generating contour lines took %s.",
+		time.Since(start))
+
+	// Store for potential later removal.
+	if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
+		return nil, err
+	}
+
+	summary := struct {
+		Bottleneck string      `json:"bottleneck"`
+		Date       models.Date `json:"date"`
+		Lat        float64     `json:"lat"`
+		Lon        float64     `json:"lon"`
+	}{
+		Bottleneck: m.Bottleneck,
+		Date:       m.Date,
+		Lat:        lat,
+		Lon:        lon,
+	}
+
+	return &summary, nil
 }
 
 func (sr *SoundingResult) multiBeamScan(
@@ -433,6 +665,7 @@
 	start = time.Now()
 
 	var clippingPolygon octree.Polygon
+
 	if err := clippingPolygon.FromWKB(hull); err != nil {
 		return nil, err
 	}
diff -r 7126fdc5779a -r b30c0b02e8fb pkg/octree/polygon.go
--- a/pkg/octree/polygon.go	Fri Jun 14 11:16:09 2019 +0200
+++ b/pkg/octree/polygon.go	Fri Jun 14 12:03:47 2019 +0200
@@ -409,6 +409,47 @@
 	return raySlope >= diagSlope
 }
 
+func (p *Polygon) NumVertices(ring int) int {
+	if ring < 0 || ring >= len(p.rings) {
+		return 0
+	}
+	return len(p.rings[ring]) / 2
+}
+
+func (p *Polygon) Vertices(ring int, fn func(float64, float64)) {
+	if ring < 0 || ring >= len(p.rings) {
+		return
+	}
+	rng := p.rings[ring]
+	for i := 0; i < len(rng); i += 2 {
+		fn(rng[i+0], rng[i+1])
+	}
+}
+
+func (p *Polygon) AsWKB() []byte {
+
+	size := 1 + 4 + 4
+	for _, r := range p.rings {
+		size += 4 + len(r)*8
+	}
+
+	buf := bytes.NewBuffer(make([]byte, 0, size))
+
+	binary.Write(buf, binary.LittleEndian, wkb.NDR)
+	binary.Write(buf, binary.LittleEndian, wkb.Polygon)
+	binary.Write(buf, binary.LittleEndian, uint32(len(p.rings)))
+
+	for _, r := range p.rings {
+		binary.Write(buf, binary.LittleEndian, uint32(len(r)/2))
+		for i := 0; i < len(r); i += 2 {
+			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0]))
+			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1]))
+		}
+	}
+
+	return buf.Bytes()
+}
+
 func (p *Polygon) FromWKB(data []byte) error {
 
 	r := bytes.NewReader(data)
@@ -477,3 +518,14 @@
 
 	return nil
 }
+
+func (p *Polygon) FromLineStringZ(ls LineStringZ) {
+	r := make([]float64, 2*len(ls))
+	var pos int
+	for i := range ls {
+		r[pos+0] = ls[i].X
+		r[pos+1] = ls[i].Y
+		pos += 2
+	}
+	p.rings = []ring{r}
+}
diff -r 7126fdc5779a -r b30c0b02e8fb pkg/octree/tree.go
--- a/pkg/octree/tree.go	Fri Jun 14 11:16:09 2019 +0200
+++ b/pkg/octree/tree.go	Fri Jun 14 12:03:47 2019 +0200
@@ -15,6 +15,10 @@
 
 import (
 	"math"
+	"runtime"
+	"sync"
+
+	"gemma.intevation.de/gemma/pkg/common"
 )
 
 // Tree is an Octree holding triangles.
@@ -316,3 +320,65 @@
 
 	return result
 }
+
+func (t *Tree) GenerateRandomVertices(n int, callback func([]Vertex)) {
+	var wg sync.WaitGroup
+
+	jobs := make(chan int)
+	out := make(chan []Vertex)
+	done := make(chan struct{})
+
+	cpus := runtime.NumCPU()
+
+	free := make(chan []Vertex, cpus)
+
+	for i := 0; i < cpus; i++ {
+		wg.Add(1)
+		go func() {
+			defer wg.Done()
+			xRange := common.Random(t.Min.X, t.Max.X)
+			yRange := common.Random(t.Min.Y, t.Max.Y)
+
+			for size := range jobs {
+				var vertices []Vertex
+				select {
+				case vertices = <-free:
+				default:
+					vertices = make([]Vertex, 0, 1000)
+				}
+				for len(vertices) < size {
+					x, y := xRange(), yRange()
+					if z, ok := t.Value(x, y); ok {
+						vertices = append(vertices, Vertex{X: x, Y: y, Z: z})
+					}
+				}
+				out <- vertices
+			}
+		}()
+	}
+
+	go func() {
+		defer close(done)
+		for vertices := range out {
+			callback(vertices)
+			select {
+			case free <- vertices[:0]:
+			default:
+			}
+		}
+	}()
+
+	for remain := n; remain > 0; {
+		if remain > 1000 {
+			jobs <- 1000
+			remain -= 1000
+		} else {
+			jobs <- remain
+			remain = 0
+		}
+	}
+	close(jobs)
+	wg.Wait()
+	close(out)
+	<-done
+}
diff -r 7126fdc5779a -r b30c0b02e8fb pkg/octree/triangulation.go
--- a/pkg/octree/triangulation.go	Fri Jun 14 11:16:09 2019 +0200
+++ b/pkg/octree/triangulation.go	Fri Jun 14 12:03:47 2019 +0200
@@ -39,11 +39,11 @@
 	return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err
 }
 
-func (t *Triangulation) ConcaveHull(tooLong float64) []int32 {
+func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) {
 
 	tooLong *= tooLong
 
-	var oldCandidates []int32
+	var candidates []int32
 
 	for i, num := 0, len(t.Triangles)/3; i < num; i++ {
 		idx := i * 3
@@ -56,17 +56,21 @@
 			}
 		}
 		if max > tooLong {
-			oldCandidates = append(oldCandidates, int32(i))
+			candidates = append(candidates, int32(i))
 		}
 	}
 
-	removed := map[int32]bool{}
+	removed := map[int32]struct{}{}
 
 	isBorder := func(n int32) bool {
 		n *= 3
 		for i := int32(0); i < 3; i++ {
 			e := n + i
-			if t.Halfedges[e] == -1 || removed[e] {
+			o := t.Halfedges[e]
+			if o < 0 {
+				return true
+			}
+			if _, found := removed[o/3]; found {
 				return true
 			}
 		}
@@ -75,15 +79,15 @@
 
 	var newCandidates []int32
 
-	for len(oldCandidates) > 0 {
-		log.Printf("candidates: %d\n", len(oldCandidates))
+	log.Printf("info: candidates: %d\n", len(candidates))
+	for len(candidates) > 0 {
 
 		oldRemoved := len(removed)
 
-		for _, i := range oldCandidates {
+		for _, i := range candidates {
 
 			if isBorder(i) {
-				removed[i] = true
+				removed[i] = struct{}{}
 			} else {
 				newCandidates = append(newCandidates, i)
 			}
@@ -93,13 +97,110 @@
 			break
 		}
 
-		oldCandidates = newCandidates
+		candidates = newCandidates
 		newCandidates = newCandidates[:0]
 	}
 
+	log.Printf("info: candidates left: %d\n", len(candidates))
+	log.Printf("info: triangles: %d\n", len(t.Triangles)/3)
 	log.Printf("info: triangles to remove: %d\n", len(removed))
 
-	return nil
+	type edge struct {
+		a, b       int32
+		prev, next *edge
+	}
+
+	isClosed := func(e *edge) bool {
+		for curr := e.next; curr != nil; curr = curr.next {
+			if curr == e {
+				return true
+			}
+		}
+		return false
+	}
+
+	open := map[int32]*edge{}
+	var rings []*edge
+
+	for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ {
+		if _, found := removed[i]; found {
+			continue
+		}
+		n := i * 3
+		for j := int32(0); j < 3; j++ {
+			e := n + j
+			f := t.Halfedges[e]
+			if f >= 0 {
+				if _, found := removed[f/3]; !found {
+					continue
+				}
+			}
+			a := t.Triangles[e]
+			b := t.Triangles[n+(j+1)%3]
+
+			curr := &edge{a: a, b: b}
+
+			if old := open[a]; old != nil {
+				delete(open, a)
+				if old.a == a {
+					old.prev = curr
+					curr.next = old
+				} else {
+					old.next = curr
+					curr.prev = old
+				}
+
+				if isClosed(old) {
+					rings = append(rings, old)
+				}
+			} else {
+				open[a] = curr
+			}
+
+			if old := open[b]; old != nil {
+				delete(open, b)
+				if old.b == b {
+					old.next = curr
+					curr.prev = old
+				} else {
+					old.prev = curr
+					curr.next = old
+				}
+
+				if isClosed(old) {
+					rings = append(rings, old)
+				}
+			} else {
+				open[b] = curr
+			}
+		}
+	}
+
+	if len(open) > 0 {
+		log.Printf("warn: open vertices left: %d\n", len(open))
+	}
+
+	if len(rings) == 0 {
+		log.Println("warn: no ring found")
+		return nil, removed
+	}
+
+	curr := rings[0]
+
+	polygon := LineStringZ{
+		t.Points[curr.a],
+		t.Points[curr.b],
+	}
+
+	for curr = curr.next; curr != rings[0]; curr = curr.next {
+		polygon = append(polygon, t.Points[curr.b])
+	}
+
+	polygon = append(polygon, t.Points[rings[0].a])
+
+	log.Printf("length of boundary: %d\n", len(polygon))
+
+	return polygon, removed
 }
 
 func (t *Triangulation) TriangleSlices() [][]int32 {
diff -r 7126fdc5779a -r b30c0b02e8fb pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Fri Jun 14 11:16:09 2019 +0200
+++ b/pkg/octree/vertex.go	Fri Jun 14 12:03:47 2019 +0200
@@ -210,6 +210,11 @@
 	return math.Hypot(v.X-w.X, v.Y-w.Y)
 }
 
+func (v Vertex) Distance(w Vertex) float64 {
+	v = v.Sub(w)
+	return math.Sqrt(v.Dot(v))
+}
+
 // Minimize adjust this vertex v to hold the minimum
 // values component-wise of v and w.
 func (v *Vertex) Minimize(w Vertex) {
@@ -413,6 +418,34 @@
 	}
 }
 
+func (ls LineStringZ) BBox() Box2D {
+
+	min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64}
+	max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64}
+
+	for _, v := range ls {
+		min.Minimize(v)
+		max.Maximize(v)
+	}
+
+	return Box2D{
+		X1: min.X,
+		Y1: min.Y,
+		X2: max.X,
+		Y2: max.Y,
+	}
+}
+
+func (ls LineStringZ) Area() float64 {
+	return polygonArea(ls)
+}
+
+func (ls LineStringZ) Reverse() {
+	for i, j := 0, len(ls)-1; i < j; i, j = i+1, j-1 {
+		ls[i], ls[j] = ls[j], ls[i]
+	}
+}
+
 func (ls LineStringZ) order(position func(Vertex) float64) {
 	type posVertex struct {
 		pos float64
@@ -738,6 +771,10 @@
 	}
 }
 
+func (a Box2D) Area() float64 {
+	return (a.X2 - a.X1) * (a.Y2 - a.Y1)
+}
+
 // NewPlane2D creates a new Plane2D from two given points.
 func NewPlane2D(x1, y1, x2, y2 float64) Plane2D {
 	b := x2 - x1