Mercurial > gemma

changeset 679:899116318d58 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 Wed, 19 Sep 2018 17:35:00 +0200
parents 7bb961d750b6 (diff) 3605af94d1ee (current diff)
children c79c7be29a7a
files
diffstat 8 files changed, 897 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/3rdpartylibs.sh	Wed Sep 19 16:37:03 2018 +0200
+++ b/3rdpartylibs.sh	Wed Sep 19 17:35:00 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/loader.go	Wed Sep 19 17:35:00 2018 +0200
@@ -0,0 +1,203 @@
+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
+}
+
+type byteReader interface {
+	io.ByteReader
+	io.Reader
+}
+
+func loadOctreeReader(r byteReader) (*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([]int32)) {
+
+	type frame struct {
+		pos int32
+		min float64
+		max float64
+	}
+
+	if h < ot.min.z || h > ot.max.z {
+		return
+	}
+
+	stack := []frame{{1, ot.min.z, ot.max.z}}
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		pos, min, max := top.pos, top.min, top.max
+		if pos == 0 {
+			continue
+		}
+
+		if pos > 0 { // node
+			if mid := (max-min)*0.5 + min; h <= mid {
+				stack = append(stack,
+					frame{ot.index[pos+0], min, mid},
+					frame{ot.index[pos+1], min, mid},
+					frame{ot.index[pos+4], min, mid},
+					frame{ot.index[pos+5], min, mid})
+			} else {
+				stack = append(stack,
+					frame{ot.index[pos+2], mid, max},
+					frame{ot.index[pos+3], mid, max},
+					frame{ot.index[pos+6], mid, max},
+					frame{ot.index[pos+7], mid, 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 {
+				tri := ot.triangles[idx]
+				v0 := ot.vertices[tri[0]]
+				v1 := ot.vertices[tri[1]]
+				v2 := ot.vertices[tri[2]]
+
+				if !(math.Min(v0.z, math.Min(v1.z, v2.z)) > h ||
+					math.Max(v0.z, math.Max(v1.z, v2.z)) < h) {
+					fn(tri)
+				}
+			}
+		}
+	}
+}
+
+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	Wed Sep 19 17:35:00 2018 +0200
@@ -0,0 +1,40 @@
+package main
+
+import (
+	"flag"
+	"log"
+	"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 poiint")
+)
+
+func process(tree *octree) {
+	if *one {
+		var triangles int
+		tree.horizontal(*step, func([]int32) {
+			triangles++
+		})
+		log.Printf("num triangles: %d\n", triangles)
+	} else {
+	}
+}
+
+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))
+		process(tree)
+	}
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/octree2contour/vertex.go	Wed Sep 19 17:35:00 2018 +0200
@@ -0,0 +1,54 @@
+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 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
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cmd/tin2octree/builder.go	Wed Sep 19 17:35:00 2018 +0200
@@ -0,0 +1,137 @@
+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("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	Wed Sep 19 17:35:00 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	Wed Sep 19 17:35:00 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	Wed Sep 19 17:35:00 2018 +0200
@@ -0,0 +1,54 @@
+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 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
+}