view pkg/octree/strtree.go @ 2512:2768f74d54ab octree-diff

Added an STRTree for the triangulation. Should be better suited for filtering out the clipped triangles.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 05 Mar 2019 13:05:50 +0100
parents
children 1534df518201
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 {
	tri       *Triangulation
	index     []int32
	triangles [][]int32
	bboxes    []Box2D
}

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

	s.tri = t

	s.triangles = t.TriangleSlices()

	all := make([]int32, len(s.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) build(items []int32) int32 {
	sort.Slice(items, func(i, j int) bool {
		ti := s.triangles[items[i]]
		tj := s.triangles[items[j]]
		return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X
	})

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

	sm := S * numEntries

	slices := make([][]int32, S)
	for i := range slices {
		end := (i + 1) * sm
		if end > len(items) {
			end = len(items)
		}
		slices[i] = items[i*sm : end]
	}

	leaves := make([]int32, 0, S*S)

	for _, slice := range slices {
		sort.Slice(slice, func(i, j int) bool {
			ti := s.triangles[slice[i]]
			tj := s.triangles[slice[j]]
			return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X
		})

		for len(slice) > 0 {
			n := numEntries
			if l := len(slice); l < numEntries {
				n = l
			}
			leaves = append(leaves, s.allocLeaf(slice[:n]))
			slice = slice[n:]
		}
	}

	return s.buildNodes(leaves)
}

func bboxIndex(idx int32) int32 {
	if idx < 0 {
		return -idx - 1
	}
	return idx
}

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

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

	sort.Slice(items, func(i, j int) bool {
		bi := s.bboxes[s.index[bboxIndex(items[i])]]
		bj := s.bboxes[s.index[bboxIndex(items[j])]]
		return bi.X1 < bj.X1
	})

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

	sm := S * numEntries

	slices := make([][]int32, S)
	for i := range slices {
		end := (i + 1) * sm
		if end > len(items) {
			end = len(items)
		}
		slices[i] = items[i*sm : end]
	}

	nodes := make([]int32, 0, S*S)

	for _, slice := range slices {
		sort.Slice(slice, func(i, j int) bool {
			bi := s.bboxes[s.index[bboxIndex(slice[i])]]
			bj := s.bboxes[s.index[bboxIndex(slice[j])]]
			return bi.X1 < bj.X1
		})

		for len(slice) > 0 {
			n := numEntries
			if l := len(slice); l < numEntries {
				n = l
			}
			nodes = append(nodes, s.allocNode(slice[:n]))
			slice = slice[n:]
		}
	}

	return s.buildNodes(nodes)
}

func (s *STRTree) allocNode(items []int32) int32 {
	pos := len(s.index)
	s.index = append(s.index, 0)
	s.index = append(s.index, items...)
	if len(items) > 0 {
		box := s.bboxes[s.index[bboxIndex(items[0])]]
		for i := 1; i < len(items); i++ {
			it := items[i]
			b := s.bboxes[s.index[bboxIndex(it)]]
			box = box.Union(b)
		}
		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)
	s.index = append(s.index, items...)
	if len(items) > 0 {
		ti := s.triangles[items[0]]
		t := Triangle{
			s.tri.Points[ti[0]],
			s.tri.Points[ti[1]],
			s.tri.Points[ti[2]],
		}
		box := t.BBox()
		for i := 1; i < len(items); i++ {
			it := items[i]
			ti := s.triangles[it]
			t := Triangle{
				s.tri.Points[ti[0]],
				s.tri.Points[ti[1]],
				s.tri.Points[ti[2]],
			}
			box = box.Union(t.BBox())
		}
		s.index[pos] = int32(s.allocBBox(box))
	}
	return -int32(pos) - 1
}