Mercurial > gemma

changeset 3744:16c3c0150030

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge concave-hull branch into default.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 24 Jun 2019 15:22:39 +0200
parents 71de5ce5a402 (current diff) fcac7787ec22 (diff)
children ae79c9546908
files
diffstat 2 files changed, 115 insertions(+), 221 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/imports/sr.go	Mon Jun 24 12:47:34 2019 +0200
+++ b/pkg/imports/sr.go	Mon Jun 24 15:22:39 2019 +0200
@@ -68,7 +68,6 @@
 )
 
 const (
-	tooLongEdge          = 50.0
 	pointsPerSquareMeter = 2
 )
 
@@ -307,20 +306,22 @@
 	var summary interface{}
 
 	if sr.isSingleBeam() {
-		summary, err = sr.singleBeamScan(
+		summary, err = sr.processScan(
 			ctx,
 			tx,
 			feedback,
+			true,
 			importID,
 			m,
 			xyz,
 			boundary,
 		)
 	} else {
-		summary, err = sr.multiBeamScan(
+		summary, err = sr.processScan(
 			ctx,
 			tx,
 			feedback,
+			false,
 			importID,
 			m,
 			xyz,
@@ -343,24 +344,30 @@
 	return summary, nil
 }
 
-func (sr *SoundingResult) singleBeamScan(
+func (sr *SoundingResult) processScan(
 	ctx context.Context,
 	tx *sql.Tx,
 	feedback Feedback,
+	singleBeam bool,
 	importID int64,
 	m *models.SoundingResultMeta,
 	xyz octree.MultiPointZ,
 	boundary polygonSlice,
 ) (interface{}, error) {
 
-	feedback.Info("Processing as single beam scan.")
+	if singleBeam {
+		feedback.Info("Processing as single beam scan.")
+	} else {
+		feedback.Info("Processing as multi beam scan.")
+	}
+
 	feedback.Info("Reproject XYZ data.")
 
 	start := time.Now()
 
 	xyzWKB := xyz.AsWKB()
 	var reproj []byte
-	var epsg uint
+	var epsg uint32
 
 	if err := tx.QueryRowContext(
 		ctx,
@@ -394,10 +401,12 @@
 		removed                 map[int32]struct{}
 		polygonArea             float64
 		clippingPolygonWKB      []byte
+		tin                     *octree.Tin
 	)
 
 	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
@@ -451,7 +460,7 @@
 		}
 
 		tin := tri.Tin()
-		tin.EPSG = uint32(epsg)
+		tin.EPSG = epsg
 
 		var str octree.STRTree
 		str.Build(tin)
@@ -459,68 +468,72 @@
 		removed = str.Clip(&clippingPolygon)
 	}
 
-	// Build the first mesh to generate random points on.
+	if singleBeam {
+
+		// Build the first mesh to generate random points on.
+
+		feedback.Info("Build virtual DEM based on original XYZ data.")
+
+		start = time.Now()
 
-	feedback.Info("Build virtual DEM based on original XYZ data.")
+		var tree *octree.Tree
+		{
+			builder := octree.NewBuilder(tri.Tin())
+			builder.Build(removed)
+			tree = builder.Tree()
+		}
+
+		feedback.Info("Building took %v", time.Since(start))
 
-	start = time.Now()
+		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))
 
-	var tree *octree.Tree
-	{
-		builder := octree.NewBuilder(tri.Tin())
-		builder.Build(removed)
-		tree = builder.Tree()
-	}
+		tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
+			generated = append(generated, vertices...)
+		})
+
+		feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
 
-	feedback.Info("Building took %v", time.Since(start))
-
-	feedback.Info("Boundary area: %.2fm²", polygonArea)
+		// 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})
+			}
+		})
 
-	numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
+		feedback.Info("Triangulate new point cloud.")
+		xyz = octree.MultiPointZ(generated)
+		start = time.Now()
 
