Mercurial > gemma

changeset 3658:1c3df921361d single-beam

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Handle th case that a boundary polygon is uploaded along side with the single beam scan.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 14 Jun 2019 11:59:14 +0200
parents a6c671abbc35
children 66f2cb789905
files pkg/imports/sr.go pkg/octree/polygon.go
diffstat 2 files changed, 215 insertions(+), 155 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/imports/sr.go	Thu Jun 13 18:33:54 2019 +0200
+++ b/pkg/imports/sr.go	Fri Jun 14 11:59:14 2019 +0200
@@ -136,8 +136,14 @@
 `
 
 	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
@@ -338,6 +344,7 @@
 ) (interface{}, error) {
 
 	feedback.Info("Processing as single beam scan.")
+	feedback.Info("Reproject XYZ data.")
 
 	start := time.Now()
 
@@ -361,6 +368,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 +378,23 @@
 	feedback.Info("Triangulation took %v.", time.Since(start))
 	feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
 
-	var clippingPolygon octree.Polygon
+	var (
+		clippingPolygon         octree.Polygon
+		clippingPolygonBuffered octree.Polygon
+		removed                 map[int32]struct{}
+		polygonArea             float64
+		clippingPolygonWKB      []byte
+	)
 
 	if boundary == nil {
-		feedback.Info("No boundary given. Calulating from XYZ data.")
+		feedback.Info("No boundary given. Calulate from XYZ data.")
 		feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
 
-		polygon, removed := tri.ConcaveHull(tooLongEdge)
+		var polygon octree.LineStringZ
+		start = time.Now()
+		polygon, removed = tri.ConcaveHull(tooLongEdge)
 
-		polygonArea := polygon.Area()
+		polygonArea = polygon.Area()
 		if polygonArea < 0.0 { // counter clockwise
 			polygonArea = -polygonArea
 			polygon.Reverse()
@@ -386,178 +402,205 @@
 
 		clippingPolygon.FromLineStringZ(polygon)
 
-		var repaired, buffered []byte
-
+		var buffered []byte
 		if err := tx.QueryRowContext(
 			ctx,
 			repairBoundarySQL,
 			clippingPolygon.AsWKB(),
 			epsg,
-		).Scan(&repaired, &buffered); err != nil {
+		).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
 		}
 
-		if err := clippingPolygon.FromWKB(repaired); err != nil {
+		feedback.Info("Calculating took %v.", time.Since(start))
+
+	} 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
 		}
-		var clippingPolygonBuffered octree.Polygon
+		if err := clippingPolygon.FromWKB(clippingPolygonWKB); err != nil {
+			return nil, err
+		}
 		if err := clippingPolygonBuffered.FromWKB(buffered); err != nil {
 			return nil, err
 		}
 
-		firstTin := tri.Tin()
-		builder := octree.NewBuilder(firstTin)
-		builder.Build(removed)
-
-		firstTree := builder.Tree()
-
-		feedback.Info("Boundary area: %.2fm²", polygonArea)
-
-		numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
-
-		feedback.Info("Generating %d random points for an average density of ~%d points/m².",
-			numPoints, pointsPerSquareMeter)
-
-		start := time.Now()
-
-		generated := make(octree.LineStringZ, 0, numPoints+len(polygon)-1)
-
-		firstTree.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.
-		generated = append(generated, polygon[:len(polygon)-1]...)
-
-		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()
-		secondTin := tri.Tin()
-		secondTin.EPSG = uint32(epsg)
-		clippingPolygonBuffered.Indexify()
+		tin := tri.Tin()
+		tin.EPSG = uint32(epsg)
 
 		var str octree.STRTree
-		str.Build(secondTin)
-		feedback.Info("Building STR tree took %v", time.Since(start))
-
-		start = time.Now()
-
-		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("Building final octree index")
-
-		builder = octree.NewBuilder(secondTin)
-		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()
+		str.Build(tin)
 
-		var (
-			id       int64
-			dummy    uint32
-			lat, lon float64
-		)
-
-		var hull []byte
+		removed = str.Clip(&clippingPolygon)
+	}
 
-		if err := tx.QueryRowContext(
-			ctx,
-			insertHullSQL,
-			m.Bottleneck,
-			m.Date.Time,
-			m.DepthReference,
-			nil,
-			repaired,
-			epsg,
-		).Scan(
-			&id,
-			&lat,
-			&lon,
-			&dummy,
-			&hull,
-		); err != nil {
-			return nil, err
-		}
+	// Build the first mesh to generate random points on.
 
-		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()
+	feedback.Info("Build virtual DEM based on original XYZ data.")
 
-		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
-		}
+	start = time.Now()
 
-		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
-
-	} 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()
+	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(
--- a/pkg/octree/polygon.go	Thu Jun 13 18:33:54 2019 +0200
+++ b/pkg/octree/polygon.go	Fri Jun 14 11:59:14 2019 +0200
@@ -409,6 +409,23 @@
 	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