view cmd/isoareas/algo.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
children
line wrap: on
line source

// This is Free Software under GNU Affero General Public License v >= 3.0
// without warranty, see README.md and license for details.
//
// SPDX-License-Identifier: AGPL-3.0-or-later
// License-Filename: LICENSES/AGPL-3.0.txt
//
// Copyright (C) 2019 by via donau
//   – Österreichische Wasserstraßen-Gesellschaft mbH
// Software engineering by Intevation GmbH
//
// Author(s):
//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>

package main

import (
	"bufio"
	"container/list"
	"log"
	"math"
	"math/bits"
	"os"
	"runtime"
	"sync"

	"gemma.intevation.de/gemma/pkg/octree"
)

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
}