-	feedback.Info("Generate %d random points for an average density of ~%d points/m².",
-		numPoints, pointsPerSquareMeter)
+		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.")
+
+	} else { // multi beam
+		// Nothing special
+	}
 
 	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)
+	tin = tri.Tin()
+	tin.EPSG = epsg
 
 	var str octree.STRTree
 	str.Build(tin)
@@ -533,8 +546,6 @@
 	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)
@@ -619,160 +630,6 @@
 	return &summary, nil
 }
 
-func (sr *SoundingResult) multiBeamScan(
-	ctx context.Context,
-	tx *sql.Tx,
-	feedback Feedback,
-	importID int64,
-	m *models.SoundingResultMeta,
-	xyz octree.MultiPointZ,
-	boundary polygonSlice,
-) (interface{}, error) {
-
-	feedback.Info("Processing as multi beam scan.")
-	var (
-		id       int64
-		epsg     uint32
-		lat, lon float64
-	)
-	start := time.Now()
-
-	var hull []byte
-
-	xyzWKB := xyz.AsWKB()
-
-	err := tx.QueryRowContext(
-		ctx,
-		insertHullSQL,
-		m.Bottleneck,
-		m.Date.Time,
-		m.DepthReference,
-		xyzWKB,
-		boundary.asWKB(),
-		m.EPSG,
-	).Scan(
-		&id,
-		&lat,
-		&lon,
-		&epsg,
-		&hull,
-	)
-	xyz, boundary = nil, nil // not need from now on.
-	feedback.Info("Calculating hull took %s.", time.Since(start))
-	switch {
-	case err == sql.ErrNoRows:
-		return nil, fmt.Errorf(
-			"No matching bottleneck of given name or time available: %v", err)
-	case err != nil:
-		return nil, err
-	}
-	feedback.Info("Best suited UTM EPSG: %d", epsg)
-
-	start = time.Now()
-
-	var clippingPolygon octree.Polygon
-
-	if err := clippingPolygon.FromWKB(hull); err != nil {
-		return nil, err
-	}
-	clippingPolygon.Indexify()
-	feedback.Info("Building clipping polygon took %v.", time.Since(start))
-
-	start = time.Now()
-
-	var reproj []byte
-
-	if err = tx.QueryRowContext(
-		ctx,
-		reprojectPointsSQL,
-		xyzWKB,
-		m.EPSG,
-		epsg,
-	).Scan(&reproj); err != nil {
-		return nil, err
-	}
-
-	if err := xyz.FromWKB(reproj); err != nil {
-		return nil, err
-	}
-
-	feedback.Info("Reprojecting points took %v.", time.Since(start))
-	feedback.Info("Number of reprojected points: %d", len(xyz))
-
-	start = time.Now()
-
-	tri, err := octree.Triangulate(xyz)
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Triangulation took %v.", time.Since(start))
-
-	start = time.Now()
-
-	tin := tri.Tin()
-
-	var str octree.STRTree
-	str.Build(tin)
-	feedback.Info("Building STR tree took %v", time.Since(start))
-
-	start = time.Now()
-
-	removed := str.Clip(&clippingPolygon)
-	feedback.Info("Clipping STR tree took %v.", time.Since(start))
-	feedback.Info("Number of triangles to clip %d.", len(removed))
-
-	start = time.Now()
-
-	tin.EPSG = epsg
-
-	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))
-
-	start = time.Now()
-	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))
-
-	tree := builder.Tree()
-
-	start = time.Now()
-	err = generateContours(ctx, tx, tree, id)
-	if err != nil {
-		return nil, err
-	}
-	feedback.Info("Generating and storing 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
-}
-
 // CleanUp removes the folder containing the ZIP file with the
 // the sounding result import.
 func (sr *SoundingResult) CleanUp() error {
--- a/pkg/octree/triangulation.go	Mon Jun 24 12:47:34 2019 +0200
+++ b/pkg/octree/triangulation.go	Mon Jun 24 15:22:39 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])