--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/raster.go	Mon Oct 21 02:01:56 2019 +0200
@@ -0,0 +1,404 @@
+// 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 octree
+
+import (
+	"log"
+	"math"
+	"runtime"
+	"sync"
+	"time"
+
+	"gemma.intevation.de/gemma/pkg/common"
+	"gemma.intevation.de/gemma/pkg/wkb"
+	"github.com/fogleman/contourmap"
+)
+
+type Raster struct {
+	BBox     Box2D
+	CellSize float64
+	XCells   int
+	YCells   int
+	Cells    []float64
+}
+
+const noData = -math.MaxFloat64
+
+func NewRaster(bbox Box2D, cellSize float64) *Raster {
+
+	width, height := bbox.Size()
+
+	log.Printf("info raster extent: %.2f / %.2f", width, height)
+
+	xCells := int(math.Ceil(width / cellSize))
+	yCells := int(math.Ceil(height / cellSize))
+
+	log.Printf("info raster size: %d / %d\n", xCells, yCells)
+
+	size := (xCells + 2) * (yCells + 2)
+	cells := make([]float64, size)
+	for i := range cells {
+		cells[i] = noData
+	}
+	return &Raster{
+		BBox:     bbox,
+		CellSize: cellSize,
+		XCells:   xCells,
+		YCells:   yCells,
+		Cells:    cells,
+	}
+}
+
+func (r *Raster) Rasterize(eval func(float64, float64) (float64, bool)) {
+	var wg sync.WaitGroup
+
+	rows := make(chan int)
+
+	rasterRow := func() {
+		defer wg.Done()
+		quat := 0.25 * r.CellSize
+		for i := range rows {
+			pos := (i+1)*(r.XCells+2) + 1
+			row := r.Cells[pos : pos+r.XCells]
+			py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2
+			px := r.BBox.X1 + r.CellSize/2
+			y1 := py - quat
+			y2 := py + quat
+			for j := range row {
+				var n int
+				var sum float64
+
+				if v, ok := eval(px-quat, y1); ok {
+					sum = v
+					n = 1
+				}
+				if v, ok := eval(px-quat, y2); ok {
+					sum += v
+					n++
+				}
+				if v, ok := eval(px+quat, y1); ok {
+					sum += v
+					n++
+				}
+				if v, ok := eval(px+quat, y2); ok {
+					sum += v
+					n++
+				}
+
+				if n > 0 {
+					row[j] = sum / float64(n)
+				}
+				px += r.CellSize
+			}
+		}
+	}
+
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go rasterRow()
+	}
+
+	for i := 0; i < r.YCells; i++ {
+		rows <- i
+	}
+	close(rows)
+}
+
+func (r *Raster) Diff(eval func(float64, float64) (float64, bool)) {
+	var wg sync.WaitGroup
+
+	rows := make(chan int)
+
+	rasterRow := func() {
+		defer wg.Done()
+		quat := 0.25 * r.CellSize
+		for i := range rows {
+			pos := (i+1)*(r.XCells+2) + 1
+			row := r.Cells[pos : pos+r.XCells]
+			py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2
+			px := r.BBox.X1 + r.CellSize/2
+			y1 := py - quat
+			y2 := py + quat
+			for j, old := range row {
+				// only diff where values
+				if old == noData {
+					px += r.CellSize
+					continue
+				}
+				var n int
+				var sum float64
+
+				if v, ok := eval(px-quat, y1); ok {
+					sum = v
+					n = 1
+				}
+				if v, ok := eval(px-quat, y2); ok {
+					sum += v
+					n++
+				}
+				if v, ok := eval(px+quat, y1); ok {
+					sum += v
+					n++
+				}
+				if v, ok := eval(px+quat, y2); ok {
+					sum += v
+					n++
+				}
+
+				if n > 0 {
+					row[j] -= sum / float64(n)
+				} else {
+					row[j] = noData
+				}
+
+				px += r.CellSize
+			}
+		}
+	}
+
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go rasterRow()
+	}
+
+	for i := 0; i < r.YCells; i++ {
+		rows <- i
+	}
+	close(rows)
+}
+
+func (r *Raster) ZExtent() (float64, float64, bool) {
+	min, max := math.MaxFloat64, -math.MaxFloat64
+	for _, v := range r.Cells {
+		if v == noData {
+			continue
+		}
+		if v < min {
+			min = v
+		}
+		if v > max {
+			max = v
+		}
+	}
+	return min, max, min != math.MaxFloat64
+}
+
+func (r *Raster) Trace(heights []float64) []wkb.MultiPolygonGeom {
+	start := time.Now()
+
+	tracer := contourmap.FromFloat64s(r.XCells+2, r.YCells+2, r.Cells)
+
+	areas := make([]wkb.MultiPolygonGeom, len(heights))
+
+	reprojX := common.Linear(0.5, r.BBox.X1, 1.5, r.BBox.X1+r.CellSize)
+	reprojY := common.Linear(0.5, r.BBox.Y1, 1.5, r.BBox.Y1+r.CellSize)
+
+	var wg sync.WaitGroup
+
+	cnts := make(chan int)
+
+	doContours := func() {
+		defer wg.Done()
+		for hIdx := range cnts {
+			if c := tracer.Contours(heights[hIdx]); len(c) > 0 {
+				areas[hIdx] = buildMultipolygon(c, reprojX, reprojY)
+			}
+		}
+	}
+
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go doContours()
+	}
+
+	for i := range heights {
+		cnts <- i
+	}
+	close(cnts)
+
+	wg.Wait()
+	log.Printf("info: Tracing areas took %v\n", time.Since(start))
+
+	return areas
+}
+
+type contour []contourmap.Point
+
+type bboxNode struct {
+	box      Box2D
+	cnt      contour
+	children []*bboxNode
+}
+
+func (cnt contour) contains(o contour) bool {
+	return contains(cnt, o[0].X, o[0].Y) ||
+		contains(cnt, o[len(o)/2].X, o[len(o)/2].Y)
+}
+
+func (cnt contour) length() int {
+	return len(cnt)
+}
+
+func (cnt contour) point(i int) (float64, float64) {
+	return cnt[i].X, cnt[i].Y
+}
+
+func (cnt contour) bbox() Box2D {
+	minX, minY := math.MaxFloat64, math.MaxFloat64
+	maxX, maxY := -math.MaxFloat64, -math.MaxFloat64
+
+	for _, p := range cnt {
+		if p.X < minX {
+			minX = p.X
+		}
+		if p.X > maxX {
+			maxX = p.X
+		}
+		if p.Y < minY {
+			minY = p.Y
+		}
+		if p.Y > maxY {
+			maxY = p.Y
+		}
+	}
+	return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY}
+}
+
+func (bn *bboxNode) insert(cnt contour, box Box2D) {
+	// check if children are inside new
+	var nr *bboxNode
+
+	for i, r := range bn.children {
+		if r.box.Inside(box) && cnt.contains(r.cnt) {
+			if nr == nil {
+				nr = &bboxNode{box: box, cnt: cnt}
+			}
+			nr.children = append(nr.children, r)
+			bn.children[i] = nil
+		}
+	}
+
+	// we have a new child
+	if nr != nil {
+		// compact the list
+		for i := len(bn.children) - 1; i >= 0; i-- {
+			if bn.children[i] == nil {
+				if i < len(bn.children)-1 {
+					copy(bn.children[i:], bn.children[i+1:])
+				}
+				bn.children[len(bn.children)-1] = nil
+				bn.children = bn.children[:len(bn.children)-1]
+			}
+		}
+		bn.children = append(bn.children, nr)
+		return
+	}
+
+	// check if new is inside an old
+	for _, r := range bn.children {
+		if box.Inside(r.box) && r.cnt.contains(cnt) {
+			r.insert(cnt, box)
+			return
+		}
+	}
+
+	// its a new child node.
+	nr = &bboxNode{box: box, cnt: cnt}
+	bn.children = append(bn.children, nr)
+}
+
+func (bn *bboxNode) insertRoot(cnt contour) {
+	bn.insert(cnt, cnt.bbox())
+}
+
+type bboxOutFunc func(contour, []contour)
+
+func (bn *bboxNode) generate(out bboxOutFunc) {
+
+	var grands []*bboxNode
+
+	holes := make([]contour, len(bn.children))
+
+	for i, ch := range bn.children {
+		holes[i] = ch.cnt
+		grands = append(grands, ch.children...)
+	}
+	out(bn.cnt, holes)
+
+	// the grand children are new polygons.
+	for _, grand := range grands {
+		grand.generate(out)
+	}
+}
+
+func (bn *bboxNode) generateRoot(out bboxOutFunc) {
+	for _, r := range bn.children {
+		r.generate(out)
+	}
+}
+
+func buildMultipolygon(
+	cnts []contourmap.Contour,
+	reprojX, reprojY func(float64) float64,
+) wkb.MultiPolygonGeom {
+
+	var forest bboxNode
+
+	for _, cnt := range cnts {
+		forest.insertRoot(contour(cnt))
+	}
+
+	//log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots))
+
+	var mp wkb.MultiPolygonGeom
+
+	out := func(sh contour, hls []contour) {
+
+		polygon := make(wkb.PolygonGeom, 1+len(hls))
+
+		// Handle shell
+		shell := make(wkb.LinearRingGeom, len(sh))
+		for j, pt := range sh {
+			shell[j] = wkb.PointGeom{
+				X: reprojX(pt.X),
+				Y: reprojY(pt.Y),
+			}
+		}
+		if shell.CCW() {
+			shell.Reverse()
+		}
+		polygon[0] = shell
+
+		// handle holes
+		for i, hl := range hls {
+			hole := make(wkb.LinearRingGeom, len(hl))
+			for j, pt := range hl {
+				hole[j] = wkb.PointGeom{
+					X: reprojX(pt.X),
+					Y: reprojY(pt.Y),
+				}
+			}
+			if !hole.CCW() {
+				hole.Reverse()
+			}
+			polygon[1+i] = hole
+		}
+
+		mp = append(mp, polygon)
+	}
+
+	forest.generateRoot(out)
+
+	return mp
+}