# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1537222292 -7200
# Node ID af1d4d44a88ab31639c646c576d6717b821fa6e2
# Parent  79db27e3a999dac5b7911ab08405b2457d107b04
Experimental tin octree indexer.

diff -r 79db27e3a999 -r af1d4d44a88a cmd/tin2octree/main.go
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/main.go	Tue Sep 18 00:11:32 2018 +0200
@@ -0,0 +1,121 @@
+package main
+
+import (
+	"bufio"
+	"database/sql"
+	"flag"
+	"io"
+	"log"
+	"os"
+	"time"
+
+	"github.com/golang/snappy"
+	"github.com/jackc/pgx"
+	"github.com/jackc/pgx/stdlib"
+)
+
+var (
+	bottleneck = flag.String("bottleneck", "", "bottleneck id")
+	date       = flag.String("date", "", "date info")
+	file       = flag.String("file", "", "save to file")
+
+	snap       = flag.Bool("snappy", false, "use snappy compression")
+	dbhost     = flag.String("dbhost", "localhost", "database host")
+	dbport     = flag.Uint("dbport", 5432, "database port")
+	dbname     = flag.String("dbname", "gemma", "database user")
+	dbuser     = flag.String("dbuser", "scott", "database user")
+	dbpassword = flag.String("dbpw", "tiger", "database password")
+	dbssl      = flag.String("dbssl", "prefer", "database SSL mode")
+)
+
+func run(fn func(*sql.DB) error) error {
+
+	// To ease SSL config ride a bit on parsing.
+	cc, err := pgx.ParseConnectionString("sslmode=" + *dbssl)
+	if err != nil {
+		return err
+	}
+
+	// Do the rest manually to allow whitespace in user/password.
+	cc.Host = *dbhost
+	cc.Port = uint16(*dbport)
+	cc.User = *dbuser
+	cc.Password = *dbpassword
+	cc.Database = *dbname
+
+	db := stdlib.OpenDB(cc)
+	defer db.Close()
+
+	return fn(db)
+}
+
+const tinSQL = `
+SELECT ST_AsBinary(ST_DelaunayTriangles(point_cloud::geometry, 0, 2))
+FROM waterway.sounding_results
+WHERE bottleneck_id = $1 AND date_info = $2
+`
+
+func main() {
+	flag.Parse()
+
+	if *bottleneck == "" || *date == "" {
+		log.Fatalln("missing bottleneck or date option.")
+	}
+
+	dateInfo, err := time.Parse("2006-01-02", *date)
+	if err != nil {
+		log.Fatalf("error: %v\n", err)
+	}
+
+	var t tin
+
+	if err := run(func(db *sql.DB) error {
+		start := time.Now()
+		err := db.QueryRow(tinSQL, *bottleneck, dateInfo).Scan(&t)
+		switch {
+		case err == sql.ErrNoRows:
+			return nil
+		case err != nil:
+			return err
+		}
+		log.Printf("query took: %s\n", time.Since(start))
+		return nil
+	}); err != nil {
+		log.Fatalf("error: %v\n", err)
+	}
+
+	tb := &treeBuilder{t: &t}
+	tb.build()
+
+	if *file != "" {
+		f, err := os.Create(*file)
+		if err != nil {
+			log.Printf("error: %v\n", err)
+		}
+		type flusher interface {
+			io.Writer
+			Flush() error
+		}
+		var out flusher
+		if *snap {
+			out = snappy.NewBufferedWriter(f)
+		} else {
+			out = bufio.NewWriter(f)
+		}
+		if err := t.Serialize(out); err != nil {
+			f.Close()
+			log.Fatalf("error: %v\n", err)
+		}
+		if err := tb.Serialize(out); err != nil {
+			f.Close()
+			log.Fatalf("error: %v\n", err)
+		}
+		if err := out.Flush(); err != nil {
+			f.Close()
+			log.Fatalf("error: %v\n", err)
+		}
+		if err := f.Close(); err != nil {
+			log.Fatalf("error: %v\n", err)
+		}
+	}
+}
diff -r 79db27e3a999 -r af1d4d44a88a cmd/tin2octree/tin.go
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/tin.go	Tue Sep 18 00:11:32 2018 +0200
@@ -0,0 +1,234 @@
+package main
+
+import (
+	"bytes"
+	"encoding/binary"
+	"errors"
+	"fmt"
+	"io"
+	"log"
+	"math"
+)
+
+var (
+	errNoByteSlice   = errors.New("Not a byte slice")
+	errTooLessPoints = errors.New("Too less points")
+)
+
+const (
+	wkbXDR       byte   = 0
+	wkbNDR       byte   = 1
+	wkbTinZ      uint32 = 1000 + 16
+	wkbTriangleZ uint32 = 1000 + 17
+)
+
+type tin struct {
+	vertices  []vertex
+	triangles [][]int32
+
+	min vertex
+	max vertex
+}
+
+func (t *tin) FromWKB(data []byte) error {
+	log.Printf("data length %d\n", len(data))
+
+	r := bytes.NewReader(data)
+
+	endian, err := r.ReadByte()
+
+	var order binary.ByteOrder
+
+	switch {
+	case err != nil:
+		return err
+	case endian == wkbNDR:
+		order = binary.LittleEndian
+	case endian == wkbXDR:
+		order = binary.BigEndian
+	default:
+		return fmt.Errorf("unknown byte order %x", endian)
+	}
+
+	var geomType uint32
+	err = binary.Read(r, order, &geomType)
+
+	switch {
+	case err != nil:
+		return err
+	case geomType != wkbTinZ:
+		return fmt.Errorf("unknown geometry type %x", geomType)
+	}
+
+	var num uint32
+	if err = binary.Read(r, order, &num); err != nil {
+		return err
+	}
+
+	vertices := make([]vertex, 0, 100000)
+
+	var v vertex
+
+	v2i := make(map[vertex]int32, 100000)
+
+	var indexPool []int32
+
+	allocIndices := func() []int32 {
+		if len(indexPool) == 0 {
+			indexPool = make([]int32, 3*8*1024)
+		}
+		ids := indexPool[:3]
+		indexPool = indexPool[3:]
+		return ids
+	}
+
+	var triangles [][]int32
+
+	min := vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64}
+	max := vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64}
+
+	for i := uint32(0); i < num; i++ {
+
+		endian, err = r.ReadByte()
+		switch {
+		case err != nil:
+			return err
+		case endian == wkbNDR:
+			order = binary.LittleEndian
+		case endian == wkbXDR:
+			order = binary.BigEndian
+		default:
+			return fmt.Errorf("unknown byte order %x", endian)
+		}
+
+		err = binary.Read(r, order, &geomType)
+		switch {
+		case err != nil:
+			return err
+		case geomType != wkbTriangleZ:
+			return fmt.Errorf("unknown geometry type %d", geomType)
+		}
+
+		var rings uint32
+		if err = binary.Read(r, order, &rings); err != nil {
+			return err
+		}
+		triangle := allocIndices()
+
+		for ring := uint32(0); ring < rings; ring++ {
+			var npoints uint32
+			if err = binary.Read(r, order, &npoints); err != nil {
+				return err
+			}
+
+			if npoints < 3 {
+				return errTooLessPoints
+			}
+
+			for p := uint32(0); p < npoints; p++ {
+				var x, y, z uint64
+				for _, addr := range []*uint64{&x, &y, &z} {
+					if err = binary.Read(r, order, addr); err != nil {
+						return err
+					}
+				}
+				if p >= 3 || ring >= 1 {
+					// Don't store the forth point.
+					continue
+				}
+				// Do this conversion later to spare reflect calls
+				// and allocs in binary.Read.
+				v.x = math.Float64frombits(x)
+				v.y = math.Float64frombits(y)
+				v.z = math.Float64frombits(z)
+				idx, found := v2i[v]
+				if !found {
+					idx = int32(len(vertices))
+					v2i[v] = idx
+					vertices = append(vertices, v)
+					min.minimize(v)
+					max.maximize(v)
+				}
+				triangle[p] = idx
+			}
+		}
+		triangles = append(triangles, triangle)
+	}
+
+	*t = tin{
+		vertices:  vertices,
+		triangles: triangles,
+		min:       min,
+		max:       max,
+	}
+
+	return nil
+}
+
+func (t *tin) Scan(raw interface{}) error {
+
+	data, ok := raw.([]byte)
+	if !ok {
+		return errNoByteSlice
+	}
+	return t.FromWKB(data)
+}
+
+func (v *vertex) write(w io.Writer) error {
+	for _, addr := range []*float64{&v.x, &v.y, &v.z} {
+		if err := binary.Write(
+			w, binary.LittleEndian, math.Float64bits(*addr)); err != nil {
+			return err
+		}
+	}
+	return nil
+}
+
+func (t *tin) Serialize(w io.Writer) error {
+
+	if err := t.min.write(w); err != nil {
+		return err
+	}
+	if err := t.max.write(w); err != nil {
+		return err
+	}
+
+	if err := binary.Write(
+		w, binary.LittleEndian, uint32(len(t.vertices))); err != nil {
+		return err
+	}
+
+	for _, v := range t.vertices {
+		if err := v.write(w); err != nil {
+			return err
+		}
+	}
+	log.Printf("vertices %d (%d)\n", len(t.vertices), len(t.vertices)*3*8)
+
+	if err := binary.Write(
+		w, binary.LittleEndian, uint32(len(t.triangles))); err != nil {
+		return err
+	}
+
+	var buf [binary.MaxVarintLen32]byte
+	var written int
+	var last int32
+	for _, triangle := range t.triangles {
+		for _, idx := range triangle {
+			value := idx - last
+			n := binary.PutVarint(buf[:], int64(value))
+			for p := buf[:n]; len(p) > 0; p = p[n:] {
+				var err error
+				if n, err = w.Write(p); err != nil {
+					return err
+				}
+				written += n
+			}
+			last = idx
+		}
+	}
+	log.Printf("compressed tin indices in bytes: %d (%d)\n",
+		written, 3*4*len(t.triangles))
+
+	return nil
+}
diff -r 79db27e3a999 -r af1d4d44a88a cmd/tin2octree/tree.go
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/tree.go	Tue Sep 18 00:11:32 2018 +0200
@@ -0,0 +1,134 @@
+package main
+
+import (
+	"encoding/binary"
+	"io"
+	"log"
+)
+
+type treeBuilder struct {
+	t      *tin
+	nodes  int
+	leaves int
+	index  []int32
+}
+
+var cubes = [8][2]vertex{
+	{{0.0, 0.0, 0.0}, {0.5, 0.5, 0.5}},
+	{{0.5, 0.0, 0.0}, {1.0, 0.5, 0.5}},
+	{{0.0, 0.0, 0.5}, {0.5, 0.5, 1.0}},
+	{{0.0, 0.5, 0.5}, {0.5, 1.0, 1.0}},
+	{{0.5, 0.0, 0.0}, {1.0, 0.5, 0.5}},
+	{{0.5, 0.5, 0.0}, {1.0, 1.0, 0.5}},
+	{{0.5, 0.0, 0.5}, {1.0, 0.5, 1.0}},
+	{{0.5, 0.5, 0.5}, {1.0, 1.0, 1.0}},
+}
+
+func (tb *treeBuilder) build() {
+
+	triangles := make([]int32, len(tb.t.triangles))
+	for i := range triangles {
+		triangles[i] = int32(i)
+	}
+
+	tb.index = append(tb.index, 0)
+
+	tb.buildRecursive(triangles, tb.t.min, tb.t.max, 0)
+	tb.index[0] = int32(len(tb.index))
+
+	log.Printf("nodes: %d leaves: %d index %d\n",
+		tb.nodes, tb.leaves, tb.index[0])
+}
+
+func (tb *treeBuilder) buildRecursive(
+	triangles []int32,
+	min, max vertex,
+	depth int,
+) int32 {
+	if len(triangles) <= 16 || depth > 8 {
+		pos := len(tb.index)
+		tb.index = append(tb.index, int32(len(tb.index)))
+		tb.index = append(tb.index, triangles...)
+		//log.Printf("leaf entries: %d (%d)\n", len(triangles), depth)
+		tb.leaves++
+		return int32(-(pos + 1))
+	}
+
+	pos := len(tb.index)
+	tb.index = append(tb.index, 0, 0, 0, 0)
+
+	bbox := interpolate(min, max)
+
+	bboxes := make([][2]vertex, len(cubes))
+
+	for i := range cubes {
+		bboxes[i] = [2]vertex{
+			bbox(cubes[i][0]),
+			bbox(cubes[i][1]),
+		}
+	}
+
+	var quandrants [8][]int32
+
+	for _, tri := range triangles {
+		triangle := tb.t.triangles[tri]
+		v0 := tb.t.vertices[triangle[0]]
+		v1 := tb.t.vertices[triangle[1]]
+		v2 := tb.t.vertices[triangle[2]]
+
+		l := v0
+		l.minimize(v1)
+		l.minimize(v2)
+
+		h := v0
+		h.maximize(v1)
+		h.maximize(v2)
+
+		for i := range bboxes {
+			if !(h.less(bboxes[i][0]) || bboxes[i][1].less(l)) {
+				quandrants[i] = append(quandrants[i], tri)
+			}
+		}
+	}
+
+	for i := range quandrants {
+		if len(quandrants[i]) > 0 {
+			child := tb.buildRecursive(
+				quandrants[i],
+				bboxes[i][0], bboxes[i][1],
+				depth+1)
+			tb.index[pos+i] = child
+		}
+	}
+	tb.nodes++
+	return int32(pos)
+}
+
+func (tb *treeBuilder) Serialize(w io.Writer) error {
+	var buf [binary.MaxVarintLen32]byte
+
+	if err := binary.Write(w, binary.LittleEndian, tb.index[0]); err != nil {
+		return err
+	}
+
+	var last int32
+	var written int
+
+	for _, x := range tb.index[1:] {
+		delta := x - last
+		n := binary.PutVarint(buf[:], int64(delta))
+		for p := buf[:n]; len(p) > 0; p = p[n:] {
+			var err error
+			if n, err = w.Write(p); err != nil {
+				return err
+			}
+			written += n
+		}
+
+		last = x
+	}
+	log.Printf("compressed octree index in bytes: %d (%d)\n",
+		written, 4*len(tb.index))
+
+	return nil
+}
diff -r 79db27e3a999 -r af1d4d44a88a cmd/tin2octree/vertex.go
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/vertex.go	Tue Sep 18 00:11:32 2018 +0200
@@ -0,0 +1,107 @@
+package main
+
+import "math"
+
+type vertex struct {
+	x float64
+	y float64
+	z float64
+}
+
+func (v *vertex) minimize(w vertex) {
+	if w.x < v.x {
+		v.x = w.x
+	}
+	if w.y < v.y {
+		v.y = w.y
+	}
+	if w.z < v.z {
+		v.z = w.z
+	}
+}
+
+func (v *vertex) maximize(w vertex) {
+	if w.x > v.x {
+		v.x = w.x
+	}
+	if w.y > v.y {
+		v.y = w.y
+	}
+	if w.z > v.z {
+		v.z = w.z
+	}
+}
+
+type plane struct {
+	a float64
+	b float64
+	c float64
+	d float64
+}
+
+func (a vertex) cross(b vertex) vertex {
+	return vertex{
+		a.y*b.z - a.z*b.y,
+		a.z*b.x - a.x*b.z,
+		a.x*b.y - a.y*b.x,
+	}
+}
+
+func (a vertex) sub(b vertex) vertex {
+	return vertex{
+		a.x - b.x,
+		a.y - b.y,
+		a.z - b.z,
+	}
+}
+
+func (a vertex) scale(s float64) vertex {
+	return vertex{
+		s * a.x,
+		s * a.y,
+		s * a.z,
+	}
+}
+
+func (a vertex) dot(b vertex) float64 {
+	return a.x*b.x + a.y*b.y + a.z*b.z
+}
+
+func (a vertex) length() float64 {
+	return math.Sqrt(a.x*a.x + a.y*a.y + a.z*a.z)
+}
+
+func newPlane(v1, v2, v3 vertex) plane {
+	w1 := v2.sub(v1)
+	w2 := v3.sub(v1)
+	p := w1.cross(w2)
+	p = p.scale(1 / p.length())
+
+	// x*p.x + y*p.y + z*p.z + d = 0
+	// d = -(x*p.x + y*p.y + z*p.z)
+	d := -p.dot(v1)
+	return plane{
+		a: p.x,
+		b: p.y,
+		c: p.z,
+		d: d,
+	}
+}
+
+func (p plane) eval(v vertex) float64 {
+	return p.a*v.x + p.b*v.y + p.c*v.z + p.d
+}
+
+func interpolate(v1, v2 vertex) func(vertex) vertex {
+	return func(s vertex) vertex {
+		return vertex{
+			(v2.x-v1.x)*s.x + v1.x,
+			(v2.y-v1.y)*s.y + v1.y,
+			(v2.z-v1.z)*s.z + v1.z,
+		}
+	}
+}
+
+func (a vertex) less(b vertex) bool {
+	return a.x < b.x || a.y < b.y || a.z < b.z
+}