# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1563205790 -7200
# Node ID 4233570de212033588052b1fd1246af47bd04595
# Parent  858b1bf1d21b583d9aede57ebe7941af8129538b
Building octrees: Don't cause havok by building ain overly large tree if the points are not really 3D.

diff -r 858b1bf1d21b -r 4233570de212 pkg/octree/builder.go
--- a/pkg/octree/builder.go	Mon Jul 15 17:09:59 2019 +0200
+++ b/pkg/octree/builder.go	Mon Jul 15 17:49:50 2019 +0200
@@ -37,7 +37,7 @@
 
 type buildStep func(chan buildStep)
 
-var cubes = [8][2]Vertex{
+var cubes = [8]Box{
 	makeCube(0),
 	makeCube(1),
 	makeCube(2),
@@ -48,7 +48,7 @@
 	makeCube(7),
 }
 
-func makeCube(i int) [2]Vertex {
+func makeCube(i int) Box {
 	var d Vertex
 	if i&1 == 1 {
 		d.X = 0.5
@@ -59,12 +59,19 @@
 	if i&4 == 4 {
 		d.Z = 0.5
 	}
-	return [2]Vertex{
+	return Box{
 		Vertex{0.0, 0.0, 0.0}.Add(d),
 		Vertex{0.5, 0.5, 0.5}.Add(d),
 	}
 }
 
+func twoElseOne(b bool) int {
+	if b {
+		return 2
+	}
+	return 1
+}
+
 // NewBuilder creates a new Builder for a TIN.
 func NewBuilder(t *Tin) *Builder {
 	return &Builder{t: t}
@@ -143,15 +150,28 @@
 			parent(tb.buildRecursive(triangles, min, max, depth))
 			return
 		}
+		box := Box{min, max}
 
-		bbox := Interpolate(min, max)
+		xLimit := twoElseOne(box.HasX())
+		yLimit := twoElseOne(box.HasY())
+		zLimit := twoElseOne(box.HasZ())
 
-		bboxes := make([][2]Vertex, len(cubes))
+		indices := make([]byte, 0, 8)
+
+		bbox := box.Interpolate()
 
-		for i := range cubes {
-			bboxes[i] = [2]Vertex{
-				bbox(cubes[i][0]),
-				bbox(cubes[i][1]),
+		var bboxes [8]Box
+
+		for x := 0; x < xLimit; x++ {
+			for y := 0; y < yLimit; y++ {
+				for z := 0; z < zLimit; z++ {
+					idx := byte(z<<2 | y<<1 | x)
+					bboxes[idx] = Box{
+						bbox(cubes[idx][0]),
+						bbox(cubes[idx][1]),
+					}
+					indices = append(indices, idx)
+				}
 			}
 		}
 
@@ -171,7 +191,7 @@
 			h.Maximize(v1)
 			h.Maximize(v2)
 
-			for i := range bboxes {
+			for _, i := range indices {
 				if !(h.Less(bboxes[i][0]) || bboxes[i][1].Less(l)) {
 					quandrants[i] = append(quandrants[i], tri)
 				}
@@ -179,7 +199,7 @@
 		}
 
 		used := new(int32)
-		for i := range quandrants {
+		for _, i := range indices {
 			if len(quandrants[i]) > 0 {
 				*used++
 			}
@@ -187,7 +207,12 @@
 
 		pos := tb.allocNode()
 
-		for i := range quandrants {
+		if *used == 0 {
+			parent(pos)
+			return
+		}
+
+		for _, i := range indices {
 			if len(quandrants[i]) > 0 {
 				j := int32(i)
 				parent := func(v int32) {
@@ -239,14 +264,28 @@
 		return int32(-(pos + 1))
 	}
 
-	bbox := Interpolate(min, max)
+	box := Box{min, max}
 
-	bboxes := make([][2]Vertex, len(cubes))
+	xLimit := twoElseOne(box.HasX())
+	yLimit := twoElseOne(box.HasY())
+	zLimit := twoElseOne(box.HasZ())
+
+	indices := make([]byte, 0, 8)
+
+	bbox := box.Interpolate()
 
-	for i := range cubes {
-		bboxes[i] = [2]Vertex{
-			bbox(cubes[i][0]),
-			bbox(cubes[i][1]),
+	var bboxes [8]Box
+
+	for x := 0; x < xLimit; x++ {
+		for y := 0; y < yLimit; y++ {
+			for z := 0; z < zLimit; z++ {
+				idx := byte(z<<2 | y<<1 | x)
+				bboxes[idx] = Box{
+					bbox(cubes[idx][0]),
+					bbox(cubes[idx][1]),
+				}
+				indices = append(indices, idx)
+			}
 		}
 	}
 
@@ -266,7 +305,7 @@
 		h.Maximize(v1)
 		h.Maximize(v2)
 
-		for i := range bboxes {
+		for _, i := range indices {
 			if !(h.Less(bboxes[i][0]) || bboxes[i][1].Less(l)) {
 				quandrants[i] = append(quandrants[i], tri)
 			}
@@ -275,7 +314,7 @@
 
 	pos := tb.allocNode()
 
-	for i := range quandrants {
+	for _, i := range indices {
 		if len(quandrants[i]) > 0 {
 			child := tb.buildRecursive(
 				quandrants[i],
diff -r 858b1bf1d21b -r 4233570de212 pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Mon Jul 15 17:09:59 2019 +0200
+++ b/pkg/octree/vertex.go	Mon Jul 15 17:49:50 2019 +0200
@@ -45,6 +45,9 @@
 	// and the second being the direction.
 	Line [2]Vertex
 
+	// Box is a 3D box.
+	Box [2]Vertex
+
 	// MultiPointZ is a set of vertices.
 	MultiPointZ []Vertex
 
@@ -288,9 +291,10 @@
 	}
 }
 
-// Interpolate returns a function that return s*v2 + v1
+// Interpolate returns a function that return s*b[1] + b[0]
 // component-wise.
-func Interpolate(v1, v2 Vertex) func(Vertex) Vertex {
+func (b Box) Interpolate() func(Vertex) Vertex {
+	v1, v2 := b[0], b[1]
 	v2 = v2.Sub(v1)
 	return func(s Vertex) Vertex {
 		return Vertex{
@@ -301,6 +305,10 @@
 	}
 }
 
+func (b Box) HasX() bool { return math.Abs(b[0].X-b[1].X) > epsPlane }
+func (b Box) HasY() bool { return math.Abs(b[0].Y-b[1].Y) > epsPlane }
+func (b Box) HasZ() bool { return math.Abs(b[0].Z-b[1].Z) > epsPlane }
+
 // Less returns if one of v component is less than the
 // corresponing component in w.
 func (v Vertex) Less(w Vertex) bool {