# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1571651512 -7200
# Node ID 26ce6fbfcd7c65da4c428ffa38976c038631e11b
# Parent  3c68d1572cab9bf012ed4bf1d16ebef24dfbdbeb# Parent  e16babbed40d31e5284c2b587722a849f3b00436
Merged direct-diff branch back into default.

diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/controllers/diff.go
--- a/pkg/controllers/diff.go	Mon Oct 21 11:31:41 2019 +0200
+++ b/pkg/controllers/diff.go	Mon Oct 21 11:51:52 2019 +0200
@@ -21,7 +21,6 @@
 	"log"
 	"net/http"
 	"runtime"
-	"sync"
 	"time"
 
 	"gemma.intevation.de/gemma/pkg/auth"
@@ -92,6 +91,22 @@
     4326
   )
 `
+	commonDiffBBoxSQL = `
+WITH joined AS (
+  SELECT
+    sr.area      AS area,
+    sr.date_info AS date_info
+  FROM waterway.sounding_results sr
+  WHERE sr.bottleneck_id = $1
+),
+bbox AS (
+  SELECT ST_Extent(ST_intersection(
+    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
+    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
+  )) AS area
+)
+SELECT ST_XMin(area), ST_YMin(area), ST_XMax(area), ST_YMax(area) FROM bbox
+`
 )
 
 type (
@@ -162,6 +177,38 @@
 			dci.Minuend.Format(common.DateFormat))
 	}
 
+	epsg := minuendTree.EPSG()
+
+	var box octree.Box2D
+
+	switch err := conn.QueryRowContext(
+		ctx,
+		commonDiffBBoxSQL,
+		dci.Bottleneck,
+		epsg,
+		dci.Minuend.Time,
+		dci.Subtrahend.Time,
+	).Scan(&box.X1, &box.Y1, &box.X2, &box.Y2); {
+	case err == sql.ErrNoRows:
+		return 0, errors.New("No such intersection")
+	case err != nil:
+		return 0, err
+	}
+
+	if box.Empty() {
+		return 0, errors.New("Intersection is empty")
+	}
+
+	log.Printf("info: bbox of intersection: (%.2f, %.2f) - (%.2f, %.2f)\n",
+		box.X1, box.Y1, box.X2, box.Y2)
+
+	start = time.Now()
+	raster := octree.NewRaster(box, isoCellSize)
+	raster.Rasterize(minuendTree.Value)
+	log.Printf("info: rasterizing minuend took %v\n", time.Since(start))
+
+	minuendTree = nil
+
 	start = time.Now()
 
 	subtrahendTree, err := octree.FromCache(
@@ -180,85 +227,15 @@
 	}
 
 	// We need a slow path implementation for this.
-	epsg := minuendTree.EPSG()
 	if epsg != subtrahendTree.EPSG() {
 		return 0, errors.New("Calculating differences between two different " +
 			"EPSG code meshes are not supported, yet.")
 	}
 
 	start = time.Now()
-	points := minuendTree.Diff(subtrahendTree)
+	raster.Diff(subtrahendTree.Value)
 	log.Printf("info: A - B took %v\n", time.Since(start))
-
-	minuendTree, subtrahendTree = nil, nil
-
-	// The Triangulation and the loading of the clipping
-	// polygon can be done concurrently.
-
-	jobs := make(chan func())
-
-	wg := new(sync.WaitGroup)
-	for i := 0; i < 2; i++ {
-		wg.Add(1)
-		go func() {
-			defer wg.Done()
-			for job := range jobs {
-				job()
-			}
-		}()
-	}
-
-	var (
-		tri     *octree.Triangulation
-		triErr  error
-		clip    *octree.Polygon
-		clipErr error
-	)
-
-	jobs <- func() {
-		start := time.Now()
-		tri, triErr = points.Triangulate()
-		log.Printf("info: triangulation took %v\n", time.Since(start))
-	}
-
-	jobs <- func() {
-		start := time.Now()
-		clip, clipErr = octree.LoadClippingPolygon(
-			ctx, conn,
-			epsg,
-			dci.Bottleneck,
-			dci.Minuend.Time,
-			dci.Subtrahend.Time)
-		log.Printf("info: loading clipping polygon took %v\n", time.Since(start))
-	}
-	close(jobs)
-	wg.Wait()
-
-	switch {
-	case triErr != nil && clipErr != nil:
-		return 0, fmt.Errorf("%v %v", triErr, clipErr)
-	case triErr != nil:
-		return 0, triErr
-	case clipErr != nil:
-		return 0, clipErr
-	}
-
-	start = time.Now()
-	tin := tri.Tin()
-	removed := tin.Clip(clip)
-	clip = nil
-	log.Printf("info: clipping TIN took %v\n", time.Since(start))
-
-	log.Printf("info: Number of triangles to clip: %d\n", len(removed))
-
-	start = time.Now()
-	var tree octree.STRTree
-
-	tree.BuildWithout(tin, removed)
-
-	log.Printf("info: Building final mesh took: %v\n", time.Since(start))
-
-	start = time.Now()
+	subtrahendTree = nil
 
 	// XXX: Maybe we should start this transaction earlier!?
 	var tx *sql.Tx
@@ -267,6 +244,13 @@
 	}
 	defer tx.Rollback()
 
+	zMin, zMax, ok := raster.ZExtent()
+	if !ok {
+		return 0, errors.New("Scans do not have common points")
+	}
+
+	log.Printf("info: z range: %.3f - %.3f\n", zMin, zMax)
+
 	var heights []float64
 
 	heights, err = octree.LoadClassBreaks(
@@ -275,12 +259,12 @@
 	if err != nil {
 		log.Printf("warn: Loading class breaks failed: %v\n", err)
 		err = nil
-		heights = octree.SampleDiffHeights(tin.Min.Z, tin.Max.Z, contourStep)
+		heights = octree.SampleDiffHeights(zMin, zMax, contourStep)
 	} else {
-		heights = octree.ExtrapolateClassBreaks(heights, tin.Min.Z, tin.Max.Z)
+		heights = octree.ExtrapolateClassBreaks(heights, zMin, zMax)
 	}
 
-	log.Printf("info: z range: %.3f - %.3f\n", tin.Min.Z, tin.Max.Z)
+	heights = common.DedupFloat64s(heights)
 
 	log.Printf("info: num heights: %d\n", len(heights))
 
@@ -302,9 +286,9 @@
 		return 0, err
 	}
 
-	heights = common.DedupFloat64s(heights)
+	areas := raster.Trace(heights)
 
-	areas := octree.TraceAreas(heights, isoCellSize, tin.Min, tin.Max, tree.Value)
+	raster = nil
 
 	var size int
 
@@ -327,7 +311,7 @@
 	log.Printf("info: Transferred WKB size: %.2fMB.\n",
 		float64(size)/(1024*1024))
 
-	log.Printf("info: calculating and storing iso lines took %v\n",
+	log.Printf("info: calculating and storing iso areas took %v\n",
 		time.Since(start))
 
 	if err = tx.Commit(); err != nil {
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/imports/isr.go
--- a/pkg/imports/isr.go	Mon Oct 21 11:31:41 2019 +0200
+++ b/pkg/imports/isr.go	Mon Oct 21 11:51:52 2019 +0200
@@ -184,7 +184,17 @@
 		}
 
 		// Calculate and store the iso areas.
-		areas := octree.TraceAreas(hs, isoCellSize, tree.Min(), tree.Max(), tree.Value)
+		box := octree.Box2D{
+			X1: tree.Min().X,
+			Y1: tree.Min().Y,
+			X2: tree.Max().X,
+			Y2: tree.Max().Y,
+		}
+
+		raster := octree.NewRaster(box, isoCellSize)
+		raster.Rasterize(tree.Value)
+		areas := raster.Trace(hs)
+
 		for i, a := range areas {
 			if len(a) == 0 {
 				continue
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/imports/sr.go
--- a/pkg/imports/sr.go	Mon Oct 21 11:31:41 2019 +0200
+++ b/pkg/imports/sr.go	Mon Oct 21 11:51:52 2019 +0200
@@ -893,7 +893,16 @@
 			time.Since(total))
 	}()
 
-	areas := octree.TraceAreas(heights, isoCellSize, tree.Min(), tree.Max(), tree.Value)
+	box := octree.Box2D{
+		X1: tree.Min().X,
+		Y1: tree.Min().Y,
+		X2: tree.Max().X,
+		Y2: tree.Max().Y,
+	}
+
+	raster := octree.NewRaster(box, isoCellSize)
+	raster.Rasterize(tree.Value)
+	areas := raster.Trace(heights)
 
 	return storeAreas(
 		ctx, tx, feedback,
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/octree/areas.go
--- a/pkg/octree/areas.go	Mon Oct 21 11:31:41 2019 +0200
+++ b/pkg/octree/areas.go	Mon Oct 21 11:51:52 2019 +0200
@@ -14,16 +14,10 @@
 package octree
 
 import (
-	"log"
-	"math"
 	"runtime"
 	"sync"
-	"time"
-
-	"github.com/fogleman/contourmap"
 
 	"gemma.intevation.de/gemma/pkg/common"
-	"gemma.intevation.de/gemma/pkg/wkb"
 )
 
 func GenerateRandomVertices(
@@ -92,300 +86,3 @@
 	close(out)
 	<-done
 }
-
-func TraceAreas(
-	heights []float64,
-	cellSize float64,
-	min, max Vertex,
-	eval func(float64, float64) (float64, bool),
-) []wkb.MultiPolygonGeom {
-
-	width := max.X - min.X
-	height := max.Y - min.Y
-
-	log.Printf("info: Width/Height: %.2f / %.2f\n", width, height)
-
-	xcells := int(math.Ceil(width / cellSize))
-	ycells := int(math.Ceil(height / cellSize))
-
-	log.Printf("info: Raster size: (%d, %d)\n", xcells, ycells)
-
-	start := time.Now()
-
-	// Add border for closing
-	raster := make([]float64, (xcells+2)*(ycells+2))
-
-	// prefill for no data
-	const nodata = -math.MaxFloat64
-	for i := range raster {
-		raster[i] = nodata
-	}
-
-	// rasterize the height model
-
-	var wg sync.WaitGroup
-
-	rows := make(chan int)
-
-	rasterRow := func() {
-		defer wg.Done()
-		quat := 0.25 * cellSize
-		for i := range rows {
-			pos := (i+1)*(xcells+2) + 1
-			row := raster[pos : pos+xcells]
-			py := min.Y + float64(i)*cellSize + cellSize/2
-			px := min.X + 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 += cellSize
-			}
-		}
-	}
-
-	for n := runtime.NumCPU(); n >= 1; n-- {
-		wg.Add(1)
-		go rasterRow()
-	}
-
-	for i := 0; i < ycells; i++ {
-		rows <- i
-	}
-	close(rows)
-
-	wg.Wait()
-	log.Printf("info: Rastering took %v\n", time.Since(start))
-
-	start = time.Now()
-
-	tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
-
-	areas := make([]wkb.MultiPolygonGeom, len(heights))
-
-	// TODO: Check if this correct!
-	reprojX := common.Linear(0.5, min.X, 1.5, min.X+cellSize)
-	reprojY := common.Linear(0.5, min.Y, 1.5, min.Y+cellSize)
-
-	cnts := make(chan int)
-
-	doContours := func() {
-		defer wg.Done()
-		for hIdx := range cnts {
-			c := tracer.Contours(heights[hIdx])
-
-			if len(c) == 0 {
-				continue
-			}
-
-			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
-}
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/octree/polygon.go
--- a/pkg/octree/polygon.go	Mon Oct 21 11:31:41 2019 +0200
+++ b/pkg/octree/polygon.go	Mon Oct 21 11:51:52 2019 +0200
@@ -15,13 +15,10 @@
 
 import (
 	"bytes"
-	"context"
-	"database/sql"
 	"encoding/binary"
 	"fmt"
 	"log"
 	"math"
-	"time"
 
 	"github.com/tidwall/rtree"
 
@@ -42,56 +39,11 @@
 )
 
 const (
-	clippingPolygonSQL = `
-WITH joined AS (
-  SELECT
-    sr.area      AS area,
-    sr.date_info AS date_info
-  FROM waterway.sounding_results sr
-  WHERE sr.bottleneck_id = $1
-)
-SELECT ST_AsBinary(
-  ST_Buffer(ST_intersection(
-    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
-    (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
-  ), 0.1)
-  ) AS area
-`
-)
-
-const (
 	IntersectionInside IntersectionType = iota
 	IntersectionOutSide
 	IntersectionOverlaps
 )
 
-func LoadClippingPolygon(
-	ctx context.Context,
-	conn *sql.Conn,
-	epsg uint32,
-	bottleneck string,
-	first, second time.Time,
-) (*Polygon, error) {
-
-	var clip []byte
-
-	if err := conn.QueryRowContext(
-		ctx, clippingPolygonSQL,
-		bottleneck,
-		epsg,
-		first, second,
-	).Scan(&clip); err != nil {
-		return nil, err
-	}
-
-	var polygon Polygon
-	if err := polygon.FromWKB(clip); err != nil {
-		return nil, err
-	}
-	polygon.Indexify()
-	return &polygon, nil
-}
-
 func (ls lineSegment) Rect(interface{}) ([]float64, []float64) {
 
 	var min, max [2]float64
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/octree/raster.go
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/raster.go	Mon Oct 21 11:51:52 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
+}
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/octree/slice.go
--- a/pkg/octree/slice.go	Mon Oct 21 11:31:41 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,105 +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, 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 (
-	"runtime"
-	"sync"
-)
-
-type PointMap map[Point]float64
-
-func (pm PointMap) Triangulate() (*Triangulation, error) {
-	vertices := make([]Vertex, len(pm))
-	var i int
-	for p, z := range pm {
-		vertices[i] = Vertex{X: p.X, Y: p.Y, Z: z}
-		i++
-	}
-	return Triangulate(vertices)
-}
-
-func sliceWork(
-	vs []Vertex,
-	dst PointMap,
-	fn func([]Vertex, func([]Vertex) []Vertex),
-) {
-	n := runtime.NumCPU()
-
-	wg := new(sync.WaitGroup)
-
-	slices := make(chan []Vertex)
-	out := make(chan []Vertex)
-
-	pool := make(chan []Vertex, n)
-
-	const pageSize = 2048
-
-	turn := func(p []Vertex) []Vertex {
-		if p != nil {
-			out <- p
-		}
-		select {
-		case p = <-pool:
-		default:
-			p = make([]Vertex, 0, pageSize)
-		}
-		return p
-	}
-
-	for i := 0; i < n; i++ {
-		wg.Add(1)
-		go func() {
-			defer wg.Done()
-			for slice := range slices {
-				fn(slice, turn)
-			}
-		}()
-	}
-	done := make(chan struct{})
-	go func() {
-		defer close(done)
-		for s := range out {
-			for i := range s {
-				v := &s[i]
-				key := Point{X: v.X, Y: v.Y}
-				if z, found := dst[key]; found {
-					dst[key] = (z + v.Z) * 0.5
-				} else {
-					dst[key] = v.Z
-				}
-			}
-			select {
-			case pool <- s[:0:pageSize]:
-			default:
-			}
-		}
-	}()
-
-	size := len(vs)/n + 1
-	for len(vs) > 0 {
-		var l int
-		if len(vs) < size {
-			l = len(vs)
-		} else {
-			l = size
-		}
-		slices <- vs[:l]
-		vs = vs[l:]
-	}
-	close(slices)
-	wg.Wait()
-	close(out)
-	<-done
-}
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/octree/strtree.go
--- a/pkg/octree/strtree.go	Mon Oct 21 11:31:41 2019 +0200
+++ b/pkg/octree/strtree.go	Mon Oct 21 11:51:52 2019 +0200
@@ -484,51 +484,3 @@
 	}
 	return int32(-(pos + 1))
 }
-
-func (s *STRTree) Diff(other *STRTree) PointMap {
-
-	firstVs, secondVs := s.tin.Vertices, other.tin.Vertices
-
-	result := make(PointMap, len(firstVs)+len(secondVs))
-
-	sliceWork(
-		firstVs,
-		result,
-		func(slice []Vertex, turn func([]Vertex) []Vertex) {
-			p := turn(nil)
-			for i := range slice {
-				v := &slice[i]
-				if z, found := other.Value(v.X, v.Y); found {
-					p = append(p, Vertex{v.X, v.Y, v.Z - z})
-					if len(p) == cap(p) {
-						p = turn(p)
-					}
-				}
-			}
-			if len(p) > 0 {
-				turn(p)
-			}
-		})
-
-	sliceWork(
-		secondVs,
-		result,
-		func(
-			slice []Vertex, turn func([]Vertex) []Vertex) {
-			p := turn(nil)
-			for i := range slice {
-				v := &slice[i]
-				if z, found := s.Value(v.X, v.Y); found {
-					p = append(p, Vertex{v.X, v.Y, z - v.Z})
-					if len(p) == cap(p) {
-						p = turn(p)
-					}
-				}
-			}
-			if len(p) > 0 {
-				turn(p)
-			}
-		})
-
-	return result
-}
diff -r 3c68d1572cab -r 26ce6fbfcd7c pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Mon Oct 21 11:31:41 2019 +0200
+++ b/pkg/octree/vertex.go	Mon Oct 21 11:51:52 2019 +0200
@@ -113,6 +113,16 @@
 		a.Y1 >= b.Y1 && a.Y2 <= b.Y2
 }
 
+func (a Box2D) Size() (float64, float64) {
+	return a.X2 - a.X1, a.Y2 - a.Y1
+}
+
+func (a Box2D) Empty() bool {
+	const eps = 0.0000001
+	return math.Abs(a.X2-a.X1) < eps &&
+		math.Abs(a.Y2-a.Y1) < eps
+}
+
 func (p Plane3D) Z(x, y float64) float64 {
 	// p.A*x + p.B*y + p.C*z + p.D = 0
 	return -(p.A*x + p.B*y + p.D) / p.C