Mercurial > gemma

changeset 4825:5eb05714353a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged remove-octree-debris into default. This removes the oct2str tool so db migrations are only possible after or equal version 1303.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 05 Nov 2019 13:01:08 +0100
parents 8bca7b744459 (current diff) c0eb491aaaa7 (diff)
children 39ee00d06309
files cmd/oct2str/main.go pkg/octree/builder.go pkg/octree/contours.go pkg/octree/tree.go
diffstat 8 files changed, 2 insertions(+), 1050 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile	Mon Nov 04 18:01:34 2019 +0100
+++ b/Makefile	Tue Nov 05 13:01:08 2019 +0100
@@ -14,7 +14,7 @@
 
 .PHONY: all gemma client clean
 
-all: gemma client tools
+all: gemma client
 
 $(ENVWARPPER):
 	@echo "Preparing go build environment:"
@@ -27,14 +27,9 @@
 	echo 'exec "$$@"' >>"$(ENVWARPPER)"
 	chmod +x "$(ENVWARPPER)"
 
-tools: oct2str
-
 gemma: $(ENVWARPPER)
 	"$(ENVWARPPER)" go build -o ./cmd/gemma/gemma ./cmd/gemma
 
-oct2str: $(ENVWARPPER)
-	"$(ENVWARPPER)" go build -o ./cmd/oct2str/oct2str ./cmd/oct2str
-
 client:
 	$(MAKE) -f Makefile.build -C client
 
@@ -45,7 +40,7 @@
 	v="gemma-$$(hg id -i)" ;\
         tar --transform "s@^@$${v}/@" \
 	    -cJf "../$${v}.tar.xz" \
-	    cmd/gemma/gemma cmd/oct2str/oct2str schema style-templates \
+	    schema style-templates \
 	    web misc example_conf.toml
 
 clean:
