# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1569514298 -7200
# Node ID 56f4e8cbfab7bcc9a1c5c29fdd28ee75b570bb94
# Parent  53b55f8116665cd3bf51d4eec4f530de71cf4190
Sew together triangles that totally inside a single class.

diff -r 53b55f811666 -r 56f4e8cbfab7 cmd/isoareas/main.go
--- a/cmd/isoareas/main.go	Thu Sep 26 15:37:56 2019 +0200
+++ b/cmd/isoareas/main.go	Thu Sep 26 18:11:38 2019 +0200
@@ -15,6 +15,7 @@
 
 import (
 	"bufio"
+	"container/list"
 	"fmt"
 	"io"
 	"log"
@@ -142,14 +143,14 @@
 	log.Println("inside classes:")
 	for i, c := range classes {
 
-		var pb polygonBuilder
+		pb := polygonBuilder{open: list.New()}
 
 		var isolated, inside int
 	allInClass:
 		for _, idx := range c {
 			ns := neighbors(tri, idx)
-			w := where(ns, c)
-			switch bits.OnesCount8(w) {
+			mask := where(ns, c)
+			switch bits.OnesCount8(mask) {
 			case 0:
 				// Totally insides do not contribute to the geometries.
 				inside++
@@ -163,11 +164,67 @@
 					tri.Points[ti[1]],
 					tri.Points[ti[2]])
 				continue allInClass
+			default:
+				ti := tri.Triangles[idx*3 : idx*3+3]
+				for j := 0; j < 3; j++ {
+					if mask&(1<<j) == 0 {
+
+						// TODO: Look into cuts to see
+						// if there are real intersections
+						k := (j + 1) % 3
+						front := tri.Points[ti[j]]
+						back := tri.Points[ti[k]]
+
+						curr := octree.LineStringZ{front, back}
+
+					segment:
+						for e := pb.open.Front(); e != nil; e = e.Next() {
+							line := e.Value.(octree.LineStringZ)
+							first, last := line[0], line[len(line)-1]
+
+							if back.EpsEquals(first) {
+								curr = append(curr, line[1:]...)
+								pb.open.Remove(e)
+								goto segment
+							}
+
+							if back.EpsEquals(last) {
+								line.Reverse()
+								curr = append(curr, line[1:]...)
+								pb.open.Remove(e)
+								goto segment
+							}
+
+							if front.EpsEquals(last) {
+								curr = append(line, curr[1:]...)
+								pb.open.Remove(e)
+								goto segment
+							}
+
+							if front.EpsEquals(first) {
+								line.Reverse()
+								curr = append(line, curr[1:]...)
+								pb.open.Remove(e)
+								goto segment
+							}
+						} // all open
+
+						// check if closed
+						if curr[0].EpsEquals(curr[len(curr)-1]) {
+							if !curr.CCW() {
+								curr.Reverse()
+							}
+							pb.polygons = append(pb.polygons, curr)
+						} else {
+							pb.open.PushFront(curr)
+						}
+					} // for all border parts
+				}
 			}
 		}
 
-		log.Printf("\t%d: inside: %d / isolated: %d\n",
-			i, inside, isolated)
+		log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d\n",
+			i, inside, isolated, pb.open.Len(), len(pb.polygons))
 
 		result[i] = pb.polygons
 	}
@@ -185,6 +242,8 @@
 
 type polygonBuilder struct {
 	polygons []octree.LineStringZ
+
+	open *list.List
 }
 
 func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
diff -r 53b55f811666 -r 56f4e8cbfab7 pkg/octree/vertex.go
--- a/pkg/octree/vertex.go	Thu Sep 26 15:37:56 2019 +0200
+++ b/pkg/octree/vertex.go	Thu Sep 26 18:11:38 2019 +0200
@@ -629,6 +629,15 @@
 	return buf.Bytes()
 }
 
+func (ls LineStringZ) CCW() bool {
+	var sum float64
+	for i, v1 := range ls {
+		v2 := ls[(i+1)%len(ls)]
+		sum += (v2.X - v1.X) * (v2.Y + v1.Y)
+	}
+	return sum > 0
+}
+
 // Join joins two lines leaving the first of the second out.
 func (ls LineStringZ) Join(other LineStringZ) LineStringZ {
 	nline := make(LineStringZ, len(ls)+len(other)-1)