gemma: pkg/imports/sr.go comparison

comparison pkg/imports/sr.go @ 3646:810b28f59b8b single-beam

Generate random points for second mesh.

author	Sascha L. Teichmann <sascha.teichmann@intevation.de>
date	Wed, 12 Jun 2019 17:30:20 +0200
parents	e1021fd60190
children	01ce3ba9b0d0

comparison

equal deleted inserted replaced

-:10471aa73ad8
+:810b28f59b8b
 	"database/sql"
 	"encoding/hex"
 	"errors"
 	"fmt"
 	"io"
+	"log"
 	"math"
 	"os"
 	"path"
 	"path/filepath"
+	"runtime"
 	"strconv"
 	"strings"
+	"sync"
 	"time"
 	shp "github.com/jonas-p/go-shp"
 	"gemma.intevation.de/gemma/pkg/common"
 }
 const (
 	contourStepWidth = 0.1
 	contourTolerance = 0.1
+)
+const (
+	tooLongEdge          = 50.0
+	pointsPerSquareMeter = 5
 )
 // SRJobKind is the unique name of a SoundingResult import job.
 const SRJobKind JobKind = "sr"
 		return nil, err
 	}
 	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)
-		tri.ConcaveHull(50.0)
+		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)
+				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!
 	return nil, errors.New("Not implemented, yet!")
 }
 	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))

Mercurial > gemma

comparison pkg/imports/sr.go @ 3646:810b28f59b8b single-beam