--- a/cmd/oct2str/main.go	Mon Nov 04 18:01:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,167 +0,0 @@
-// This is Free Software under GNU Affero General Public License v >= 3.0
-// without warranty, see README.md and license for details.
-//
-// SPDX-License-Identifier: AGPL-3.0-or-later
-// License-Filename: LICENSES/AGPL-3.0.txt
-//
-// Copyright (C) 2019 by via donau
-//   – Österreichische Wasserstraßen-Gesellschaft mbH
-// Software engineering by Intevation GmbH
-//
-// Author(s):
-//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
-
-package main
-
-import (
-	"context"
-	"crypto/sha1"
-	"database/sql"
-	"encoding/hex"
-	"flag"
-	"log"
-
-	"gemma.intevation.de/gemma/pkg/octree"
-	"github.com/jackc/pgx"
-	"github.com/jackc/pgx/stdlib"
-)
-
-const (
-	fetchOneOctreeSQL = `
-SELECT
-  id,
-  octree_index
-FROM waterway.sounding_results
-WHERE mesh_index IS NULL
-LIMIT 1`
-
-	storeSTRTreeSQL = `
-UPDATE waterway.sounding_results
-SET mesh_index = $1, mesh_checksum = $2
-WHERE id = $3`
-)
-
-func process(ctx context.Context, conn *sql.Conn, count int) error {
-
-	fetch, err := conn.PrepareContext(ctx, fetchOneOctreeSQL)
-	if err != nil {
-		return err
-	}
-	defer fetch.Close()
-
-	store, err := conn.PrepareContext(ctx, storeSTRTreeSQL)
-	if err != nil {
-		return err
-	}
-	defer store.Close()
-
-	var next func() bool
-	if count < 0 {
-		next = func() bool { return true }
-	} else {
-		var c int
-		next = func() bool {
-			if c < count {
-				c++
-				return true
-			}
-			return false
-		}
-	}
-
-loop:
-	for next() {
-		switch err := func() error {
-			tx, err := conn.BeginTx(ctx, nil)
-			if err != nil {
-				return err
-			}
-			defer tx.Rollback()
-
-			var id int64
-			var data []byte
-
-			if err = tx.Stmt(fetch).QueryRowContext(ctx).Scan(&id, &data); err != nil {
-				return err
-			}
-
-			otree, err := octree.Deserialize(data)
-			if err != nil {
-				return err
-			}
-
-			unused := otree.FindUnused()
-
-			log.Printf("unused: %d\n", len(unused))
-
-			str := octree.STRTree{Entries: 16}
-			str.BuildWithout(otree.Tin(), unused)
-
-			out, err := str.Bytes()
-			if err != nil {
-				return err
-			}
-			h := sha1.New()
-			h.Write(out)
-			checksum := hex.EncodeToString(h.Sum(nil))
-			log.Printf("hash: %s\n", checksum)
-
-			if _, err := tx.Stmt(store).ExecContext(ctx, out, checksum, id); err != nil {
-				return err
-			}
-
-			return tx.Commit()
-		}(); {
-		case err == sql.ErrNoRows:
-			break loop
-		case err != nil:
-			return err
-		}
-	}
-
-	return nil
-}
-
-func connect(cc pgx.ConnConfig, count int) error {
-	db := stdlib.OpenDB(cc)
-	defer db.Close()
-
-	ctx := context.Background()
-
-	conn, err := db.Conn(ctx)
-	if err != nil {
-		return err
-	}
-	defer conn.Close()
-
-	return process(ctx, conn, count)
-}
-
-func main() {
-	var (
-		db       = flag.String("database", "gemma", "database name")
-		user     = flag.String("user", "sophie", "database user")
-		host     = flag.String("host", "localhost", "database host")
-		password = flag.String("password", "so2Phie4", "database password")
-		port     = flag.Uint("port", 5432, "database port")
-		ssl      = flag.String("ssl", "prefer", "SSL mode")
-		count    = flag.Int("count", -1, "how many sounding results to convert")
-	)
-
-	flag.Parse()
-	cc, err := pgx.ParseConnectionString("sslmode=" + *ssl)
-	if err != nil {
-		log.Fatalf("error: %v\n", err)
-	}
-
-	// Do the rest manually to allow whitespace in user/password.
-	cc.Host = *host
-	cc.Port = uint16(*port)
-	cc.User = *user
-	cc.Password = *password
-	cc.Database = *db
-
-	if err := connect(cc, *count); err != nil {
-		log.Fatalf("error: %v\n", err)
-	}
-}
--- a/go.mod	Mon Nov 04 18:01:34 2019 +0100
+++ b/go.mod	Tue Nov 05 13:01:08 2019 +0100
@@ -7,7 +7,6 @@
 	github.com/etcd-io/bbolt v1.3.3
 	github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199
 	github.com/gofrs/uuid v3.2.0+incompatible // indirect
-	github.com/golang/snappy v0.0.1
 	github.com/gorilla/mux v1.7.3
 	github.com/jackc/fake v0.0.0-20150926172116-812a484cc733 // indirect
 	github.com/jackc/pgx v3.6.0+incompatible
--- a/go.sum	Mon Nov 04 18:01:34 2019 +0100
+++ b/go.sum	Tue Nov 05 13:01:08 2019 +0100
@@ -46,8 +46,6 @@
 github.com/golang/mock v1.1.1/go.mod h1:oTYuIxOrZwtPieC+H1uAHpcLFnEyAGVDL/k47Jfbm0A=
 github.com/golang/protobuf v1.2.0/go.mod h1:6lQm79b+lXiMfvg/cZm0SGofjICqVBUtrP5yJMmIC1U=
 github.com/golang/protobuf v1.3.1/go.mod h1:6lQm79b+lXiMfvg/cZm0SGofjICqVBUtrP5yJMmIC1U=
-github.com/golang/snappy v0.0.1 h1:Qgr9rKW7uDUkrbSmQeiDsGa8SjGyCOGtuasMWwvp2P4=
-github.com/golang/snappy v0.0.1/go.mod h1:/XxbfmMg8lxefKM7IXC3fBNl/7bRcc72aCRzEWrmP2Q=
 github.com/google/btree v1.0.0/go.mod h1:lNA+9X1NB3Zf8V7Ke586lFgjr2dZNuvo3lPJSGZ5JPQ=
 github.com/google/go-cmp v0.2.0/go.mod h1:oXzfMopK8JAjlY9xF4vHSVASa0yLyX7SntLO5aqRK0M=
 github.com/gorilla/mux v1.7.3 h1:gnP5JzjVOuiZD07fKKToCAOjS0yOpj/qPETTXCCS6hw=
