Mercurial > gemma

changeset 704:fe0889460e97 octree

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged default into octree branch.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 21 Sep 2018 09:50:46 +0200
parents 18f6f2d89da1 (diff) e9c28c42c927 (current diff)
children 9db4ae29ded9
files
diffstat 10 files changed, 1186 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/3rdpartylibs.sh	Fri Sep 21 00:24:41 2018 +0200
+++ b/3rdpartylibs.sh	Fri Sep 21 09:50:46 2018 +0200
@@ -8,3 +8,4 @@
 go get -u -v gopkg.in/gomail.v2
 go get -u -v github.com/rs/cors
 go get -u -v golang.org/x/net/html/charset
+go get -u -v github.com/golang/snappy
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/octree2contour/db.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,39 @@
+package main
+
+import (
+	"database/sql"
+	"flag"
+
+	"github.com/jackc/pgx"
+	"github.com/jackc/pgx/stdlib"
+)
+
+var (
+	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)
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/octree2contour/loader.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,205 @@
+package main
+
+import (
+	"bufio"
+	"encoding/binary"
+	"io"
+	"log"
+	"math"
+	"os"
+
+	"github.com/golang/snappy"
+)
+
+type octree struct {
+	epsg uint32
+
+	vertices  []vertex
+	triangles [][]int32
+	index     []int32
+
+	min vertex
+	max vertex
+}
+
+func (v *vertex) read(r io.Reader) error {
+	var buf [8]byte
+	b := buf[:]
+	if _, err := io.ReadFull(r, b); err != nil {
+		return nil
+	}
+	v.x = math.Float64frombits(binary.LittleEndian.Uint64(b))
+	if _, err := io.ReadFull(r, b); err != nil {
+		return nil
+	}
+	v.y = math.Float64frombits(binary.LittleEndian.Uint64(b))
+	if _, err := io.ReadFull(r, b); err != nil {
+		return nil
+	}
+	v.z = math.Float64frombits(binary.LittleEndian.Uint64(b))
+	return nil
+}
+
+func loadOctreeReader(r *bufio.Reader) (*octree, error) {
+	tree := new(octree)
+
+	if err := binary.Read(r, binary.LittleEndian, &tree.epsg); err != nil {
+		return nil, err
+	}
+
+	log.Printf("EPSG: %d\n", tree.epsg)
+
+	if err := tree.min.read(r); err != nil {
+		return nil, err
+	}
+
+	if err := tree.max.read(r); err != nil {
+		return nil, err
+	}
+
+	log.Printf("BBOX: [[%f, %f, %f], [%f, %f, %f]]\n",
+		tree.min.x, tree.min.y, tree.min.z,
+		tree.max.x, tree.max.y, tree.max.z)
+
+	var numVertices uint32
+	if err := binary.Read(r, binary.LittleEndian, &numVertices); err != nil {
+		return nil, err
+	}
+
+	log.Printf("vertices: %d\n", numVertices)
+
+	vertices := make([]vertex, numVertices)
+	tree.vertices = vertices
+
+	for i := range vertices {
+		if err := vertices[i].read(r); err != nil {
+			return nil, err
+		}
+	}
+
+	var numTriangles uint32
+	if err := binary.Read(r, binary.LittleEndian, &numTriangles); err != nil {
+		return nil, err
+	}
+
+	log.Printf("triangles: %d\n", numTriangles)
+
+	indices := make([]int32, 3*numTriangles)
+	triangles := make([][]int32, numTriangles)
+	tree.triangles = triangles
+
+	var last int32
+
+	for i := range triangles {
+		tri := indices[:3]
+		indices = indices[3:]
+		triangles[i] = tri
+		for j := range tri {
+			v, err := binary.ReadVarint(r)
+			if err != nil {
+				return nil, err
+			}
+			value := int32(v) + last
+			tri[j] = value
+			last = value
+		}
+	}
+
+	var numNodes uint32
+	if err := binary.Read(r, binary.LittleEndian, &numNodes); err != nil {
+		return nil, err
+	}
+
+	log.Printf("num nodes: %d\n", numNodes)
+
+	tree.index = make([]int32, numNodes)
+	entries := tree.index[1:]
+
+	last = 0
+	for i := range entries {
+		v, err := binary.ReadVarint(r)
+		if err != nil {
+			return nil, err
+		}
+		value := int32(v) + last
+		entries[i] = value
+		last = value
+	}
+
+	return tree, nil
+}
+
+func (ot *octree) horizontal(h float64, fn func(*triangle)) {
+
+	type frame struct {
+		pos int32
+		min float64
+		max float64
+	}
+
+	if h < ot.min.z || ot.max.z < h {
+		return
+	}
+
+	stack := []frame{{1, ot.min.z, ot.max.z}}
+
+	dupes := map[int32]struct{}{}
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		pos := top.pos
+		if pos == 0 {
+			continue
+		}
+		min, max := top.min, top.max
+
+		if pos > 0 { // node
+			if mid := (max-min)*0.5 + min; h >= mid {
+				pos += 4 // nodes with z-bit set
+				min = mid
+			} else {
+				max = mid
+			}
+			stack = append(stack,
+				frame{ot.index[pos+0], min, max},
+				frame{ot.index[pos+1], min, max},
+				frame{ot.index[pos+2], min, max},
+				frame{ot.index[pos+3], min, max})
+		} else { // leaf
+			pos = -pos - 1
+			n := ot.index[pos]
+			//log.Printf("%d %d %d\n", pos, n, len(ot.index))
+			indices := ot.index[pos+1 : pos+1+n]
+
+			for _, idx := range indices {
+				if _, found := dupes[idx]; found {
+					continue
+				}
+				tri := ot.triangles[idx]
+				t := triangle{
+					ot.vertices[tri[0]],
+					ot.vertices[tri[1]],
+					ot.vertices[tri[2]],
+				}
+
+				if !(math.Min(t[0].z, math.Min(t[1].z, t[2].z)) > h ||
+					math.Max(t[0].z, math.Max(t[1].z, t[2].z)) < h) {
+					dupes[idx] = struct{}{}
+					fn(&t)
+				}
+			}
+		}
+	}
+}
+
+func loadOctree(fname string) (*octree, error) {
+
+	f, err := os.Open(fname)
+	if err != nil {
+		return nil, err
+	}
+	defer f.Close()
+	return loadOctreeReader(bufio.NewReader(snappy.NewReader(f)))
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/octree2contour/main.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,106 @@
+package main
+
+import (
+	"flag"
+	"log"
+	"runtime"
+	"sort"
+	"sync"
+	"time"
+)
+
+var (
+	one  = flag.Bool("o", false, "create only a single contour")
+	step = flag.Float64("s", 0.5, "step with")
+	max  = flag.Float64("m", 10, "max height from lowest point")
+)
+
+func processLevels(
+	tree *octree,
+	jobs <-chan float64,
+	results chan<- result,
+	wg *sync.WaitGroup,
+) {
+	defer wg.Done()
+	for h := range jobs {
+		var lines multiLineStringZ
+		tree.horizontal(h, func(t *triangle) {
+			line := t.intersectH(h)
+			if len(line) > 1 {
+				lines = append(lines, line)
+			}
+		})
+		results <- result{h, lines}
+	}
+}
+
+func process(tree *octree) []result {
+
+	if *one {
+		var lines multiLineStringZ
+		tree.horizontal(*step, func(t *triangle) {
+			line := t.intersectH(*step)
+			if len(line) > 0 {
+				lines = append(lines, line)
+			}
+		})
+		return []result{{*step, lines}}
+	}
+
+	results := make(chan result)
+	done := make(chan struct{})
+
+	var all []result
+	go func() {
+		for {
+			select {
+			case r := <-results:
+				all = append(all, r)
+			case <-done:
+				return
+			}
+		}
+	}()
+
+	var wg sync.WaitGroup
+	jobs := make(chan float64)
+	for i, n := 0, runtime.NumCPU(); i < n; i++ {
+		wg.Add(1)
+		go processLevels(tree, jobs, results, &wg)
+	}
+	for h := tree.min.z; h <= tree.max.z; h += *step {
+		jobs <- h
+	}
+	close(jobs)
+	wg.Wait()
+	done <- struct{}{}
+
+	sort.Slice(all, func(i, j int) bool {
+		return all[i].h < all[j].h
+	})
+
+	return all
+}
+
+func main() {
+	flag.Parse()
+
+	for _, fname := range flag.Args() {
+		log.Printf("processing %s\n", fname)
+		start := time.Now()
+		tree, err := loadOctree(fname)
+		if err != nil {
+			log.Printf("error: %v\n", err)
+			continue
+		}
+		log.Printf("loading took: %v\n", time.Since(start))
+		start = time.Now()
+		all := process(tree)
+		log.Printf("processing took: %v\n", time.Since(start))
+		start = time.Now()
+		if err = store(all, tree.epsg); err != nil {
+			log.Printf("error: %v\n", err)
+		}
+		log.Printf("storing took: %v\n", time.Since(start))
+	}
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/octree2contour/store.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,75 @@
+package main
+
+import (
+	"bytes"
+	"database/sql"
+	"encoding/binary"
+	"math"
+)
+
+type multiLineStringZ [][]vertex
+
+const (
+	wkbNDR              byte   = 1
+	wkbPointZ           uint32 = 1000 + 1
+	wkbLineStringZ      uint32 = 1000 + 2
+	wkbMultiLineStringZ uint32 = 1000 + 5
+)
+
+type result struct {
+	h     float64
+	lines multiLineStringZ
+}
+
+const insertSQL = `
+INSERT INTO waterway.contour_lines (height, geom)
+VALUES ($1, ST_Transform(
+	ST_SetSRID(ST_GeomFromWKB($2), $3), 4326)::geography)
+`
+
+func store(all []result, epsg uint32) error {
+
+	return run(func(db *sql.DB) error {
+
+		tx, err := db.Begin()
+		if err != nil {
+			return err
+		}
+		defer tx.Rollback()
+
+		stmt, err := tx.Prepare(insertSQL)
+		if err != nil {
+			return err
+		}
+
+		for _, r := range all {
+			if _, err := stmt.Exec(r.h, r.lines.asWKB(), epsg); err != nil {
+				return err
+			}
+		}
+
+		return tx.Commit()
+	})
+}
+
+func (mls multiLineStringZ) asWKB() []byte {
+
+	var buf bytes.Buffer
+
+	binary.Write(&buf, binary.LittleEndian, wkbNDR)
+	binary.Write(&buf, binary.LittleEndian, wkbMultiLineStringZ)
+	binary.Write(&buf, binary.LittleEndian, uint32(len(mls)))
+
+	for _, ml := range mls {
+		binary.Write(&buf, binary.LittleEndian, wkbNDR)
+		binary.Write(&buf, binary.LittleEndian, wkbLineStringZ)
+		binary.Write(&buf, binary.LittleEndian, uint32(len(ml)))
+		for _, p := range ml {
+			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(p.z))
+		}
+	}
+
+	return buf.Bytes()
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/octree2contour/vertex.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,136 @@
+package main
+
+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
+	}
+}
+
+func (v vertex) sub(w vertex) vertex {
+	return vertex{
+		v.x - w.x,
+		v.y - w.y,
+		v.z - w.z,
+	}
+}
+
+func (v vertex) add(w vertex) vertex {
+	return vertex{
+		v.x + w.x,
+		v.y + w.y,
+		v.z + w.z,
+	}
+}
+
+func (v vertex) scale(s float64) vertex {
+	return vertex{
+		s * v.x,
+		s * v.y,
+		s * v.z,
+	}
+}
+
+func interpolate(v1, v2 vertex) func(vertex) vertex {
+	v2 = v2.sub(v1)
+	return func(s vertex) vertex {
+		return vertex{
+			v2.x*s.x + v1.x,
+			v2.y*s.y + v1.y,
+			v2.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
+}
+
+type line [2]vertex
+
+func newLine(p1, p2 vertex) line {
+	return line{
+		p2.sub(p1),
+		p1,
+	}
+}
+
+func (l line) eval(t float64) vertex {
+	return l[0].scale(t).add(l[1])
+}
+
+func (l line) intersectH(h float64) vertex {
+	t := (h - l[1].z) / l[0].z
+	return l.eval(t)
+}
+
+type triangle [3]vertex
+
+func side(z, h float64) int {
+	switch {
+	case z < h:
+		return -1
+	case z > h:
+		return +1
+	}
+	return 0
+}
+
+func (t *triangle) intersectH(h float64) []vertex {
+	sides := [3]int{
+		side(t[0].z, h),
+		side(t[1].z, h),
+		side(t[2].z, h),
+	}
+
+	var points []vertex
+
+	for i := 0; i < 3; i++ {
+		j := (i + 1) % 3
+		si := sides[i]
+		sj := sides[j]
+
+		switch {
+		case si == 0 && sj == 0:
+			// both on plane
+			points = append(points, t[i], t[j])
+		case si == 0 && sj != 0:
+			// first on plane
+			points = append(points, t[i])
+		case si != 0 && sj == 0:
+			// second on plane
+			points = append(points, t[j])
+		case si == sj:
+			// both on same side
+		default:
+			// real intersection
+			v := newLine(t[i], t[j]).intersectH(h)
+			points = append(points, v)
+		}
+	}
+
+	return points
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/builder.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,154 @@
+package main
+
+import (
+	"encoding/binary"
+	"io"
+	"log"
+)
+
+type treeBuilder struct {
+	t      *tin
+	nodes  int
+	leaves int
+	index  []int32
+}
+
+var cubes = [8][2]vertex{
+	makeCube(0),
+	makeCube(1),
+	makeCube(2),
+	makeCube(3),
+	makeCube(4),
+	makeCube(5),
+	makeCube(6),
+	makeCube(7),
+}
+
+func makeCube(i int) [2]vertex {
+	var d vertex
+	if i&1 == 1 {
+		d.x = 0.5
+	}
+	if i&2 == 2 {
+		d.y = 0.5
+	}
+	if i&4 == 4 {
+		d.z = 0.5
+	}
+	return [2]vertex{
+		vertex{0.0, 0.0, 0.0}.add(d),
+		vertex{0.5, 0.5, 0.5}.add(d),
+	}
+}
+
+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("num nodes: %d\n", tb.index[0])
+
+	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(triangles)))
+		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,
+		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
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/main.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,163 @@
+package main
+
+import (
+	"database/sql"
+	"flag"
+	"log"
+	"math"
+	"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")
+
+	utm = flag.Bool("utm", false, "fetch in matchin UTM zone")
+
+	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
+`
+	tinUTMSQL = `
+SELECT ST_AsBinary(
+  ST_DelaunayTriangles(
+	  ST_Transform(point_cloud::geometry, $3::int), 0, 2))
+FROM waterway.sounding_results
+WHERE bottleneck_id = $1 AND date_info = $2
+`
+	centroidSQL = `
+SELECT ST_X(ST_Centroid(point_cloud::geometry)), ST_Y(ST_Centroid(point_cloud::geometry))
+FROM waterway.sounding_results
+WHERE bottleneck_id = $1 AND date_info = $2
+`
+)
+
+func utmZone(x, y float64) int {
+	var pref int
+	if y > 0 {
+		pref = 32600
+	} else {
+		pref = 32700
+	}
+	zone := int(math.Floor((x+180)/6)) + 1
+	return zone + pref
+}
+
+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 {
+		var utmZ int
+
+		if *utm {
+			var cx, cy float64
+			err := db.QueryRow(centroidSQL, *bottleneck, dateInfo).Scan(&cx, &cy)
+			switch {
+			case err == sql.ErrNoRows:
+				return nil
+			case err != nil:
+				return err
+			}
+			log.Printf("lat/lon: [%f, %f]\n", cx, cy)
+			utmZ = utmZone(cx, cy)
+			log.Printf("UTM zone: %d\n", utmZ)
+		}
+
+		start := time.Now()
+		var err error
+		if *utm {
+			err = db.QueryRow(tinUTMSQL, *bottleneck, dateInfo, utmZ).Scan(&t)
+		} else {
+			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))
+
+		if *utm {
+			t.epsg = uint32(utmZ)
+		}
+
+		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)
+		}
+		out := snappy.NewBufferedWriter(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)
+		}
+	}
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/tin.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,245 @@
+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
+)
+
+const wgs84 = 4326
+
+type tin struct {
+	epsg      uint32
+	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)
+	}
+
+	log.Printf("bbox: [[%f, %f], [%f, %f]]\n",
+		min.x, min.y, max.x, max.y)
+
+	*t = tin{
+		epsg:      wgs84,
+		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 := binary.Write(w, binary.LittleEndian, t.epsg); err != nil {
+		return err
+	}
+
+	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
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/vertex.go	Fri Sep 21 09:50:46 2018 +0200
@@ -0,0 +1,62 @@
+package main
+
+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
+	}
+}
+
+func (v vertex) sub(w vertex) vertex {
+	return vertex{
+		v.x - w.x,
+		v.y - w.y,
+		v.z - w.z,
+	}
+}
+
+func (v vertex) add(w vertex) vertex {
+	return vertex{
+		v.x + w.x,
+		v.y + w.y,
+		v.z + w.z,
+	}
+}
+
+func interpolate(v1, v2 vertex) func(vertex) vertex {
+	v2 = v2.sub(v1)
+	return func(s vertex) vertex {
+		return vertex{
+			v2.x*s.x + v1.x,
+			v2.y*s.y + v1.y,
+			v2.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
+}