# HG changeset patch # User Sascha L. Teichmann # Date 1560443634 -7200 # Node ID a6c671abbc350244588550d1312886df46ec6f2e # Parent 26e02fa5b992905f4340fe3d63f866bda5e0ced0 Started with cleaning up the single beam import. diff -r 26e02fa5b992 -r a6c671abbc35 pkg/imports/sr.go --- 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)) diff -r 26e02fa5b992 -r a6c671abbc35 pkg/octree/tree.go --- 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 +} diff -r 26e02fa5b992 -r a6c671abbc35 pkg/octree/vertex.go --- 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