diff cmd/octree2contour/loader.go @ 678:7bb961d750b6 octree

octree: traverse horizontally over tree to find out which triangles to process.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 19 Sep 2018 17:34:19 +0200
parents 120a82bd9953
children c79c7be29a7a
line wrap: on
line diff
--- a/cmd/octree2contour/loader.go	Wed Sep 19 15:37:37 2018 +0200
+++ b/cmd/octree2contour/loader.go	Wed Sep 19 17:34:19 2018 +0200
@@ -52,6 +52,8 @@
 		return nil, err
 	}
 
+	log.Printf("EPSG: %d\n", tree.epsg)
+
 	if err := tree.min.read(r); err != nil {
 		return nil, err
 	}
@@ -60,6 +62,10 @@
 		return nil, err
 	}
 
+	log.Printf("BBOX: [[%f, %f, %f], [%f, %f, %f]]\n",
+		tree.min.x, tree.min.y, tree.min.z,
+		tree.max.x, tree.max.y, tree.max.z)
+
 	var numVertices uint32
 	if err := binary.Read(r, binary.LittleEndian, &numVertices); err != nil {
 		return nil, err
@@ -128,6 +134,64 @@
 	return tree, nil
 }
 
+func (ot *octree) horizontal(h float64, fn func([]int32)) {
+
+	type frame struct {
+		pos int32
+		min float64
+		max float64
+	}
+
+	if h < ot.min.z || h > ot.max.z {
+		return
+	}
+
+	stack := []frame{{1, ot.min.z, ot.max.z}}
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		pos, min, max := top.pos, top.min, top.max
+		if pos == 0 {
+			continue
+		}
+
+		if pos > 0 { // node
+			if mid := (max-min)*0.5 + min; h <= mid {
+				stack = append(stack,
+					frame{ot.index[pos+0], min, mid},
+					frame{ot.index[pos+1], min, mid},
+					frame{ot.index[pos+4], min, mid},
+					frame{ot.index[pos+5], min, mid})
+			} else {
+				stack = append(stack,
+					frame{ot.index[pos+2], mid, max},
+					frame{ot.index[pos+3], mid, max},
+					frame{ot.index[pos+6], mid, max},
+					frame{ot.index[pos+7], mid, max})
+			}
+		} else { // leaf
+			pos = -pos - 1
+			n := ot.index[pos]
+			//log.Printf("%d %d %d\n", pos, n, len(ot.index))
+			indices := ot.index[pos+1 : pos+1+n]
+
+			for idx := range indices {
+				tri := ot.triangles[idx]
+				v0 := ot.vertices[tri[0]]
+				v1 := ot.vertices[tri[1]]
+				v2 := ot.vertices[tri[2]]
+
+				if !(math.Min(v0.z, math.Min(v1.z, v2.z)) > h ||
+					math.Max(v0.z, math.Max(v1.z, v2.z)) < h) {
+					fn(tri)
+				}
+			}
+		}
+	}
+}
+
 func loadOctree(fname string) (*octree, error) {
 
 	f, err := os.Open(fname)