# HG changeset patch
# User Sascha L. Teichmann <teichmann@intevation.de>
# Date 1570819750 -7200
# Node ID 5ef04ae34872f5f8d2ba9dbec86f2d98fd866d1c
# Parent  fb09a43b062e7d3005f9ff0034fb7b46106acc3d
Used STRTree for caluculating the diff.

diff -r fb09a43b062e -r 5ef04ae34872 pkg/controllers/diff.go
--- a/pkg/controllers/diff.go	Fri Oct 11 20:09:37 2019 +0200
+++ b/pkg/controllers/diff.go	Fri Oct 11 20:49:10 2019 +0200
@@ -171,7 +171,8 @@
 	}
 
 	// We need a slow path implementation for this.
-	if minuendTree.EPSG != subtrahendTree.EPSG {
+	epsg := minuendTree.EPSG
+	if epsg != subtrahendTree.EPSG {
 		err = mw.JSONError{
 			Code: http.StatusInternalServerError,
 			Message: "Calculating differences between two different " +
@@ -184,6 +185,8 @@
 	points := minuendTree.Diff(subtrahendTree)
 	log.Printf("info: A - B took %v\n", time.Since(start))
 
+	minuendTree, subtrahendTree = nil, nil
+
 	// The Triangulation and the loading of the clipping
 	// polygon can be done concurrently.
 
@@ -217,7 +220,7 @@
 		start := time.Now()
 		clip, clipErr = octree.LoadClippingPolygon(
 			ctx, conn,
-			minuendTree.EPSG,
+			epsg,
 			dci.Bottleneck,
 			dci.Minuend.Time,
 			dci.Subtrahend.Time)
@@ -241,17 +244,17 @@
 	start = time.Now()
 	tin := tri.Tin()
 	removed := tin.Clip(clip)
+	clip = nil
 	log.Printf("info: clipping TIN took %v\n", time.Since(start))
 
+	log.Printf("info: Number of triangles to clip: %d\n", len(removed))
+
 	start = time.Now()
-	log.Printf("info: Number of triangles to clip: %d\n", len(removed))
-	builder := octree.NewBuilder(tin)
-	builder.Build(removed)
-	log.Printf("info: building octree took %v\n", time.Since(start))
+	var tree octree.STRTree
 
-	tree := builder.Tree()
+	tree.BuildWithout(tin, removed)
 
-	log.Printf("info: min/max: %f %f\n", tree.Min.Z, tree.Max.Z)
+	log.Printf("info: Building final STRTree took: %v\n", time.Since(start))
 
 	start = time.Now()
 
@@ -269,10 +272,9 @@
 		"morphology_classbreaks_compare")
 	if err != nil {
 		log.Printf("warn: Loading class breaks failed: %v\n", err)
-		heights = octree.SampleDiffHeights(
-			tree.Min.Z, tree.Max.Z, contourStep)
+		heights = octree.SampleDiffHeights(tin.Min.Z, tin.Max.Z, contourStep)
 	} else {
-		heights = octree.ExtrapolateClassBreaks(heights, tree.Min.Z, tree.Max.Z)
+		heights = octree.ExtrapolateClassBreaks(heights, tin.Min.Z, tin.Max.Z)
 		// heights = octree.InBetweenClassBreaks(heights, 0.05, 2)
 	}
 
@@ -296,7 +298,7 @@
 
 	heights = common.DedupFloat64s(heights)
 
-	areas := octree.TraceAreas(heights, isoCellSize, tree.Min, tree.Max, tree.Value)
+	areas := octree.TraceAreas(heights, isoCellSize, tin.Min, tin.Max, tree.Value)
 
 	var size int
 
@@ -308,7 +310,7 @@
 		size += len(wkb)
 		if _, err = isoStmt.ExecContext(
 			ctx,
-			id, heights[i], minuendTree.EPSG,
+			id, heights[i], epsg,
 			wkb,
 			contourTolerance,
 		); err != nil {
diff -r fb09a43b062e -r 5ef04ae34872 pkg/octree/strtree.go
--- a/pkg/octree/strtree.go	Fri Oct 11 20:09:37 2019 +0200
+++ b/pkg/octree/strtree.go	Fri Oct 11 20:49:10 2019 +0200
@@ -27,6 +27,40 @@
 	bboxes  []Box2D
 }
 
+func (s *STRTree) Value(x, y float64) (float64, bool) {
+
+	stack := []int32{s.index[0]}
+
+	vertices := s.tin.Vertices
+
+	for len(stack) > 0 {
+		top := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		if top > 0 { // node
+			if s.bbox(top).Contains(x, y) {
+				n := s.index[top+1]
+				stack = append(stack, s.index[top+2:top+2+n]...)
+			}
+		} else { // leaf
+			top = -top - 1
+			for i, n := int32(0), s.index[top+1]; i < n; i++ {
+				idx := s.index[top+2+i]
+				ti := s.tin.Triangles[idx]
+				t := Triangle{
+					vertices[ti[0]],
+					vertices[ti[1]],
+					vertices[ti[2]],
+				}
+				if t.Contains(x, y) {
+					return t.Plane3D().Z(x, y), true
+				}
+			}
+		}
+	}
+	return 0, false
+}
+
 func (s *STRTree) Build(t *Tin) {
 
 	if s.Entries == 0 {
@@ -94,6 +128,9 @@
 				s.allTriangles(top, removed)
 			default:
 				// mixed bag
+				// TODO slicing like:
+				// n := s.index[top+1]
+				// stack = append(stack, s.index[top+2:top+2+n]...)
 				for i, n := int32(0), s.index[top+1]; i < n; i++ {
 					stack = append(stack, s.index[top+2+i])
 				}