# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1537840545 -7200
# Node ID 5e14000829d152e9d34c8fa9a6f09682b5dafa7d
# Parent  1a1a8b5f2d023d75928d17e8aa733502d2e96fc4
Complete vertical octree traversal. Seems not very selective.

diff -r 1a1a8b5f2d02 -r 5e14000829d1 pkg/controllers/octreecross.go
--- a/pkg/controllers/octreecross.go	Mon Sep 24 18:41:29 2018 +0200
+++ b/pkg/controllers/octreecross.go	Tue Sep 25 03:55:45 2018 +0200
@@ -79,7 +79,14 @@
 	}
 
 	for i := 0; i < len(coords)-1; i++ {
-		// TODO: Query segmentwise (i, i+1).
+		var count int
+		c1 := &coords[i]
+		c2 := &coords[i+1]
+		tree.Vertical(c1.Lat, c1.Lon, c2.Lat, c2.Lon, func(*octree.Triangle) {
+			// TODO: Do real intersection.
+			count++
+		})
+		log.Printf("triangles to intersect: %d\n", count)
 	}
 
 	return
diff -r 1a1a8b5f2d02 -r 5e14000829d1 pkg/octree/tree.go
--- a/pkg/octree/tree.go	Mon Sep 24 18:41:29 2018 +0200
+++ b/pkg/octree/tree.go	Tue Sep 25 03:55:45 2018 +0200
@@ -1,6 +1,7 @@
 package octree
 
 import (
+	"log"
 	"math"
 )
 
@@ -45,13 +46,36 @@
 	return a.y2
 }
 
-var scale = [][]float64{
+var scale = [4][4]float64{
 	{0.0, 0.0, 0.5, 0.5},
 	{0.5, 0.0, 1.0, 0.5},
 	{0.0, 0.5, 0.5, 1.0},
 	{0.5, 0.5, 1.0, 1.0},
 }
 
+type plane2d struct {
+	a float64
+	b float64
+	c float64
+}
+
+func newPlane2d(x1, y1, x2, y2 float64) plane2d {
+	a := x2 - x1
+	b := y2 - y1
+
+	l := math.Sqrt(a*a + b*b)
+	a /= l
+	b /= l
+
+	// a*x1 + b*y1 + c = 0
+	c := -(a*x1 + b*y1)
+	return plane2d{a, b, c}
+}
+
+func (p plane2d) eval(x, y float64) float64 {
+	return p.a*x + p.b*y + p.c
+}
+
 func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {
 
 	box := box2d{
@@ -67,13 +91,19 @@
 		return
 	}
 
+	line := newPlane2d(x1, y1, x2, y2)
+
 	type frame struct {
 		pos int32
 		box2d
 	}
 
+	dupes := map[int32]struct{}{}
+
 	stack := []frame{{1, box2d{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}}}
 
+	var lineRejected int
+
 	for len(stack) > 0 {
 		top := stack[len(stack)-1]
 		stack = stack[:len(stack)-1]
@@ -83,13 +113,13 @@
 
 			var all uint
 			for i := 0; i < 2; i++ {
+				var xcode uint
+				if box.xi(i) > midx {
+					xcode = 1
+				}
 				for j := 0; j < 2; j++ {
-					x, y := box.xi(i), box.yi(j)
-					var code uint
-					if x > midx {
-						code |= 1
-					}
-					if y > midy {
+					code := xcode
+					if box.yi(j) > midy {
 						code |= 2
 					}
 					all |= 1 << code
@@ -98,24 +128,73 @@
 
 			dx := top.x2 - top.x1
 			dy := top.y2 - top.y1
-			for i := 0; all != 0; i++ {
+			for i := int32(0); all != 0; i++ {
 				if all&1 == 1 {
-					s := scale[i]
-					nbox := box2d{
-						dx*s[0] + top.x1,
-						dy*s[1] + top.y1,
-						dx*s[2] + top.x1,
-						dy*s[3] + top.y1,
+					a := ot.index[top.pos+i]
+					b := ot.index[top.pos+i+4]
+					if b != 0 || b != 0 {
+						nbox := box2d{
+							dx*scale[i][0] + top.x1,
+							dy*scale[i][1] + top.y1,
+							dx*scale[i][2] + top.x1,
+							dy*scale[i][3] + top.y1,
+						}
+						if a != 0 {
+							stack = append(stack, frame{a, nbox})
+						}
+						if b != 0 {
+							stack = append(stack, frame{b, nbox})
+						}
 					}
-					// TODO: Push two frames to stack.
-					_ = nbox
 				}
 				all >>= 1
 			}
+		} else { // leaf
+			pos := -top.pos - 1
+			n := ot.index[pos]
+			indices := ot.index[pos+1 : pos+1+n]
 
-		} else { // leaf
+			for _, idx := range indices {
+				if _, found := dupes[idx]; found {
+					continue
+				}
+				tri := ot.triangles[idx]
+				t := Triangle{
+					ot.vertices[tri[0]],
+					ot.vertices[tri[1]],
+					ot.vertices[tri[2]],
+				}
+
+				v0 := line.eval(t[0].X, t[0].Y)
+				v1 := line.eval(t[1].X, t[1].Y)
+				v2 := line.eval(t[2].X, t[2].Y)
+
+				sides := func(s int, b bool) int {
+					if b {
+						return s | 2
+					}
+					return s | 1
+				}
+
+				const eps = 1e-05
+
+				if math.Abs(v0) < eps || math.Abs(v1) < eps ||
+					math.Abs(v2) < eps ||
+					sides(sides(sides(0,
+						math.Signbit(v0)),
+						math.Signbit(v1)),
+						math.Signbit(v2)) == 3 {
+					dupes[idx] = struct{}{}
+					fn(&t)
+				} else {
+					lineRejected++
+				}
+			}
 		}
 	}
+
+	// XXX: This value seems to high!
+	log.Printf("rejected by line test: %d\n", lineRejected)
 }
 
 func (ot *Tree) Horizontal(h float64, fn func(*Triangle)) {
@@ -130,9 +209,9 @@
 		max float64
 	}
 
-	stack := []frame{{1, ot.Min.Z, ot.Max.Z}}
+	dupes := map[int32]struct{}{}
 
-	dupes := map[int32]struct{}{}
+	stack := []frame{{1, ot.Min.Z, ot.Max.Z}}
 
 	for len(stack) > 0 {
 		top := stack[len(stack)-1]