Mercurial > gemma
diff cmd/octreediff/main.go @ 2484:4fa92d468164 octree-diff
Use fogleman's triangulation.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 01 Mar 2019 17:42:46 +0100 |
parents | 620038ade708 |
children | cb55d7eaaa36 |
line wrap: on
line diff
--- a/cmd/octreediff/main.go Fri Mar 01 15:33:27 2019 +0100 +++ b/cmd/octreediff/main.go Fri Mar 01 17:42:46 2019 +0100 @@ -14,22 +14,18 @@ package main import ( - "bytes" "context" "database/sql" - "encoding/binary" "errors" "flag" "fmt" "log" - "math" "runtime" "sync" "time" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/octree" - "gemma.intevation.de/gemma/pkg/wkb" ) var ( @@ -48,7 +44,28 @@ WHERE bn.bottleneck_id = $1 AND sr.date_info = $2::date AND sr.octree_index IS NOT NULL` - triangulateSQL = ` + insertContourSQL = ` +INSERT INTO diff_contour_lines ( + height, + lines +) +SELECT + $4, + ST_Transform( + ST_Multi( + ST_CollectionExtract( + ST_SimplifyPreserveTopology( + ST_Multi(ST_Collectionextract( + ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 2)), + $3 + ), + 2 + ) + ), + 4326 + ) +` + insertContourSQLClipped = ` WITH joined AS ( SELECT sr.area AS area, @@ -65,40 +82,30 @@ (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date) ), 0.001) AS area -), -triangles AS ( - SELECT t.geom AS geom, ST_MakePolygon(ST_ExteriorRing(t.geom)) AS poly FROM ( - SELECT (ST_Dump( - ST_DelaunayTriangles(ST_GeomFromWKB($5, $2::int), 0, 2))).geom - ) t ) -SELECT ST_AsBinary(ST_Collect(triangles.geom)) FROM triangles, inter -WHERE ST_Covers(inter.area, triangles.poly) -` - - //INSERT INTO redis_diff_countour_lines ( - //INSERT INTO diff_contour_lines ( - //INSERT INTO diff_contour_lines_extern ( - insertContourSQL = ` INSERT INTO diff_contour_lines ( height, lines ) SELECT - $1, + $5, ST_Transform( ST_Multi( - ST_CollectionExtract( - ST_SimplifyPreserveTopology( - ST_Multi(ST_Collectionextract( - ST_MakeValid(ST_GeomFromWKB($2, $3::integer)), 2)), - $4 - ), - 2 - ) - ), + ST_intersection( + ST_CollectionExtract( + ST_SimplifyPreserveTopology( + ST_Multi(ST_Collectionextract( + ST_MakeValid(ST_GeomFromWKB($6, $2::integer)), 2)), + $7 + ), + 2 + ), + area + ) + ), 4326 - )` + ) + FROM inter` ) func check(err error) { @@ -114,43 +121,14 @@ type pointMap map[point]float64 -func (pm pointMap) triangulate() { - start := time.Now() - points := make([]octree.Vertex, len(pm)) +func (pm pointMap) triangulate() (*octree.Triangulation, error) { + vertices := make([]octree.Vertex, len(pm)) var i int for p, z := range pm { - points[i] = octree.Vertex{X: p.x, Y: p.y, Z: z} + vertices[i] = octree.Vertex{X: p.x, Y: p.y, Z: z} i++ } - _, err := octree.Triangulate(points) - if err != nil { - log.Printf("triangulate error: %v\n", err) - } - log.Printf("in memory triangulation (%d points) took %s\n", i, time.Since(start)) -} - -func (pm pointMap) asWKB() []byte { - size := 1 + 4 + 4 + len(pm)*(1+4+3*8) - - buf := bytes.NewBuffer(make([]byte, 0, size)) - - binary.Write(buf, binary.LittleEndian, wkb.NDR) - binary.Write(buf, binary.LittleEndian, wkb.MultiPointZ) - binary.Write(buf, binary.LittleEndian, uint32(len(pm))) - - perPoint := bytes.NewBuffer(make([]byte, 0, 1+4)) - binary.Write(perPoint, binary.LittleEndian, wkb.NDR) - binary.Write(perPoint, binary.LittleEndian, wkb.PointZ) - hdr := perPoint.Bytes() - - for p, z := range pm { - buf.Write(hdr) - binary.Write(buf, binary.LittleEndian, math.Float64bits(p.x)) - binary.Write(buf, binary.LittleEndian, math.Float64bits(p.y)) - binary.Write(buf, binary.LittleEndian, math.Float64bits(z)) - } - - return buf.Bytes() + return octree.Triangulate(vertices) } func sliceWork( @@ -245,6 +223,7 @@ if err != nil { return err } + defer tx.Rollback() type load struct { date time.Time @@ -368,27 +347,8 @@ last = now log.Printf("num points: %d\n", len(result)) - result.triangulate() - - data := result.asWKB() - - now = time.Now() - log.Printf("turing into WKB took %v\n", now.Sub(last)) - last = now - - log.Printf("WKB size %.3fMB\n", float64(len(data))/(1024*1024)) - - var tin octree.Tin - - if err := tx.QueryRowContext( - ctx, - triangulateSQL, - bottleneck, - first.EPSG, - firstDate, - secondDate, - data, - ).Scan(&tin); err != nil { + tri, err := result.triangulate() + if err != nil { return err } @@ -396,7 +356,7 @@ log.Printf("triangulation took %v\n", now.Sub(last)) last = now - builder := octree.NewBuilder(&tin) + builder := octree.NewBuilder(tri.Tin()) builder.Build() now = time.Now() @@ -434,24 +394,30 @@ } } + log.Printf("num heights: %d\n", len(heights)) + var dataSize int - stmt, err := tx.PrepareContext(ctx, insertContourSQL) + stmt, err := tx.PrepareContext(ctx, insertContourSQLClipped) if err != nil { return err } defer stmt.Close() + log.Println("do contours") octree.DoContours(tree, heights, func(res *octree.ContourResult) { if err == nil && len(res.Lines) > 0 { - log.Printf("%f: lines: %d\n", res.Height, len(res.Lines)) + log.Printf("%.2f: lines: %d\n", res.Height, len(res.Lines)) wkb := res.Lines.AsWKB2D() dataSize += len(wkb) _, err = stmt.ExecContext( ctx, + bottleneck, + first.EPSG, + firstDate, + secondDate, res.Height, wkb, - first.EPSG, contourTolerance, ) }