--- a/pkg/octree/areas.go	Sat Oct 19 23:07:59 2019 +0200
+++ b/pkg/octree/areas.go	Mon Oct 21 02:01:56 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
-}