--- a/pkg/controllers/diff.go	Thu Oct 17 16:36:58 2019 +0200
+++ b/pkg/controllers/diff.go	Thu Oct 17 18:33:46 2019 +0200
@@ -84,15 +84,8 @@
   $2,
   ST_Transform(
     ST_Multi(
-      ST_CollectionExtract(
-        ST_MakeValid(
-          ST_Multi(
-             ST_Collectionextract(
-                ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5),
-                3
-             )
-            )
-        ),
+      ST_Collectionextract(
+        ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5),
         3
       )
     ),

--- a/pkg/octree/areas.go	Thu Oct 17 16:36:58 2019 +0200
+++ b/pkg/octree/areas.go	Thu Oct 17 18:33:46 2019 +0200
@@ -189,6 +189,10 @@
 	reprojX := common.Linear(0.5, min.X, 1.5, min.X+cellSize)
 	reprojY := common.Linear(0.5, min.Y, 1.5, min.Y+cellSize)
 
+	reproj := func(x, y float64) (float64, float64) {
+		return reprojX(x), reprojY(y)
+	}
+
 	cnts := make(chan int)
 
 	doContours := func() {
@@ -200,28 +204,7 @@
 				continue
 			}
 
-			// We need to bring it back to the
-			// none raster coordinate system.
-			a := make(wkb.MultiPolygonGeom, len(c))
-
-			for i, pl := range c {
-				shell := make(wkb.LinearRingGeom, len(pl))
-				for j, pt := range pl {
-					shell[j] = wkb.PointGeom{
-						X: reprojX(pt.X),
-						Y: reprojY(pt.Y),
-					}
-				}
-				/*
-					if !shell.CCW() {
-						log.Println("not ccw")
-						shell.Reverse()
-					}
-				*/
-				a[i] = wkb.PolygonGeom{shell}
-			}
-
-			areas[hIdx] = a
+			areas[hIdx] = buildMultipolygon(c, reproj)
 		}
 	}
 
@@ -240,3 +223,197 @@
 
 	return areas
 }
+
+func contourBBox(cnt contourmap.Contour) 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}
+}
+
+type bboxNode struct {
+	box      Box2D
+	cnt      contourmap.Contour
+	children []*bboxNode
+}
+
+type bboxForest struct {
+	roots []*bboxNode
+}
+
+func (bn *bboxNode) insert(cnt contourmap.Contour, box Box2D) {
+	// check if roots are inside new
+	var nr *bboxNode
+
+	for i, r := range bn.children {
+		if r.box.Inside(box) {
+			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
+	}
+
+	var inserted bool
+
+	// check if new is inside an old
+	for _, r := range bn.children {
+		if box.Inside(r.box) {
+			r.insert(cnt, box)
+			inserted = true
+		}
+	}
+
+	// its a new root node.
+	if !inserted {
+		nr = &bboxNode{box: box, cnt: cnt}
+		bn.children = append(bn.children, nr)
+	}
+}
+
+func (bf *bboxForest) insert(cnt contourmap.Contour) {
+	box := contourBBox(cnt)
+
+	// check if roots are inside new
+	var nr *bboxNode
+
+	for i, r := range bf.roots {
+		if r.box.Inside(box) {
+			if nr == nil {
+				nr = &bboxNode{box: box, cnt: cnt}
+			}
+			nr.children = append(nr.children, r)
+			bf.roots[i] = nil
+		}
+	}
+
+	// we have a new root
+	if nr != nil {
+		// compact the list
+		for i := len(bf.roots) - 1; i >= 0; i-- {
+			if bf.roots[i] == nil {
+				if i < len(bf.roots)-1 {
+					copy(bf.roots[i:], bf.roots[i+1:])
+				}
+				bf.roots[len(bf.roots)-1] = nil
+				bf.roots = bf.roots[:len(bf.roots)-1]
+			}
+		}
+		bf.roots = append(bf.roots, nr)
+		return
+	}
+
+	var inserted bool
+
+	// check if new is inside an old
+	for _, r := range bf.roots {
+		if box.Inside(r.box) {
+			r.insert(cnt, box)
+			inserted = true
+		}
+	}
+
+	// its a new root node.
+	if !inserted {
+		nr = &bboxNode{box: box, cnt: cnt}
+		bf.roots = append(bf.roots, nr)
+	}
+}
+
+type bboxOutFunc func(contourmap.Contour, []contourmap.Contour)
+
+func (bn *bboxNode) generate(shell bool, out bboxOutFunc) {
+	if shell && len(bn.children) == 0 {
+		out(bn.cnt, nil)
+		return
+	}
+}
+
+func (bf *bboxForest) generate(out bboxOutFunc) {
+	for _, r := range bf.roots {
+		r.generate(true, out)
+	}
+}
+
+func buildMultipolygon(
+	cnts []contourmap.Contour,
+	reproj func(float64, float64) (float64, float64),
+) wkb.MultiPolygonGeom {
+
+	var bf bboxForest
+
+	for _, cnt := range cnts {
+		bf.insert(cnt)
+	}
+
+	log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots))
+
+	var mp wkb.MultiPolygonGeom
+
+	out := func(sh contourmap.Contour, hls []contourmap.Contour) {
+
+		polygon := make(wkb.PolygonGeom, 1+len(hls))
+
+		// Handle shell
+		shell := make(wkb.LinearRingGeom, len(sh))
+		for j, pt := range sh {
+			y, x := reproj(pt.X, pt.Y)
+			shell[j] = wkb.PointGeom{X: x, Y: 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 {
+				y, x := reproj(pt.X, pt.Y)
+				hole[j] = wkb.PointGeom{X: x, Y: y}
+			}
+			if hole.CCW() {
+				hole.Reverse()
+			}
+			polygon[1+i] = hole
+		}
+
+		mp = append(mp, polygon)
+	}
+
+	bf.generate(out)
+	log.Printf("generated: %d\n", len(mp))
+
+	return mp
+}

--- a/pkg/octree/vertex.go	Thu Oct 17 16:36:58 2019 +0200
+++ b/pkg/octree/vertex.go	Thu Oct 17 18:33:46 2019 +0200
@@ -108,6 +108,11 @@
 	}
 }
 
+func (a Box2D) Inside(b Box2D) bool {
+	return a.X1 >= b.X1 && a.X2 <= b.X2 &&
+		a.Y1 >= b.Y1 && a.Y2 <= b.Y2
+}
+
 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