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,
 				)
 			}