Mercurial > gemma

changeset 3655:a6c671abbc35 single-beam

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Started with cleaning up the single beam import.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 13 Jun 2019 18:33:54 +0200
parents 26e02fa5b992
children 1c3df921361d
files pkg/imports/sr.go pkg/octree/tree.go pkg/octree/vertex.go
diffstat 3 files changed, 90 insertions(+), 131 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/imports/sr.go	Thu Jun 13 17:07:28 2019 +0200
+++ b/pkg/imports/sr.go	Thu Jun 13 18:33:54 2019 +0200
@@ -24,15 +24,12 @@
 	"errors"
 	"fmt"
 	"io"
-	"log"
 	"math"
 	"os"
 	"path"
 	"path/filepath"
-	"runtime"
 	"strconv"
 	"strings"
-	"sync"
 	"time"
 
 	shp "github.com/jonas-p/go-shp"
@@ -151,14 +148,6 @@
   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))`
 
-	whatsWrongSQL = `
-SELECT
-  ST_IsValidReason(
-    CAST(ST_MakeValid(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326))::geography AS geometry)),
-  ST_IsValid(
-    CAST(ST_MakeValid(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326))::geography AS geometry))
-`
-
 	insertContourSQL = `
 INSERT INTO waterway.sounding_results_contour_lines (
   sounding_result_id,
@@ -384,7 +373,8 @@
 	var clippingPolygon octree.Polygon
 
 	if boundary == nil {
-		feedback.Info("No boundary given.")
+		feedback.Info("No boundary given. Calulating from XYZ data.")
+		feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
 
 		polygon, removed := tri.ConcaveHull(tooLongEdge)
 
@@ -396,19 +386,6 @@
 
 		clippingPolygon.FromLineStringZ(polygon)
 
-		var wrong string
-		var valid bool
-		if err := tx.QueryRowContext(
-			ctx,
-			whatsWrongSQL,
-			clippingPolygon.AsWKB(),
-			epsg,
-		).Scan(&wrong, &valid); err != nil {
-			return nil, err
-		}
-
-		log.Printf("!!!!!!!!: whats wrong (%t): %s\n", valid, wrong)
-
 		var repaired, buffered []byte
 
 		if err := tx.QueryRowContext(
@@ -434,115 +411,58 @@
 
 		firstTree := builder.Tree()
 
-		bbox := polygon.BBox()
-		bboxArea := bbox.Area()
-
-		log.Printf("bbox area: %.2f\n", bboxArea)
-		log.Printf("polygon area: %.2f\n", polygonArea)
+		feedback.Info("Boundary area: %.2fm²", polygonArea)
 
 		numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
 
-		log.Printf("random points for density of ~%d points/m²: %d\n",
-			pointsPerSquareMeter, numPoints)
+		feedback.Info("Generating %d random points for an average density of ~%d points/m².",
+			numPoints, pointsPerSquareMeter)
 
 		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)
-
-				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:
-				}
-			}
-		}()
+		firstTree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
+			generated = append(generated, vertices...)
+		})
 
-		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))
+		feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
 
 		// Add the boundary to new point cloud.
 		generated = append(generated, polygon[:len(polygon)-1]...)
 
+		feedback.Info("Triangulate new point cloud.")
+
+		start = time.Now()
+
 		xyz = octree.MultiPointZ(generated)
 
-		start = time.Now()
 		tri, err := octree.Triangulate(xyz)
 		if err != nil {
 			return nil, err
 		}
-		log.Printf("Second triangulation took %v.", time.Since(start))
-		log.Printf("Number triangles: %d.", len(tri.Triangles)/3)
+		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()
 		secondTin := tri.Tin()
-
-		//polygon.Buffer(1.0)
-
+		secondTin.EPSG = uint32(epsg)
 		clippingPolygonBuffered.Indexify()
 
 		var str octree.STRTree
 		str.Build(secondTin)
-		log.Printf("Building STR tree took %v", time.Since(start))
+		feedback.Info("Building STR tree took %v", time.Since(start))
 
 		start = time.Now()
 
 		removed = str.Clip(&clippingPolygonBuffered)
-		log.Printf("Clipping STR tree took %v.", time.Since(start))
-		log.Printf("Number of triangles to clip %d.", len(removed))
+		feedback.Info("Clipping STR tree took %v.", time.Since(start))
+		feedback.Info("Number of triangles to clip %d.", len(removed))
 
 		start = time.Now()
 
-		secondTin.EPSG = uint32(epsg)
+		feedback.Info("Building final octree index")
 
 		builder = octree.NewBuilder(secondTin)
 		builder.Build(removed)
@@ -552,6 +472,8 @@
 		}
 		feedback.Info("Building octree took %v.", time.Since(start))
 
+		start = time.Now()
+
 		var (
 			id       int64
 			dummy    uint32
@@ -579,7 +501,6 @@
 			return nil, err
 		}
 
-		start = time.Now()
 		h := sha1.New()
 		h.Write(octreeIndex)
 		checksum := hex.EncodeToString(h.Sum(nil))
--- a/pkg/octree/tree.go	Thu Jun 13 17:07:28 2019 +0200
+++ b/pkg/octree/tree.go	Thu Jun 13 18:33:54 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
+}
--- a/pkg/octree/vertex.go	Thu Jun 13 17:07:28 2019 +0200
+++ b/pkg/octree/vertex.go	Thu Jun 13 18:33:54 2019 +0200
@@ -446,34 +446,6 @@
 	}
 }
 
-func (ls LineStringZ) Buffer(radius float64) {
-	if len(ls) == 0 {
-		return
-	}
-	var cx, cy, cz float64
-	for i := range ls {
-		cx += ls[i].X
-		cy += ls[i].Y
-		cz += ls[i].Z
-	}
-
-	s := 1.0 / float64(len(ls))
-
-	cx *= s
-	cy *= s
-	cz *= s
-	c := Vertex{X: cx, Y: cy, Z: cz}
-	for i := range ls {
-		n := c.Sub(ls[i]).Normalize()
-		l := ls[i].Distance(c)
-		// l = 1
-		// l + r = s
-		// (l+r)/l = s
-		s := (l + radius) / l
-		ls[i] = ls[i].Add(n.Scale(s))
-	}
-}
-
 func (ls LineStringZ) order(position func(Vertex) float64) {
 	type posVertex struct {
 		pos float64