# HG changeset patch # User Sascha L. Teichmann # Date 1560506354 -7200 # Node ID 1c3df921361d66efbe2890c4e3aec69683f999e2 # Parent a6c671abbc350244588550d1312886df46ec6f2e Handle th case that a boundary polygon is uploaded along side with the single beam scan. diff -r a6c671abbc35 -r 1c3df921361d pkg/imports/sr.go --- 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( diff -r a6c671abbc35 -r 1c3df921361d pkg/octree/polygon.go --- 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