--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/mesh/polygon.go	Tue Nov 05 14:30:22 2019 +0100
@@ -0,0 +1,505 @@
+// 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"
+	"fmt"
+	"log"
+	"math"
+
+	"github.com/tidwall/rtree"
+
+	"gemma.intevation.de/gemma/pkg/wkb"
+)
+
+type (
+	ring []float64
+
+	Polygon struct {
+		rings   []ring
+		indices []*rtree.RTree
+	}
+
+	IntersectionType byte
+
+	lineSegment []float64
+)
+
+const (
+	IntersectionInside IntersectionType = iota
+	IntersectionOutSide
+	IntersectionOverlaps
+)
+
+func (ls lineSegment) Rect(interface{}) ([]float64, []float64) {
+
+	var min, max [2]float64
+
+	if ls[0] < ls[2] {
+		min[0] = ls[0]
+		max[0] = ls[2]
+	} else {
+		min[0] = ls[2]
+		max[0] = ls[0]
+	}
+
+	if ls[1] < ls[3] {
+		min[1] = ls[1]
+		max[1] = ls[3]
+	} else {
+		min[1] = ls[3]
+		max[1] = ls[1]
+	}
+
+	return min[:], max[:]
+}
+
+func (p *Polygon) Indexify() {
+	indices := make([]*rtree.RTree, len(p.rings))
+
+	for i := range indices {
+		index := rtree.New(nil)
+		indices[i] = index
+
+		rng := p.rings[i]
+
+		for i := 0; i < len(rng); i += 2 {
+			var ls lineSegment
+			if i+4 <= len(rng) {
+				ls = lineSegment(rng[i : i+4])
+			} else {
+				ls = []float64{rng[i], rng[i+1], rng[0], rng[1]}
+			}
+			index.Insert(ls)
+		}
+	}
+
+	p.indices = indices
+}
+
+func (ls lineSegment) intersects(a Box2D) bool {
+	p1x := ls[0]
+	p1y := ls[1]
+	p2x := ls[2]
+	p2y := ls[3]
+
+	left := a.X1
+	right := a.X2
+	top := a.Y1
+	bottom := a.Y2
+
+	// The direction of the ray
+	dx := p2x - p1x
+	dy := p2y - p1y
+
+	min, max := 0.0, 1.0
+
+	var t0, t1 float64
+
+	// Left and right sides.
+	// - If the line is parallel to the y axis.
+	if dx == 0 {
+		if p1x < left || p1x > right {
+			return false
+		}
+	} else {
+		// - Make sure t0 holds the smaller value by checking the direction of the line.
+		if dx > 0 {
+			t0 = (left - p1x) / dx
+			t1 = (right - p1x) / dx
+		} else {
+			t1 = (left - p1x) / dx
+			t0 = (right - p1x) / dx
+		}
+
+		if t0 > min {
+			min = t0
+		}
+		if t1 < max {
+			max = t1
+		}
+		if min > max || max < 0 {
+			return false
+		}
+	}
+
+	// The top and bottom side.
+	// - If the line is parallel to the x axis.
+	if dy == 0 {
+		if p1y < top || p1y > bottom {
+			return false
+		}
+	} else {
+		// - Make sure t0 holds the smaller value by checking the direction of the line.
+		if dy > 0 {
+			t0 = (top - p1y) / dy
+			t1 = (bottom - p1y) / dy
+		} else {
+			t1 = (top - p1y) / dy
+			t0 = (bottom - p1y) / dy
+		}
+
+		if t0 > min {
+			min = t0
+		}
+		if t1 < max {
+			max = t1
+		}
+		if min > max || max < 0 {
+			return false
+		}
+	}
+
+	// The point of intersection
+	// ix = p1x + dx*min
+	// iy = p1y + dy*min
+	return true
+}
+
+func (ls lineSegment) intersectsLineSegment(o lineSegment) bool {
+
+	p0 := ls[:2]
+	p1 := ls[2:4]
+	p2 := o[:2]
+	p3 := o[2:4]
+
+	s10x := p1[0] - p0[0]
+	s10y := p1[1] - p0[1]
+	s32x := p3[0] - p2[0]
+	s32y := p3[1] - p2[1]
+
+	den := s10x*s32y - s32x*s10y
+
+	if den == 0 {
+		return false
+	}
+
+	denPos := den > 0
+
+	s02x := p0[0] - p2[0]
+	s02y := p0[1] - p2[1]
+
+	sNum := s10x*s02y - s10y*s02x
+	if sNum < 0 == denPos {
+		return false
+	}
+
+	tNum := s32x*s02y - s32y*s02x
+	if tNum < 0 == denPos {
+		return false
+	}
+
+	if sNum > den == denPos || tNum > den == denPos {
+		return false
+	}
+
+	// t := tNum / den
+	// intersection at( p0[0] + (t * s10x), p0[1] + (t * s10y) )
+	return true
+}
+
+func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType {
+
+	if len(p.rings) == 0 {
+		return IntersectionOutSide
+	}
+
+	for _, index := range p.indices {
+		var intersects bool
+		index.Search(box, func(item rtree.Item) bool {
+			if item.(lineSegment).intersects(box) {
+				intersects = true
+				return false
+			}
+			return true
+		})
+		if intersects {
+			return IntersectionOverlaps
+		}
+	}
+
+	// No intersection -> check inside or outside
+	// if an abritrary point  is inside or not.
+
+	// Check holes first: inside a hole means outside.
+	if len(p.rings) > 1 {
+		for _, hole := range p.rings[1:] {
+			if contains(hole, box.X1, box.Y1) {
+				return IntersectionOutSide
+			}
+		}
+	}
+
+	// Check shell
+	if contains(p.rings[0], box.X1, box.Y1) {
+		return IntersectionInside
+	}
+	return IntersectionOutSide
+}
+
+func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType {
+	box := t.BBox()
+	for _, index := range p.indices {
+		var intersects bool
+		index.Search(box, func(item rtree.Item) bool {
+			ls := item.(lineSegment)
+			other := make(lineSegment, 4)
+			for i := range t {
+				other[0] = t[i].X
+				other[1] = t[i].Y
+				other[2] = t[(i+1)%len(t)].X
+				other[3] = t[(i+1)%len(t)].Y
+				if ls.intersectsLineSegment(other) {
+					intersects = true
+					return false
+				}
+			}
+			return true
+		})
+		if intersects {
+			return IntersectionOverlaps
+		}
+	}
+	// No intersection -> check inside or outside
+	// if an abritrary point  is inside or not.
+	pX, pY := t[0].X, t[0].Y
+
+	// Check holes first: inside a hole means outside.
+	if len(p.rings) > 1 {
+		for _, hole := range p.rings[1:] {
+			if contains(hole, pX, pY) {
+				return IntersectionOutSide
+			}
+		}
+	}
+
+	// Check shell
+	if contains(p.rings[0], pX, pY) {
+		return IntersectionInside
+	}
+	return IntersectionOutSide
+}
+
+func (rng ring) length() int {
+	return len(rng) / 2
+}
+
+func (rng ring) point(i int) (float64, float64) {
+	i *= 2
+	return rng[i], rng[i+1]
+}
+
+type segments interface {
+	length() int
+	point(int) (float64, float64)
+}
+
+func contains(s segments, pX, pY float64) bool {
+
+	n := s.length()
+	if n < 3 {
+		return false
+	}
+
+	sX, sY := s.point(0)
+	eX, eY := s.point(n - 1)
+
+	const eps = 0.0000001
+
+	if math.Abs(sX-eX) > eps || math.Abs(sY-eY) > eps {
+		// It's not closed!
+		return false
+	}
+
+	var inside bool
+
+	for i := 1; i < n; i++ {
+		eX, eY := s.point(i)
+		if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) {
+			inside = !inside
+		}
+		sX, sY = eX, eY
+	}
+
+	return inside
+}
+
+// Using the raycast algorithm, this returns whether or not the passed in point
+// Intersects with the edge drawn by the passed in start and end points.
+// Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go
+func intersectsWithRaycast(pX, pY, sX, sY, eX, eY float64) bool {
+
+	// Always ensure that the the first point
+	// has a y coordinate that is less than the second point
+	if sY > eY {
+		// Switch the points if otherwise.
+		sX, sY, eX, eY = eX, eY, sX, sY
+	}
+
+	// Move the point's y coordinate
+	// outside of the bounds of the testing region
+	// so we can start drawing a ray
+	for pY == sY || pY == eY {
+		pY = math.Nextafter(pY, math.Inf(1))
+	}
+
+	// If we are outside of the polygon, indicate so.
+	if pY < sY || pY > eY {
+		return false
+	}
+
+	if sX > eX {
+		if pX > sX {
+			return false
+		}
+		if pX < eX {
+			return true
+		}
+	} else {
+		if pX > eX {
+			return false
+		}
+		if pX < sX {
+			return true
+		}
+	}
+
+	raySlope := (pY - sY) / (pX - sX)
+	diagSlope := (eY - sY) / (eX - sX)
+
+	return raySlope >= diagSlope
+}
+
+func (p *Polygon) NumVertices(ring int) int {
+	if ring < 0 || ring >= len(p.rings) {
+		return 0
+	}
+	return len(p.rings[ring]) / 2
+}
+
+func (p *Polygon) Vertices(ring int, fn func(float64, float64)) {
+	if ring < 0 || ring >= len(p.rings) {
+		return
+	}
+	rng := p.rings[ring]
+	for i := 0; i < len(rng); i += 2 {
+		fn(rng[i+0], rng[i+1])
+	}
+}
+
+func (p *Polygon) AsWKB() []byte {
+
+	size := 1 + 4 + 4
+	for _, r := range p.rings {
+		size += 4 + len(r)*8
+	}
+
+	buf := bytes.NewBuffer(make([]byte, 0, size))
+
+	binary.Write(buf, binary.LittleEndian, wkb.NDR)
+	binary.Write(buf, binary.LittleEndian, wkb.Polygon)
+	binary.Write(buf, binary.LittleEndian, uint32(len(p.rings)))
+
+	for _, r := range p.rings {
+		binary.Write(buf, binary.LittleEndian, uint32(len(r)/2))
+		for i := 0; i < len(r); i += 2 {
+			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0]))
+			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1]))
+		}
+	}
+
+	return buf.Bytes()
+}
+
+func (p *Polygon) FromWKB(data []byte) error {
+
+	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.Polygon:
+		return fmt.Errorf("unknown geometry type %x", geomType)
+	}
+
+	var numRings uint32
+	if err = binary.Read(r, order, &numRings); err != nil {
+		return err
+	}
+
+	rngs := make([]ring, numRings)
+
+	log.Printf("info: Number of rings: %d\n", len(rngs))
+
+	for rng := uint32(0); rng < numRings; rng++ {
+		var numVertices uint32
+		if err = binary.Read(r, order, &numVertices); err != nil {
+			return err
+		}
+
+		log.Printf("info: Number of vertices in ring %d: %d\n", rng, numVertices)
+
+		numVertices *= 2
+		vertices := make([]float64, numVertices)
+
+		for v := uint32(0); v < numVertices; v += 2 {
+			var lat, lon uint64
+			if err = binary.Read(r, order, &lat); err != nil {
+				return err
+			}
+			if err = binary.Read(r, order, &lon); err != nil {
+				return err
+			}
+			vertices[v] = math.Float64frombits(lat)
+			vertices[v+1] = math.Float64frombits(lon)
+		}
+
+		rngs[rng] = vertices
+	}
+
+	p.rings = rngs
+
+	return nil
+}
+
+func (p *Polygon) FromLineStringZ(ls LineStringZ) {
+	r := make([]float64, 2*len(ls))
+	var pos int
+	for i := range ls {
+		r[pos+0] = ls[i].X
+		r[pos+1] = ls[i].Y
+		pos += 2
+	}
+	p.rings = []ring{r}
+}