--- a/cmd/isoareas/main.go	Tue Oct 01 11:07:33 2019 +0200
+++ b/cmd/isoareas/main.go	Tue Oct 01 13:18:16 2019 +0200
@@ -15,17 +15,21 @@
 
 import (
 	"bufio"
+	"flag"
 	"fmt"
 	"io"
 	"log"
 	"math"
 	"math/rand"
 	"os"
+	"runtime"
 	"strconv"
 	"strings"
+	"sync"
 	"time"
 
 	svg "github.com/ajstarks/svgo"
+	"github.com/fogleman/contourmap"
 
 	"gemma.intevation.de/gemma/pkg/octree"
 )
@@ -66,21 +70,6 @@
 	return points, nil
 }
 
-func minMax(points []octree.Vertex) (octree.Vertex, octree.Vertex) {
-	if len(points) == 0 {
-		return octree.Vertex{}, octree.Vertex{}
-	}
-
-	min, max := points[0], points[0]
-
-	for _, v := range points {
-		min.Minimize(v)
-		max.Maximize(v)
-	}
-
-	return min, max
-}
-
 func check(err error) {
 	if err != nil {
 		log.Fatalf("error: %v\n", err)
@@ -108,6 +97,77 @@
 	}
 }
 
+func dumpContoursSVG(
+	w io.Writer,
+	contours [][]contourmap.Contour,
+) (err error) {
+
+	var (
+		minX = math.MaxFloat64
+		minY = math.MaxFloat64
+		maxX = -math.MaxFloat64
+		maxY = -math.MaxFloat64
+	)
+
+	for _, c := range contours {
+		for _, r := range c {
+			for _, p := range r {
+				if p.X < minX {
+					minX = p.X
+				}
+				if p.Y < minY {
+					minY = p.Y
+				}
+				if p.X > maxX {
+					maxX = p.X
+				}
+				if p.Y > maxY {
+					maxY = p.Y
+				}
+			}
+		}
+	}
+	ratio := (maxX - minX) / (maxY - minY)
+
+	log.Printf("ratio: %.2f\n", ratio)
+
+	const width = 50000
+	height := int(math.Ceil(width * ratio))
+
+	px := linear(minX, 0, maxX, width)
+	py := linear(minY, float64(height), maxY, 0)
+
+	out := bufio.NewWriter(w)
+	defer func() { err = out.Flush() }()
+
+	canvas := svg.New(out)
+	canvas.Start(width, height)
+
+	rnd := rand.New(rand.NewSource(42))
+
+	for _, c := range contours {
+
+		r := byte(rnd.Intn(256))
+		g := byte(rnd.Intn(256))
+		b := byte(rnd.Intn(256))
+		style := fmt.Sprintf("fill:#%02x%02x%02x", r, g, b)
+
+		for _, r := range c {
+			x := make([]int, len(r))
+			y := make([]int, len(r))
+			for i, p := range r {
+				x[i] = int(math.Round(px(p.X)))
+				y[i] = int(math.Round(py(p.Y)))
+			}
+
+			canvas.Polygon(x, y, style)
+		}
+	}
+
+	canvas.End()
+	return
+}
+
 func dumpPolygonsSVG(
 	w io.Writer,
 	min, max octree.Vertex,
@@ -218,6 +278,12 @@
 
 func main() {
 
+	cellSize := float64(0.5)
+
+	flag.Float64Var(&cellSize, "cellsize", 0.5, "width/height of raster cell [m]")
+
+	flag.Parse()
+
 	heights, err := octree.ParseClassBreaks(classBreaks)
 	check(err)
 	log.Printf("num classes: %d\n", len(heights))
@@ -229,22 +295,139 @@
 
 	log.Printf("num vertices: %d\n", len(xyz))
 
-	min, max := minMax(xyz)
+	start = time.Now()
+	tri, err := octree.Triangulate(xyz)
+	check(err)
+	log.Printf("triangulation took %v\n", time.Since(start))
+	tooLongEdge := tri.EstimateTooLong()
+	log.Printf("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
+
+	start = time.Now()
+	polygon, removed := tri.ConcaveHull(tooLongEdge)
+
+	polygonArea := polygon.Area()
+	if polygonArea < 0.0 { // counter clockwise
+		polygonArea = -polygonArea
+		polygon.Reverse()
+	}
+
+	var clippingPolygon octree.Polygon
+	clippingPolygon.FromLineStringZ(polygon)
+	clippingPolygon.Indexify()
+
+	tin := tri.Tin()
+	// tin.EPSG = epsg
+
+	var str octree.STRTree
+	str.Build(tin)
+
+	removed = str.Clip(&clippingPolygon)
+
+	builder := octree.NewBuilder(tin)
+	builder.Build(removed)
+
+	tree := builder.Tree()
+
+	min, max := tin.Min, tin.Max
 
 	log.Printf("(%.2f, %.2f, %.2f) - (%.2f, %.2f, %.2f)\n",
 		min.X, min.Y, min.Z,
 		max.X, max.Y, max.Z)
 
+	width := max.X - min.X
+	height := max.Y - min.Y
+
+	log.Printf("width/height: %.2f / %.2f\n", width, height)
+
+	xcells := int(math.Ceil(width / cellSize))
+	ycells := int(math.Ceil(height / cellSize))
+
+	log.Printf("raster size: (%d, %d)\n", xcells, ycells)
+
+	raster := make([]float64, xcells*ycells)
+
+	const closed = -math.MaxFloat64
+	for i := range raster {
+		raster[i] = closed
+	}
+
 	start = time.Now()
-	tri, err := octree.Triangulate(xyz)
-	check(err)
-	log.Printf("triangulation took %v\n", time.Since(start))
+
+	var wg sync.WaitGroup
+
+	rows := make(chan int)
 
-	log.Printf("number of triangles: %d\n", len(tri.Triangles)/3)
+	rasterRow := func() {
+		defer wg.Done()
+		for i := range rows {
+			pos := i * xcells
+			row := raster[pos : pos+xcells]
+			//log.Printf("len: %d\n", len(row))
+			py := min.Y + float64(i)*cellSize + cellSize/2
+			px := min.X + cellSize/2
+			for j := range row {
+				v, ok := tree.Value(px, py)
+				if ok {
+					row[j] = v
+				}
+				px += cellSize
+			}
+		}
+	}
 
 	start = time.Now()
-	polygons := intersectTriangles(tri, heights, min, max)
-	log.Printf("intersecting triangles took %v.\n", time.Since(start))
+
+	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("Rasterizing took %v.\n", time.Since(start))
+
+	start = time.Now()
+	cm := contourmap.FromFloat64s(xcells, ycells, raster).Closed()
+
+	start = time.Now()
+
+	contours := make([][]contourmap.Contour, len(heights))
+
+	cnts := make(chan int)
 
-	check(dumpPolygonsSVG(os.Stdout, min, max, polygons))
+	doContours := func() {
+		defer wg.Done()
+		for i := range cnts {
+			contours[i] = cm.Contours(heights[i])
+		}
+	}
+
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go doContours()
+	}
+
+	for i := range heights {
+		cnts <- i
+	}
+	close(cnts)
+
+	wg.Wait()
+
+	log.Printf("Calculating contours took %v.\n", time.Since(start))
+	check(dumpContoursSVG(os.Stdout, contours))
+
+	/*
+
+		start = time.Now()
+		polygons := intersectTriangles(tri, heights, min, max)
+		log.Printf("intersecting triangles took %v.\n", time.Since(start))
+
+		check(dumpPolygonsSVG(os.Stdout, min, max, polygons))
+	*/
 }

--- a/go.mod	Tue Oct 01 11:07:33 2019 +0200
+++ b/go.mod	Tue Oct 01 13:18:16 2019 +0200
@@ -6,6 +6,7 @@
 	github.com/ajstarks/svgo v0.0.0-20180226025133-644b8db467af
 	github.com/cockroachdb/apd v1.1.0 // indirect
 	github.com/etcd-io/bbolt v1.3.3
+	github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199
 	github.com/gofrs/uuid v3.2.0+incompatible // indirect
 	github.com/golang/snappy v0.0.1
 	github.com/gorilla/mux v1.7.3

--- a/go.sum	Tue Oct 01 11:07:33 2019 +0200
+++ b/go.sum	Tue Oct 01 13:18:16 2019 +0200
@@ -26,6 +26,8 @@
 github.com/dgryski/go-sip13 v0.0.0-20181026042036-e10d5fee7954/go.mod h1:vAd38F8PWV+bWy6jNmig1y/TA+kYO4g3RSRF0IAv0no=
 github.com/etcd-io/bbolt v1.3.3 h1:gSJmxrs37LgTqR/oyJBWok6k6SvXEUerFTbltIhXkBM=
 github.com/etcd-io/bbolt v1.3.3/go.mod h1:ZF2nL25h33cCyBtcyWeZ2/I3HQOfTP+0PIEvHjkjCrw=
+github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199 h1:kufr0u0RIG5ACpjFsPRbbuHa0FhMWsS3tnSFZ2hf07s=
+github.com/fogleman/contourmap v0.0.0-20190814184649-9f61d36c4199/go.mod h1:mqaaaP4j7nTF8T/hx5OCljA7BYWHmrH2uh+Q023OchE=
 github.com/fogleman/gg v1.2.1-0.20190220221249-0403632d5b90/go.mod h1:R/bRT+9gY/C5z7JzPU0zXsXHKM4/ayA+zqcVNZzPa1k=
 github.com/fsnotify/fsnotify v1.4.7 h1:IXs+QLmnXW2CcXuY+8Mzv/fWEsPGWxqefPtCP5CnV9I=
 github.com/fsnotify/fsnotify v1.4.7/go.mod h1:jwhsz4b93w/PPRr/qN1Yymfu8t87LnFCMoQvtojpjFo=