# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1551719738 -3600
# Node ID e0a7571ac13a348464c48596709b8677a3e57b1f
# Parent  5c3e63cfd50d974c84b736af4739ee0bf6a54d7a
Eliminate bad triangles. Not really the right solution.

diff -r 5c3e63cfd50d -r e0a7571ac13a cmd/octreediff/main.go
--- a/cmd/octreediff/main.go	Mon Mar 04 15:16:39 2019 +0100
+++ b/cmd/octreediff/main.go	Mon Mar 04 18:15:38 2019 +0100
@@ -75,12 +75,11 @@
   WHERE bn.bottleneck_id = $1
 )
 SELECT ST_AsBinary(
-  ST_Buffer(
   ST_intersection(
     (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date),
     (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date)
-  ),
-  0.001)) AS area
+  )
+  ) AS area
 `
 	insertContourSQLClipped = `
 WITH joined AS (
@@ -403,6 +402,7 @@
 
 		tree := builder.Tree()
 		tree.Clip(&clippingPolygon)
+		//tree.Dump()
 
 		now = time.Now()
 		log.Printf("clipping octree took %v\n", now.Sub(last))
diff -r 5c3e63cfd50d -r e0a7571ac13a pkg/octree/tree.go
--- a/pkg/octree/tree.go	Mon Mar 04 15:16:39 2019 +0100
+++ b/pkg/octree/tree.go	Mon Mar 04 18:15:38 2019 +0100
@@ -14,8 +14,10 @@
 package octree
 
 import (
+	"fmt"
 	"log"
 	"math"
+	"os"
 )
 
 // Tree is an Octree holding triangles.
@@ -49,8 +51,156 @@
 	{0.5, 0.5, 1.0, 1.0},
 }
 
+func (ot *Tree) Dump() error {
+	fake, err := os.Create("/tmp/dump.geojson")
+	if err != nil {
+		return err
+	}
+	defer fake.Close()
+
+	fmt.Fprintln(fake,
+		`{
+"crs": {
+    "type": "name",
+    "properties": {
+      "name": "urn:ogc:def:crs:EPSG::32633"
+    }
+  },
+
+  "type": "FeatureCollection",
+  "features": [`)
+
+	first := true
+
+	already := make(map[int32]struct{})
+
+	stack := []int32{1}
+	for len(stack) > 0 {
+		pos := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		if pos > 0 { // node
+			index := ot.index[pos:]
+			if len(index) > 8 {
+				index = index[:8]
+			}
+			for i := range index {
+				if index[i] != 0 {
+					stack = append(stack, index[i])
+				}
+			}
+		} else {
+			pos = -pos - 1
+			n := ot.index[pos]
+			indices := ot.index[pos+1 : pos+1+n]
+			for _, index := range indices {
+				if _, found := already[index]; found {
+					continue
+				}
+				already[index] = struct{}{}
+				tri := ot.triangles[index]
+				t := Triangle{
+					ot.vertices[tri[0]],
+					ot.vertices[tri[1]],
+					ot.vertices[tri[2]],
+				}
+				if first {
+					first = false
+				} else {
+					fmt.Fprintln(fake, ",")
+				}
+				fmt.Fprintf(fake,
+					`{ "type": "Feature",
+  "properties": {
+     "id": %d
+  },
+  "geometry": {
+     "type": "Polygon",
+     "coordinates": [[
+       [ %.6f, %.6f ],
+       [ %.6f, %.6f ],
+       [ %.6f, %.6f ],
+       [ %.6f, %.6f ]]
+    ]
+  }
+}`, index,
+					t[0].X, t[0].Y,
+					t[1].X, t[1].Y,
+					t[2].X, t[2].Y,
+					t[0].X, t[0].Y)
+			}
+		}
+	}
+
+	fmt.Fprintln(fake,
+		`]
+}`)
+	return nil
+}
+
+func (ot *Tree) eliminate(removed map[int32]IntersectionType) {
+
+	var bad int
+
+	stack := []int32{1}
+	for len(stack) > 0 {
+		pos := stack[len(stack)-1]
+		stack = stack[:len(stack)-1]
+
+		if pos > 0 { // node
+			index := ot.index[pos:]
+			if len(index) > 8 {
+				index = index[:8]
+			}
+			for i := range index {
+				if index[i] != 0 {
+					stack = append(stack, index[i])
+				}
+			}
+		} else {
+			pos = -pos - 1
+			n := ot.index[pos]
+			indices := ot.index[pos+1 : pos+1+n]
+
+			for i := len(indices) - 1; i >= 0; i-- {
+				index := indices[i]
+				if t, found := removed[index]; found && t != IntersectionInside {
+					if i < len(indices)-1 {
+						copy(indices[i:], indices[i+1:])
+					}
+					indices[len(indices)-1] = 0
+					indices = indices[:len(indices)-1]
+					bad++
+				}
+			}
+			ot.index[pos] = int32(len(indices))
+		}
+	}
+
+	log.Printf("bad triangles: %d\n", bad)
+}
+
 func (ot *Tree) Clip(p *Polygon) {
 
+	/*
+			fake, _ := os.Create("/tmp/removed.geojson")
+			defer fake.Close()
+
+			fmt.Fprintln(fake,
+				`{
+		"crs": {
+		    "type": "name",
+		    "properties": {
+		      "name": "urn:ogc:def:crs:EPSG::32633"
+		    }
+		  },
+
+		  "type": "FeatureCollection",
+		  "features": [`)
+
+			first := true
+	*/
+
 	log.Printf("info: num triangles: %d\n", len(ot.triangles))
 
 	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
@@ -59,7 +209,6 @@
 
 	triChecks := make(map[int32]IntersectionType)
 
-	var triangleTests int
 	var nodeTests int
 	var nodesClipped int
 	var trianglesClipped int
@@ -151,7 +300,36 @@
 					}
 					what = p.IntersectionWithTriangle(&t)
 					triChecks[triIndex] = what
-					triangleTests++
+
+					/*
+
+											if what != IntersectionInside {
+												if first {
+													first = false
+												} else {
+													fmt.Fprintln(fake, ",")
+												}
+												fmt.Fprintf(fake,
+													`{ "type": "Feature",
+						  "properties": {
+						    "id": %d
+						  },
+						  "geometry": {
+						     "type": "Polygon",
+						     "coordinates": [[
+						       [ %.6f, %.6f ],
+						       [ %.6f, %.6f ],
+						       [ %.6f, %.6f ],
+						       [ %.6f, %.6f ]]
+						    ]
+						  }
+						}`, triIndex,
+													t[0].X, t[0].Y,
+													t[1].X, t[1].Y,
+													t[2].X, t[2].Y,
+													t[0].X, t[0].Y)
+											}
+					*/
 				}
 				switch what {
 				case IntersectionInside:
@@ -170,11 +348,18 @@
 			ot.index[pos] = int32(len(indices))
 		}
 	}
+	/*
+			fmt.Fprintln(fake,
+				`]
+		}`)
+	*/
 	log.Printf("info: node tests: %d\n", nodeTests)
 	log.Printf("info: nodes clipped: %d\n", nodesClipped)
 	log.Printf("info: nodes all inside: %d\n", nodesAllInside)
-	log.Printf("info: triangle tests: %d\n", triangleTests)
 	log.Printf("info: triangle clipped: %d\n", trianglesClipped)
+	log.Printf("info: tested triangles: %d\n", len(triChecks))
+
+	ot.eliminate(triChecks)
 }
 
 func (ot *Tree) Value(x, y float64) (float64, bool) {