--- a/pkg/octree/builder.go	Mon Nov 04 18:01:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,389 +0,0 @@
-// This is Free Software under GNU Affero General Public License v >= 3.0
-// without warranty, see README.md and license for details.
-//
-// SPDX-License-Identifier: AGPL-3.0-or-later
-// License-Filename: LICENSES/AGPL-3.0.txt
-//
-// Copyright (C) 2018 by via donau
-//   – Österreichische Wasserstraßen-Gesellschaft mbH
-// Software engineering by Intevation GmbH
-//
-// Author(s):
-//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
-
-package octree
-
-import (
-	"bytes"
-	"encoding/binary"
-	"io"
-	"log"
-	"runtime"
-	"sync"
-	"sync/atomic"
-
-	"github.com/golang/snappy"
-)
-
-// Builder is used to turn a TIN into an Octree.
-type Builder struct {
-	t      *Tin
-	nodes  int
-	leaves int
-	index  []int32
-
-	mu sync.Mutex
-}
-
-type buildStep func(chan buildStep)
-
-var cubes = [8]Box{
-	makeCube(0),
-	makeCube(1),
-	makeCube(2),
-	makeCube(3),
-	makeCube(4),
-	makeCube(5),
-	makeCube(6),
-	makeCube(7),
-}
-
-func makeCube(i int) Box {
-	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 Box{
-		Vertex{0.0, 0.0, 0.0}.Add(d),
-		Vertex{0.5, 0.5, 0.5}.Add(d),
-	}
-}
-
-func twoElseOne(b bool) int {
-	if b {
-		return 2
-	}
-	return 1
-}
-
-// NewBuilder creates a new Builder for a TIN.
-func NewBuilder(t *Tin) *Builder {
-	return &Builder{t: t}
-}
-
-// Build builds the Octree.
-func (tb *Builder) Build(removed map[int32]struct{}) {
-
-	var triangles []int32
-
-	if len(removed) > 0 {
-		triangles = make([]int32, len(tb.t.Triangles)-len(removed))
-		idx := 0
-		for i := range tb.t.Triangles {
-			if _, found := removed[int32(i)]; !found {
-				triangles[idx] = int32(i)
-				idx++
-			}
-		}
-	} else {
-		triangles = make([]int32, len(tb.t.Triangles))
-		for i := range triangles {
-			triangles[i] = int32(i)
-		}
-	}
-
-	n := runtime.NumCPU()
-
-	steps := make(chan buildStep)
-
-	var wg sync.WaitGroup
-	for i := 0; i < n; i++ {
-		wg.Add(1)
-		go func() {
-			defer wg.Done()
-			for step := range steps {
-				step(steps)
-			}
-		}()
-	}
-
-	tb.index = append(tb.index, 0)
-
-	root := func(int32) {
-		close(steps)
-	}
-
-	steps <- tb.buildConcurrent(
-		triangles,
-		tb.t.Min, tb.t.Max,
-		0,
-		root)
-
-	wg.Wait()
-
-	/*
-		tb.buildRecursive(triangles, tb.t.Min, tb.t.Max, 0)
-	*/
-	tb.index[0] = int32(len(tb.index))
-	log.Printf("info: num nodes: %d\n", tb.index[0])
-	log.Printf("info: nodes: %d leaves: %d index %d\n",
-		tb.nodes, tb.leaves, tb.index[0])
-}
-
-func (tb *Builder) buildConcurrent(
-	triangles []int32,
-	min, max Vertex,
-	depth int,
-	parent func(int32),
-) buildStep {
-
-	return func(steps chan buildStep) {
-
-		// none concurrent for small parts.
-		if len(triangles) <= 1024 || depth > 8 {
-			parent(tb.buildRecursive(triangles, min, max, depth))
-			return
-		}
-		box := Box{min, max}
-
-		xLimit := twoElseOne(box.HasX())
-		yLimit := twoElseOne(box.HasY())
-		zLimit := twoElseOne(box.HasZ())
-
-		indices := make([]byte, 0, 8)
-
-		bbox := box.Interpolate()
-
-		var bboxes [8]Box
-
-		for x := 0; x < xLimit; x++ {
-			for y := 0; y < yLimit; y++ {
-				for z := 0; z < zLimit; z++ {
-					idx := byte(z<<2 | y<<1 | x)
-					bboxes[idx] = Box{
-						bbox(cubes[idx][0]),
-						bbox(cubes[idx][1]),
-					}
-					indices = append(indices, idx)
-				}
-			}
-		}
-
-		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 indices {
-				if !(h.Less(bboxes[i][0]) || bboxes[i][1].Less(l)) {
-					quandrants[i] = append(quandrants[i], tri)
-				}
-			}
-		}
-
-		used := new(int32)
-		for _, i := range indices {
-			if len(quandrants[i]) > 0 {
-				*used++
-			}
-		}
-
-		pos := tb.allocNode()
-
-		if *used == 0 {
-			parent(pos)
-			return
-		}
-
-		for _, i := range indices {
-			if len(quandrants[i]) > 0 {
-				j := int32(i)
-				parent := func(v int32) {
-					tb.index[pos+j] = v
-					if atomic.AddInt32(used, -1) == 0 {
-						parent(pos)
-					}
-				}
-				step := tb.buildConcurrent(
-					quandrants[i],
-					bboxes[i][0], bboxes[i][1],
-					depth+1,
-					parent)
-				select {
-				case steps <- step:
-				default:
-					// all slots busy -> execute directly.
-					step(steps)
-				}
-			}
-		}
-	}
-}
-
-func (tb *Builder) allocNode() int32 {
-	tb.mu.Lock()
-	pos := int32(len(tb.index))
-	tb.index = append(tb.index,
-		0, 0, 0, 0,
-		0, 0, 0, 0)
-	tb.nodes++
-	tb.mu.Unlock()
-	return pos
-}
-
-func (tb *Builder) buildRecursive(
-	triangles []int32,
-	min, max Vertex,
-	depth int,
-) int32 {
-	if len(triangles) <= 16 || depth > 8 {
-		tb.mu.Lock()
-		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++
-		tb.mu.Unlock()
-		return int32(-(pos + 1))
-	}
-
-	box := Box{min, max}
-
-	xLimit := twoElseOne(box.HasX())
-	yLimit := twoElseOne(box.HasY())
-	zLimit := twoElseOne(box.HasZ())
-
-	indices := make([]byte, 0, 8)
-
-	bbox := box.Interpolate()
-
-	var bboxes [8]Box
-
-	for x := 0; x < xLimit; x++ {
-		for y := 0; y < yLimit; y++ {
-			for z := 0; z < zLimit; z++ {
-				idx := byte(z<<2 | y<<1 | x)
-				bboxes[idx] = Box{
-					bbox(cubes[idx][0]),
-					bbox(cubes[idx][1]),
-				}
-				indices = append(indices, idx)
-			}
-		}
-	}
-
-	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 indices {
-			if !(h.Less(bboxes[i][0]) || bboxes[i][1].Less(l)) {
-				quandrants[i] = append(quandrants[i], tri)
-			}
-		}
-	}
-
-	pos := tb.allocNode()
-
-	for _, i := range indices {
-		if len(quandrants[i]) > 0 {
-			child := tb.buildRecursive(
-				quandrants[i],
-				bboxes[i][0], bboxes[i][1],
-				depth+1)
-			tb.index[pos+int32(i)] = child
-		}
-	}
-
-	return pos
-}
-
-func (tb *Builder) 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("info: compressed octree index in bytes: %d (%d)\n",
-		written, 4*len(tb.index))
-
-	return nil
-}
-
-func (tb *Builder) writeTo(w io.Writer) error {
-	out := snappy.NewBufferedWriter(w)
-	if err := tb.t.serialize(out); err != nil {
-		return err
-	}
-	if err := tb.serialize(out); err != nil {
-		return err
-	}
-	return out.Flush()
-}
-
-// Bytes serializes an Octree into a byte slice.
-func (tb *Builder) Bytes() ([]byte, error) {
-	var buf bytes.Buffer
-	if err := tb.writeTo(&buf); err != nil {
-		return nil, err
-	}
-	return buf.Bytes(), nil
-}
-
-// Tree returns an Octree from the Builder.
-func (tb *Builder) Tree() *Tree {
-	return &Tree{
-		EPSG:      tb.t.EPSG,
-		vertices:  tb.t.Vertices,
-		triangles: tb.t.Triangles,
-		index:     tb.index,
-		Min:       tb.t.Min,
-		Max:       tb.t.Max,
-	}
-}
--- a/pkg/octree/contours.go	Mon Nov 04 18:01:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,122 +0,0 @@
-// This is Free Software under GNU Affero General Public License v >= 3.0
-// without warranty, see README.md and license for details.
-//
-// SPDX-License-Identifier: AGPL-3.0-or-later
-// License-Filename: LICENSES/AGPL-3.0.txt
-//
-// Copyright (C) 2018 by via donau
-//   – Österreichische Wasserstraßen-Gesellschaft mbH
-// Software engineering by Intevation GmbH
-//
-// Author(s):
-//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
-//  * Tom Gottfried <tom.gottfried@intevation.de>
-
-package octree
-
-import (
-	"runtime"
-	"sync"
-)
-
-// ContourResult stores an calculated iso line for a given height.
-// Is used as a future variable in the concurrent iso line calculation.
-type ContourResult struct {
-	Height float64
-	Lines  MultiLineStringZ
-
-	done bool
-	mu   sync.Mutex
-	cond *sync.Cond
-}
-
-// NewContourResult prepares a future variable to later hold
-// the result of the iso line calculation.
-func NewContourResult(height float64) *ContourResult {
-	cr := ContourResult{Height: height}
-	cr.cond = sync.NewCond(&cr.mu)
-	return &cr
-}
-
-func (cr *ContourResult) wait() {
-	cr.cond.L.Lock()
-	for !cr.done {
-		cr.cond.Wait()
-	}
-	cr.cond.L.Unlock()
-}
-
-func (cr *ContourResult) get() float64 {
-	cr.cond.L.Lock()
-	defer cr.cond.L.Unlock()
-	return cr.Height
-}
-
-func (cr *ContourResult) set(lines MultiLineStringZ) {
-	cr.cond.L.Lock()
-	defer cr.cond.L.Unlock()
-	cr.Lines = lines
-	cr.done = true
-	cr.cond.Signal()
-}
-
-// DoContours calculates the iso line for the given heights.
-// This is done concurrently.
-// It is guaranteed that the results are given to the store
-// function in order of the original heights values.
-func DoContours(tree *Tree, heights []float64, store func(*ContourResult)) {
-
-	if len(heights) == 0 {
-		return
-	}
-
-	contours := make([]*ContourResult, len(heights))
-
-	for i, h := range heights {
-		contours[i] = NewContourResult(h)
-	}
-
-	jobs := make(chan *ContourResult)
-
-	var wg sync.WaitGroup
-	for i, n := 0, runtime.NumCPU(); i < n; i++ {
-		wg.Add(1)
-		go processLevels(tree, jobs, &wg)
-	}
-
-	done := make(chan struct{})
-	go func() {
-		defer close(done)
-		for _, cr := range contours {
-			cr.wait()
-			store(cr)
-		}
-	}()
-
-	for _, cr := range contours {
-		jobs <- cr
-	}
-	close(jobs)
-
-	wg.Wait()
-	<-done
-}
-
-func processLevels(
-	tree *Tree,
-	jobs <-chan *ContourResult,
-	wg *sync.WaitGroup,
-) {
-	defer wg.Done()
-	for cr := range jobs {
-		var lines MultiLineStringZ
-		h := cr.get()
-		tree.Horizontal(h, func(t *Triangle) {
-			line := t.IntersectHorizontal(h)
-			if len(line) > 1 {
-				lines = append(lines, line)
-			}
-		})
-		cr.set(lines.Merge())
-	}
-}
--- a/pkg/octree/loader.go	Mon Nov 04 18:01:34 2019 +0100
+++ b/pkg/octree/loader.go	Tue Nov 05 13:01:08 2019 +0100
@@ -19,8 +19,6 @@
 	"compress/gzip"
 	"encoding/binary"
 	"log"
