diff cmd/isoareas/main.go @ 4551:b5b23b6d76f1 iso-areas

Move own algorith to separate file.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 01 Oct 2019 11:07:33 +0200
parents aa2d0006e742
children 7ed5a4a94efc
line wrap: on
line diff
--- a/cmd/isoareas/main.go	Mon Sep 30 16:34:55 2019 +0200
+++ b/cmd/isoareas/main.go	Tue Oct 01 11:07:33 2019 +0200
@@ -15,18 +15,14 @@
 
 import (
 	"bufio"
-	"container/list"
 	"fmt"
 	"io"
 	"log"
 	"math"
-	"math/bits"
 	"math/rand"
 	"os"
-	"runtime"
 	"strconv"
 	"strings"
-	"sync"
 	"time"
 
 	svg "github.com/ajstarks/svgo"
@@ -91,398 +87,6 @@
 	}
 }
 
-type indexedArc struct {
-	arc   int32
-	index int32
-}
-
-func glue(a, b octree.LineStringZ) octree.LineStringZ {
-
-	a1, a2 := a[0], a[len(a)-1]
-	b1, b2 := b[0], b[len(b)-1]
-
-	if a1.EpsEquals2D(b2) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, b)
-		copy(c[len(b):], a[1:])
-		return c
-	}
-
-	if a2.EpsEquals2D(b1) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, a)
-		copy(c[len(a):], b[1:])
-		return c
-	}
-
-	if a1.EpsEquals2D(b1) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, b)
-		c[:len(b)].Reverse()
-		copy(c[len(b):], a[1:])
-		return c
-	}
-
-	if a2.EpsEquals2D(b2) {
-		c := make(octree.LineStringZ, len(a)-1+len(b))
-		copy(c, a)
-		copy(c[len(a):], b[:len(b)-1])
-		c[len(a):].Reverse()
-		return c
-	}
-
-	return nil
-}
-
-func connectArcs(tri *octree.Triangulation, cuts []indexedArc, arcs *[]octree.LineStringZ) {
-
-	unique := map[int32]struct{}{}
-	for i := range cuts {
-		unique[cuts[i].arc] = struct{}{}
-	}
-	before := len(unique)
-
-	origLen := int32(len(*arcs))
-
-	for i := range cuts {
-		if cuts[i].arc >= origLen {
-			// already has a connected arc assigned.
-			continue
-		}
-
-		line := (*arcs)[cuts[i].arc]
-
-		var modified []int
-		visited := map[int32]struct{}{}
-
-		var recursive func(int32)
-		recursive = func(idx int32) {
-			visited[idx] = struct{}{}
-
-			ns := neighbors(tri, idx)
-			for _, n := range ns {
-				n /= 3
-				if _, already := visited[n]; already {
-					continue
-				}
-
-				arcIdx := findArc(n, cuts)
-				if arcIdx == -1 {
-					continue
-				}
-				arc := cuts[arcIdx].arc
-				if arc >= origLen {
-					// already a new arc.
-					continue
-				}
-				oline := (*arcs)[arc]
-
-				nline := glue(line, oline)
-				if nline == nil {
-					// not joinable
-					continue
-				}
-				line = nline
-				modified = append(modified, arcIdx)
-				recursive(cuts[arcIdx].index)
-			}
-		}
-		recursive(cuts[i].index)
-		if len(modified) > 0 {
-			// alloc a new line an fix the references.
-			nidx := int32(len(*arcs))
-			*arcs = append(*arcs, line)
-			cuts[i].arc = nidx
-			for _, j := range modified {
-				cuts[j].arc = nidx
-			}
-		}
-	}
-
-	unique = map[int32]struct{}{}
-	for i := range cuts {
-		unique[cuts[i].arc] = struct{}{}
-	}
-	log.Printf("unique arcs: before: %d after %d\n",
-		before, len(unique))
-}
-
-func intersectTriangles(
-	tri *octree.Triangulation,
-	heights []float64,
-	min, max octree.Vertex,
-) [][]octree.LineStringZ {
-
-	cuts := make([][]indexedArc, len(heights))
-
-	classes := make([][]int32, len(heights)+1)
-
-	var arcs []octree.LineStringZ
-
-all:
-	for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
-		ti := tri.Triangles[i : i+3]
-		v0 := tri.Points[ti[0]]
-		v1 := tri.Points[ti[1]]
-		v2 := tri.Points[ti[2]]
-		min := math.Min(v0.Z, math.Min(v1.Z, v2.Z))
-
-		if min > heights[len(heights)-1] {
-			classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
-			continue
-		}
-		max := math.Max(v0.Z, math.Max(v1.Z, v2.Z))
-
-		for j, h := range heights {
-
-			var l float64
-			if j > 0 {
-				l = heights[j-1]
-			} else {
-				l = -math.MaxFloat64
-			}
-
-			if l < min && max < h {
-				classes[j] = append(classes[j], i/3)
-				continue all
-			}
-			if min > h || max < h {
-				continue
-			}
-			t := octree.Triangle{v0, v1, v2}
-			cut := t.IntersectHorizontal(h)
-			if len(cut) >= 2 {
-				arc := int32(len(arcs))
-				arcs = append(arcs, cut)
-				cuts[j] = append(cuts[j], indexedArc{arc: arc, index: i / 3})
-			}
-		}
-	}
-
-	// connect the arcs in a cut list to longer arcs.
-
-	for _, c := range cuts {
-		connectArcs(tri, c, &arcs)
-	}
-
-	func() {
-		out, err := os.Create("arcs.svg")
-		if err != nil {
-			log.Printf("err: %v\n", err)
-			return
-		}
-		defer func() {
-			out.Close()
-			log.Println("writing arcs done")
-		}()
-
-		buf := bufio.NewWriter(out)
-
-		dumpArcsSVG(
-			buf,
-			min, max,
-			cuts,
-			arcs)
-		buf.Flush()
-	}()
-
-	result := make([][]octree.LineStringZ, len(classes))
-
-	jobs := make(chan int)
-
-	var wg sync.WaitGroup
-
-	worker := func() {
-		defer wg.Done()
-
-		pb := polygonBuilder{open: list.New()}
-
-		for i := range jobs {
-			usedArcs := map[int32]struct{}{}
-			var dupes int
-			var isolated, inside, found int
-			c := classes[i]
-
-		allInClass:
-			for _, idx := range c {
-				ns := neighbors(tri, idx)
-				mask := where(ns, c)
-				switch bits.OnesCount8(mask) {
-				case 3:
-					// Totally insides do not contribute to the geometries.
-					inside++
-					continue allInClass
-				case 0:
-					// Isolated are areas w/o connections to the rest.
-					isolated++
-					ti := tri.Triangles[idx*3 : idx*3+3]
-					pb.addTriangle(
-						tri.Points[ti[0]],
-						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 {
-
-							var curr octree.LineStringZ
-
-							for l := i - 1; l <= i; l++ {
-								if l < 0 || l >= len(cuts) {
-									continue
-								}
-								arcIdx := findArc(ns[j]/3, cuts[l])
-								if arcIdx == -1 {
-									continue
-								}
-								found++
-								aIdx := cuts[l][arcIdx].arc
-								if _, already := usedArcs[aIdx]; already {
-									dupes++
-									continue
-								}
-								usedArcs[aIdx] = struct{}{}
-
-								curr = arcs[aIdx]
-								break
-							}
-
-							if curr == nil {
-								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)
-								if nline := glue(curr, line); nline != nil {
-									curr = nline
-									pb.open.Remove(e)
-									goto segment
-								}
-							} // all open
-
-							// check if closed
-							if curr[0].EpsEquals2D(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 open: %d closed: %d dupes: %d found: %d\n",
-				i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found)
-			/*
-				for e := pb.open.Front(); e != nil; e = e.Next() {
-					line := e.Value.(octree.LineStringZ)
-					if !line.CCW() {
-						line.Reverse()
-					}
-					pb.polygons = append(pb.polygons, line)
-				}
-			*/
-
-			result[i] = pb.polygons
-			pb.reset()
-		}
-	}
-
-	for n := runtime.NumCPU(); n >= 0; n-- {
-		wg.Add(1)
-		go worker()
-	}
-
-	log.Println("inside classes:")
-	for i := range classes {
-		jobs <- i
-	}
-	close(jobs)
-
-	wg.Wait()
-
-	log.Println("cuts:")
-	for i, c := range cuts {
-		log.Printf("\t%.3f: %d\n", heights[i], len(c))
-	}
-
-	return result
-
-	// TODO: sew segments together.
-
-}
-
-type polygonBuilder struct {
-	polygons []octree.LineStringZ
-
-	open *list.List
-}
-
-func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
-	polygon := octree.LineStringZ{v0, v1, v2, v0}
-	pb.polygons = append(pb.polygons, polygon)
-}
-
-func (pb *polygonBuilder) reset() {
-	pb.polygons = pb.polygons[:0]
-	pb.open.Init()
-}
-
-func neighbors(t *octree.Triangulation, idx int32) []int32 {
-	idx *= 3
-	return t.Halfedges[idx : idx+3]
-}
-
-func findArc(needle int32, haystack []indexedArc) int {
-	lo, hi := 0, len(haystack)-1
-	for lo <= hi {
-		mid := (hi-lo)/2 + lo
-		switch v := haystack[mid].index; {
-		case v < needle:
-			lo = mid + 1
-		case v > needle:
-			hi = mid - 1
-		default:
-			return mid
-		}
-	}
-	return -1
-}
-
-func contains(needle int32, haystack []int32) bool {
-	lo, hi := 0, len(haystack)-1
-	for lo <= hi {
-		mid := (hi-lo)/2 + lo
-		switch v := haystack[mid]; {
-		case v < needle:
-			lo = mid + 1
-		case v > needle:
-			hi = mid - 1
-		default:
-			return true
-		}
-	}
-	return false
-}
-
-func where(neighbors, indices []int32) byte {
-	var mask byte
-	for i, n := range neighbors {
-		if n >= 0 && contains(n/3, indices) {
-			mask |= 1 << i
-		}
-	}
-	return mask
-}
-
 func linear(x1, y1, x2, y2 float64) func(float64) float64 {
 	// f(x1) = y1
 	// f(x2) = y2