--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/mesh/tin.go	Tue Nov 05 14:30:22 2019 +0100
@@ -0,0 +1,271 @@
+// 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 mesh
+
+import (
+	"bytes"
+	"encoding/binary"
+	"errors"
+	"fmt"
+	"io"
+	"log"
+	"math"
+
+	"gemma.intevation.de/gemma/pkg/models"
+	"gemma.intevation.de/gemma/pkg/wkb"
+)
+
+var (
+	errNoByteSlice   = errors.New("Not a byte slice")
+	errTooLessPoints = errors.New("Too less points")
+)
+
+// Tin stores a mesh of triangles with common vertices.
+type Tin struct {
+	// EPSG holds the projection.
+	EPSG uint32
+	// Vertices are the shared vertices.
+	Vertices []Vertex
+	// Triangles are the triangles.
+	Triangles [][]int32
+
+	// Min is the lower left corner of the bbox.
+	Min Vertex
+	// Max is the upper right corner of the bbox.
+	Max Vertex
+}
+
+func (t *Tin) Clip(polygon *Polygon) map[int32]struct{} {
+	var tree STRTree
+	tree.Build(t)
+	return tree.Clip(polygon)
+}
+
+// FromWKB constructs the TIN from a WKB representation.
+// Shared vertices are identified and referenced by the
+// same index.
+func (t *Tin) FromWKB(data []byte) error {
+	log.Printf("info: 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 == wkb.NDR:
+		order = binary.LittleEndian
+	case endian == wkb.XDR:
+		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 != wkb.TinZ:
+		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 == wkb.NDR:
+			order = binary.LittleEndian
+		case endian == wkb.XDR:
+			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 != wkb.TriangleZ:
+			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("info: bbox: [[%f, %f], [%f, %f]]\n",
+		min.X, min.Y, max.X, max.Y)
+
+	*t = Tin{
+		EPSG:      models.WGS84,
+		Vertices:  vertices,
+		Triangles: triangles,
+		Min:       min,
+		Max:       max,
+	}
+
+	return nil
+}
+
+// Scan implements the sql.Scanner interface.
+func (t *Tin) Scan(raw interface{}) error {
+	if raw == nil {
+		return nil
+	}
+	data, ok := raw.([]byte)
+	if !ok {
+		return errNoByteSlice
+	}
+	return t.FromWKB(data)
+}
+
+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
+	}
+
+	var err error
+	vwrite := func(v float64) {
+		if err == nil {
+			err = binary.Write(w, binary.LittleEndian, math.Float64bits(v))
+		}
+	}
+
+	for _, v := range t.Vertices {
+		vwrite(v.X)
+		vwrite(v.Y)
+		vwrite(v.Z)
+	}
+
+	if err != nil {
+		return err
+	}
+	log.Printf("info: 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("info: compressed tin indices in bytes: %d (%d)\n",
+		written, 3*4*len(t.Triangles))
+
+	return nil
+}