-
-	"github.com/golang/snappy"
 )
 
 func (s *STRTree) deserializeIndex(r *bufio.Reader) error {
@@ -166,49 +164,3 @@
 
 	return nil
 }
-
-func loadReader(r *bufio.Reader) (*Tree, error) {
-	tree := new(Tree)
-
-	var tin Tin
-
-	if err := tin.Deserialize(r); err != nil {
-		return nil, err
-	}
-
-	tree.EPSG = tin.EPSG
-	tree.vertices = tin.Vertices
-	tree.triangles = tin.Triangles
-	tree.Min = tin.Min
-	tree.Max = tin.Max
-
-	var numNodes uint32
-	if err := binary.Read(r, binary.LittleEndian, &numNodes); err != nil {
-		return nil, err
-	}
-
-	log.Printf("info: num nodes: %d\n", numNodes)
-
-	tree.index = make([]int32, numNodes)
-	entries := tree.index[1:]
-
-	var last int32
-	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 Deserialize(data []byte) (*Tree, error) {
-	return loadReader(
-		bufio.NewReader(
-			snappy.NewReader(
-				bytes.NewReader(data))))
-}
--- a/pkg/octree/tree.go	Mon Nov 04 18:01:34 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,314 +0,0 @@
-// This is Free Software under GNU Affero General Public License v >= 3.0
-// without warranty, see README.md and license for details.
-//
-// SPDX-License-Identifier: AGPL-3.0-or-later
-// License-Filename: LICENSES/AGPL-3.0.txt
-//
-// Copyright (C) 2018 by via donau
-//   – Österreichische Wasserstraßen-Gesellschaft mbH
-// Software engineering by Intevation GmbH
-//
-// Author(s):
-//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
-
-package octree
-
-import (
-	"math"
-)
-
-// Tree is an Octree holding triangles.
-type Tree struct {
-	// EPSG is the projection.
-	EPSG uint32
-
-	vertices  []Vertex
-	triangles [][]int32
-	index     []int32
-
-	// Min is the lower left corner of the bbox.
-	Min Vertex
-	// Max is the upper right corner of the bbox.
-	Max Vertex
-}
-
-type boxFrame struct {
-	pos int32
-	Box2D
-}
-
-func (ot *Tree) Vertices() []Vertex {
-	return ot.vertices
-}
-
-var scale = [4][4]float64{
-	{0.0, 0.0, 0.5, 0.5},
-	{0.5, 0.0, 1.0, 0.5},
-	{0.0, 0.5, 0.5, 1.0},
-	{0.5, 0.5, 1.0, 1.0},
-}
-
-func (ot *Tree) Tin() *Tin {
-	return &Tin{
-		EPSG:      ot.EPSG,
-		Vertices:  ot.vertices,
-		Triangles: ot.triangles,
-		Min:       ot.Min,
-		Max:       ot.Max,
-	}
-}
-
-func (ot *Tree) FindUnused() map[int32]struct{} {
-
-	used := make(map[int32]struct{}, len(ot.triangles))
-	for i := int32(0); i < int32(len(ot.triangles)); i++ {
-		used[i] = struct{}{}
-	}
-
-	stack := []int32{1}
-
-	for len(stack) > 0 {
-		top := stack[len(stack)-1]
-		stack = stack[:len(stack)-1]
-
-		if top > 0 { // node
-			if index := ot.index[top:]; len(index) > 7 {
-				for _, idx := range index[:8] {
-					if idx != 0 {
-						stack = append(stack, idx)
-					}
-				}
-			}
-		} else { // leaf
-			pos := -top - 1
-			n := ot.index[pos]
-			indices := ot.index[pos+1 : pos+1+n]
-			for _, idx := range indices {
-				delete(used, idx)
-			}
-		}
-	}
-
-	return used
-}
-
-func (ot *Tree) Value(x, y float64) (float64, bool) {
-
-	// out of bounding box
-	if x < ot.Min.X || ot.Max.X < x ||
-		y < ot.Min.Y || ot.Max.Y < y {
-		return 0, false
-	}
-
-	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
-
-	stack := []boxFrame{{1, all}}
-
-	for len(stack) > 0 {
-		top := stack[len(stack)-1]
-		stack = stack[:len(stack)-1]
-
-		if top.pos > 0 { // node
-			if index := ot.index[top.pos:]; len(index) > 7 {
-				for i := 0; i < 4; i++ {
-					a := index[i]
-					b := index[i+4]
-					if a == 0 && b == 0 {
-						continue
-					}
-					dx := top.X2 - top.X1
-					dy := top.Y2 - top.Y1
-					nbox := Box2D{
-						dx*scale[i][0] + top.X1,
-						dy*scale[i][1] + top.Y1,
-						dx*scale[i][2] + top.X1,
-						dy*scale[i][3] + top.Y1,
-					}
-					if nbox.Contains(x, y) {
-						if a != 0 {
-							stack = append(stack, boxFrame{a, nbox})
-						}
-						if b != 0 {
-							stack = append(stack, boxFrame{b, nbox})
-						}
-						break
-					}
-				}
-			}
-		} else { // leaf
-			pos := -top.pos - 1
-			n := ot.index[pos]
-			indices := ot.index[pos+1 : pos+1+n]
-
-			for _, idx := range indices {
-				tri := ot.triangles[idx]
-				t := Triangle{
-					ot.vertices[tri[0]],
-					ot.vertices[tri[1]],
-					ot.vertices[tri[2]],
-				}
-				if t.Contains(x, y) {
-					return t.Plane3D().Z(x, y), true
-				}
-			}
-		}
-	}
-
-	return 0, false
-}
-
-// Vertical does a vertical cross cut from (x1, y1) to (x2, y2).
-func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {
-
-	box := Box2D{
-		X1: math.Min(x1, x2),
-		Y1: math.Min(y1, y2),
-		X2: math.Max(x1, x2),
-		Y2: math.Max(y1, y2),
-	}
-
-	// out of bounding box
-	if box.X2 < ot.Min.X || ot.Max.X < box.X1 ||
-		box.Y2 < ot.Min.Y || ot.Max.Y < box.Y1 {
-		return
-	}
-
-	line := NewPlane2D(x1, y1, x2, y2)
-
-	dupes := map[int32]struct{}{}
-
-	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
-	//log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1))
-	//log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1))
-
-	stack := []boxFrame{{1, all}}
-
-	for len(stack) > 0 {
-		top := stack[len(stack)-1]
-		stack = stack[:len(stack)-1]
-
-		if top.pos > 0 { // node
-			if index := ot.index[top.pos:]; len(index) > 7 {
-				for i := 0; i < 4; i++ {
-					a := index[i]
-					b := index[i+4]
-					if a == 0 && b == 0 {
-						continue
-					}
-					dx := top.X2 - top.X1
-					dy := top.Y2 - top.Y1
-					nbox := Box2D{
-						dx*scale[i][0] + top.X1,
-						dy*scale[i][1] + top.Y1,
-						dx*scale[i][2] + top.X1,
-						dy*scale[i][3] + top.Y1,
-					}
-					if nbox.Intersects(box) && nbox.IntersectsPlane(line) {
-						if a != 0 {
-							stack = append(stack, boxFrame{a, nbox})
-						}
-						if b != 0 {
-							stack = append(stack, boxFrame{b, nbox})
-						}
-					}
-				}
-			}
-		} else { // leaf
-			pos := -top.pos - 1
-			n := ot.index[pos]
-			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]],
-				}
-
-				v0 := line.Eval(t[0].X, t[0].Y)
-				v1 := line.Eval(t[1].X, t[1].Y)
-				v2 := line.Eval(t[2].X, t[2].Y)
-
-				if onPlane(v0) || onPlane(v1) || onPlane(v2) ||
-					sides(sides(sides(0, v0), v1), v2) == 3 {
-					fn(&t)
-				}
-				dupes[idx] = struct{}{}
-			}
-		}
-	}
-}
-
-// Horizontal does a horizontal cross cut.
-func (ot *Tree) Horizontal(h float64, fn func(*Triangle)) {
-
-	if h < ot.Min.Z || ot.Max.Z < h {
-		return
-	}
-
-	type frame struct {
-		pos int32
-		min float64
-		max float64
-	}
-
-	dupes := map[int32]struct{}{}
-
-	stack := []frame{{1, ot.Min.Z, ot.Max.Z}}
-
-	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
-			}
-			if pos < int32(len(ot.index)) {
-				if index := ot.index[pos:]; len(index) > 3 {
-					stack = append(stack,
-						frame{index[0], min, max},
-						frame{index[1], min, max},
-						frame{index[2], min, max},
-						frame{index[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)
-				}
-			}
-		}
-	}
-}