view pkg/octree/strtree.go @ 4017:639bdb17c3f2

Fixed offset for fairway box This was broken by changeset: 4080:bf86f9a08733 user: Thomas Junk <thomas.junk@intevation.de> Date: Thu Jul 18 15:04:30 2019 +0200 summary: improve fairwaydiagram printing positioning For the record: I think the current implementation exceptionally flawed. Instead of adding extra offset parameters to the diagram elements the whole building block with all contained elements should be translated in one step, that would be less cluttered and less error prone...
author Sascha Wilde <wilde@intevation.de>
date Fri, 19 Jul 2019 16:59:25 +0200
parents 114979e97a6c
children f456ce0a6a0e
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) 2018 by via donau
//   – Österreichische Wasserstraßen-Gesellschaft mbH
// Software engineering by Intevation GmbH
//
// Author(s):
//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>

package octree

import (
	"math"
	"sort"
)

const numEntries = 8

type STRTree struct {
	tin    *Tin
	index  []int32
	bboxes []Box2D
}

func (s *STRTree) Build(t *Tin) {

	s.tin = t

	all := make([]int32, len(t.Triangles))

	for i := range all {
		all[i] = int32(i)
	}

	s.index = append(s.index, 0)

	root := s.build(all)

	s.index[0] = root
}

func (s *STRTree) Clip(p *Polygon) map[int32]struct{} {

	removed := make(map[int32]struct{})

	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
			switch p.IntersectionBox2D(s.bbox(top)) {
			case IntersectionInside:
				// all triangles are inside polygon
			case IntersectionOutSide:
				// all triangles are outside polygon
				s.allTriangles(top, removed)
			default:
				// mixed bag
				for i, n := int32(0), s.index[top+1]; i < n; i++ {
					stack = append(stack, s.index[top+2+i])
				}
			}
		} 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 p.IntersectionWithTriangle(&t) != IntersectionInside {
					removed[idx] = struct{}{}
				}
			}
		}
	}

	return removed
}

func (s *STRTree) allTriangles(pos int32, tris map[int32]struct{}) {
	stack := []int32{pos}
	for len(stack) > 0 {
		top := stack[len(stack)-1]
		stack = stack[:len(stack)-1]
		if top > 0 { // node
			for i, n := int32(0), s.index[top+1]; i < n; i++ {
				stack = append(stack, s.index[top+2+i])
			}
		} else { // leaf
			top = -top - 1
			for i, n := int32(0), s.index[top+1]; i < n; i++ {
				tris[s.index[top+2+i]] = struct{}{}
			}
		}
	}
}

func (s *STRTree) build(items []int32) int32 {
	sort.Slice(items, func(i, j int) bool {
		ti := s.tin.Triangles[items[i]]
		tj := s.tin.Triangles[items[j]]
		return s.tin.Vertices[ti[0]].X < s.tin.Vertices[tj[0]].X
	})

	P := int(math.Ceil(float64(len(items)) / float64(numEntries)))
	S := int(math.Ceil(math.Sqrt(float64(P))))

	slices := strSplit(items, S)

	leaves := strJoin(
		slices, S,
		func(i, j int32) bool {
			ti := s.tin.Triangles[i]
			tj := s.tin.Triangles[j]
			return s.tin.Vertices[ti[0]].Y < s.tin.Vertices[tj[0]].Y
		},
		s.allocLeaf,
	)

	return s.buildNodes(leaves)
}

func (s *STRTree) buildNodes(items []int32) int32 {

	if len(items) <= numEntries {
		return s.allocNode(items)
	}

	sort.Slice(items, func(i, j int) bool {
		return s.bbox(items[i]).X1 < s.bbox(items[j]).X1
	})

	P := int(math.Ceil(float64(len(items)) / float64(numEntries)))
	S := int(math.Ceil(math.Sqrt(float64(P))))

	slices := strSplit(items, S)

	nodes := strJoin(
		slices, S,
		func(i, j int32) bool { return s.bbox(i).Y1 < s.bbox(j).Y1 },
		s.allocNode,
	)

	return s.buildNodes(nodes)
}

func (s *STRTree) bbox(idx int32) Box2D {
	if idx < 0 { // Don't care if leaf or node.
		idx = -idx - 1
	}
	return s.bboxes[s.index[idx]]
}

func strSplit(items []int32, S int) [][]int32 {
	sm := S * numEntries
	slices := make([][]int32, S)
	for i := range slices {
		var n int
		if l := len(items); l < sm {
			n = l
		} else {
			n = sm
		}
		slices[i] = items[:n]
		items = items[n:]
	}
	return slices
}

func strJoin(
	slices [][]int32, S int,
	less func(int32, int32) bool,
	alloc func([]int32) int32,
) []int32 {
	nodes := make([]int32, 0, S*S)

	for _, slice := range slices {
		sort.Slice(slice, func(i, j int) bool {
			return less(slice[i], slice[j])
		})

		for len(slice) > 0 {
			var n int
			if l := len(slice); l >= numEntries {
				n = numEntries
			} else {
				n = l
			}
			nodes = append(nodes, alloc(slice[:n]))
			slice = slice[n:]
		}
	}
	return nodes
}

func (s *STRTree) allocNode(items []int32) int32 {
	pos := len(s.index)
	s.index = append(s.index, 0, int32(len(items)))
	s.index = append(s.index, items...)
	if len(items) > 0 {
		box := s.bbox(items[0])
		for i := 1; i < len(items); i++ {
			box = box.Union(s.bbox(items[i]))
		}
		s.index[pos] = int32(s.allocBBox(box))
	}
	return int32(pos)
}

func (s *STRTree) allocBBox(box Box2D) int {
	pos := len(s.bboxes)
	s.bboxes = append(s.bboxes, box)
	return pos
}

func (s *STRTree) allocLeaf(items []int32) int32 {
	pos := len(s.index)
	s.index = append(s.index, 0, int32(len(items)))
	s.index = append(s.index, items...)
	if len(items) > 0 {
		vertices := s.tin.Vertices
		ti := s.tin.Triangles[items[0]]
		t := Triangle{
			vertices[ti[0]],
			vertices[ti[1]],
			vertices[ti[2]],
		}
		box := t.BBox()
		for i := 1; i < len(items); i++ {
			it := items[i]
			ti := s.tin.Triangles[it]
			t := Triangle{
				vertices[ti[0]],
				vertices[ti[1]],
				vertices[ti[2]],
			}
			box = box.Union(t.BBox())
		}
		s.index[pos] = int32(s.allocBBox(box))
	}
	return int32(-(pos + 1))
}