Mercurial > gemma

changeset 2484:4fa92d468164 octree-diff

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use fogleman's triangulation.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 01 Mar 2019 17:42:46 +0100
parents 620038ade708
children 86173ac7f222
files cmd/octreediff/main.go pkg/octree/node.go pkg/octree/triangulation.go pkg/octree/triangulator.go
diffstat 4 files changed, 108 insertions(+), 112 deletions(-) [+]
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,
 				)
 			}
--- 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 {
--- 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")
 		}
 	}
--- 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