# HG changeset patch # User Sascha L. Teichmann # Date 1551458566 -3600 # Node ID 4fa92d4681642a3bf034d02cd923eb3924f734ce # Parent 620038ade708923c4411adb2f7ab36cd7604c75d Use fogleman's triangulation. diff -r 620038ade708 -r 4fa92d468164 cmd/octreediff/main.go --- 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, ) } diff -r 620038ade708 -r 4fa92d468164 pkg/octree/node.go --- a/pkg/octree/node.go Fri Mar 01 15:33:27 2019 +0100 +++ b/pkg/octree/node.go Fri Mar 01 17:42:46 2019 +0100 @@ -20,13 +20,13 @@ package octree type node struct { - i int - t int + i int32 + t int32 prev *node next *node } -func newNode(nodes []node, i int, prev *node) *node { +func newNode(nodes []node, i int32, prev *node) *node { n := &nodes[i] n.i = i if prev == nil { diff -r 620038ade708 -r 4fa92d468164 pkg/octree/triangulation.go --- a/pkg/octree/triangulation.go Fri Mar 01 15:33:27 2019 +0100 +++ b/pkg/octree/triangulation.go Fri Mar 01 17:42:46 2019 +0100 @@ -27,8 +27,8 @@ type Triangulation struct { Points []Vertex ConvexHull []Vertex - Triangles []int - Halfedges []int + Triangles []int32 + Halfedges []int32 } // Triangulate returns a Delaunay triangulation of the provided points. @@ -38,6 +38,36 @@ return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err } +func (t *Triangulation) TriangleSlices() [][]int32 { + sl := make([][]int32, len(t.Triangles)/3) + var j int + for i := range sl { + sl[i] = t.Triangles[j : j+3] + j += 3 + } + return sl +} + +func (t *Triangulation) Tin() *Tin { + + min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} + max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} + + vertices := t.Points + + for _, v := range vertices { + min.Minimize(v) + max.Maximize(v) + } + + return &Tin{ + Vertices: vertices, + Triangles: t.TriangleSlices(), + Min: min, + Max: max, + } +} + func (t *Triangulation) area() float64 { var result float64 points := t.Points @@ -60,7 +90,7 @@ if i1 != -1 && t.Halfedges[i1] != i2 { return fmt.Errorf("invalid halfedge connection") } - if i2 != -1 && t.Halfedges[i2] != i1 { + if i2 != -1 && t.Halfedges[i2] != int32(i1) { return fmt.Errorf("invalid halfedge connection") } } diff -r 620038ade708 -r 4fa92d468164 pkg/octree/triangulator.go --- a/pkg/octree/triangulator.go Fri Mar 01 15:33:27 2019 +0100 +++ b/pkg/octree/triangulator.go Fri Mar 01 17:42:46 2019 +0100 @@ -28,10 +28,10 @@ type triangulator struct { points []Vertex squaredDistances []float64 - ids []int + ids []int32 center Vertex - triangles []int - halfedges []int + triangles []int32 + halfedges []int32 trianglesLen int hull *node hash []*node @@ -73,7 +73,7 @@ return nil } - tri.ids = make([]int, n) + tri.ids = make([]int32, n) // compute bounds x0 := points[0].X @@ -93,10 +93,10 @@ if p.Y > y1 { y1 = p.Y } - tri.ids[i] = i + tri.ids[i] = int32(i) } - var i0, i1, i2 int + var i0, i1, i2 int32 // pick a seed point close to midpoint m := Vertex{X: (x0 + x1) / 2, Y: (y0 + y1) / 2} @@ -104,7 +104,7 @@ for i, p := range points { d := p.SquaredDistance2D(m) if d < minDist { - i0 = i + i0 = int32(i) minDist = d } } @@ -112,12 +112,12 @@ // find point closest to seed point minDist = infinity for i, p := range points { - if i == i0 { + if int32(i) == i0 { continue } d := p.SquaredDistance2D(points[i0]) if d > 0 && d < minDist { - i1 = i + i1 = int32(i) minDist = d } } @@ -125,12 +125,12 @@ // find the third point which forms the smallest circumcircle minRadius := infinity for i, p := range points { - if i == i0 || i == i1 { + if int32(i) == i0 || int32(i) == i1 { continue } r := circumradius(points[i0], points[i1], p) if r < minRadius { - i2 = i + i2 = int32(i) minRadius = r } } @@ -174,8 +174,8 @@ tri.hull = e maxTriangles := 2*n - 5 - tri.triangles = make([]int, maxTriangles*3) - tri.halfedges = make([]int, maxTriangles*3) + tri.triangles = make([]int32, maxTriangles*3) + tri.halfedges = make([]int32, maxTriangles*3) tri.addTriangle(i0, i1, i2, -1, -1, -1) @@ -191,7 +191,7 @@ pp = p // skip seed triangle points - if i == i0 || i == i1 || i == i2 { + if i == int32(i0) || i == int32(i1) || i == int32(i2) { continue } @@ -273,8 +273,8 @@ t.hash[t.hashKey(t.points[e.i])] = e } -func (t *triangulator) addTriangle(i0, i1, i2, a, b, c int) int { - i := t.trianglesLen +func (t *triangulator) addTriangle(i0, i1, i2, a, b, c int32) int32 { + i := int32(t.trianglesLen) t.triangles[i] = i0 t.triangles[i+1] = i1 t.triangles[i+2] = i2 @@ -285,14 +285,14 @@ return i } -func (t *triangulator) link(a, b int) { +func (t *triangulator) link(a, b int32) { t.halfedges[a] = b if b >= 0 { t.halfedges[b] = a } } -func (t *triangulator) legalize(a int) int { +func (t *triangulator) legalize(a int32) int32 { // if the pair of triangles doesn't satisfy the Delaunay condition // (p1 is inside the circumcircle of [p0, pl, pr]), flip them, // then do the same check/flip recursively for the new pair of triangles