view pkg/octree/tree.go @ 1234:1a5564655f2a

refac: Sidebar reorganized In order to make context switches between administrative tasks which are map related and those which are system related, we now have a category "administration" and "systemadministration". The Riverbedmorphology does nothing than display the map, so it is renamed to that (map). In case the context of "systemadministration" is chosen, the "map" brings you just back to the map.
author Thomas Junk <thomas.junk@intevation.de>
date Tue, 20 Nov 2018 09:54:53 +0100
parents a244b18cb916
children f4dcbe8941a1
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"
)

type Tree struct {
	EPSG uint32

	vertices  []Vertex
	triangles [][]int32
	index     []int32

	Min Vertex
	Max Vertex
}

var scale = [4][4]float64{
	{0.0, 0.0, 0.5, 0.5},
	{0.5, 0.0, 1.0, 0.5},
	{0.0, 0.5, 0.5, 1.0},
	{0.5, 0.5, 1.0, 1.0},
}

func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {

	box := Box2D{
		X1: math.Min(x1, x2),
		Y1: math.Min(y1, y2),
		X2: math.Max(x1, x2),
		Y2: math.Max(y1, y2),
	}

	// out of bounding box
	if box.X2 < ot.Min.X || ot.Max.X < box.X1 ||
		box.Y2 < ot.Min.Y || ot.Max.Y < box.Y1 {
		return
	}

	line := NewPlane2D(x1, y1, x2, y2)

	type frame struct {
		pos int32
		Box2D
	}

	dupes := map[int32]struct{}{}

	all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
	//log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1))
	//log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1))

	stack := []frame{{1, all}}

	for len(stack) > 0 {
		top := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		if top.pos > 0 { // node
			for i := int32(0); i < 4; i++ {
				a := ot.index[top.pos+i]
				b := ot.index[top.pos+i+4]
				if a == 0 && b == 0 {
					continue
				}
				dx := top.X2 - top.X1
				dy := top.Y2 - top.Y1
				nbox := Box2D{
					dx*scale[i][0] + top.X1,
					dy*scale[i][1] + top.Y1,
					dx*scale[i][2] + top.X1,
					dy*scale[i][3] + top.Y1,
				}
				if nbox.Intersects(box) && nbox.IntersectsPlane(line) {
					if a != 0 {
						stack = append(stack, frame{a, nbox})
					}
					if b != 0 {
						stack = append(stack, frame{b, nbox})
					}
				}
			}
		} else { // leaf
			pos := -top.pos - 1
			n := ot.index[pos]
			indices := ot.index[pos+1 : pos+1+n]

			for _, idx := range indices {
				if _, found := dupes[idx]; found {
					continue
				}
				tri := ot.triangles[idx]
				t := Triangle{
					ot.vertices[tri[0]],
					ot.vertices[tri[1]],
					ot.vertices[tri[2]],
				}

				v0 := line.Eval(t[0].X, t[0].Y)
				v1 := line.Eval(t[1].X, t[1].Y)
				v2 := line.Eval(t[2].X, t[2].Y)

				if onPlane(v0) || onPlane(v1) || onPlane(v2) ||
					sides(sides(sides(0, v0), v1), v2) == 3 {
					fn(&t)
				}
				dupes[idx] = struct{}{}
			}
		}
	}
}

func (ot *Tree) Horizontal(h float64, fn func(*Triangle)) {

	if h < ot.Min.Z || ot.Max.Z < h {
		return
	}

	type frame struct {
		pos int32
		min float64
		max float64
	}

	dupes := map[int32]struct{}{}

	stack := []frame{{1, ot.Min.Z, ot.Max.Z}}

	for len(stack) > 0 {
		top := stack[len(stack)-1]
		stack = stack[:len(stack)-1]

		pos := top.pos
		if pos == 0 {
			continue
		}
		min, max := top.min, top.max

		if pos > 0 { // node
			if mid := (max-min)*0.5 + min; h >= mid {
				pos += 4 // nodes with z-bit set
				min = mid
			} else {
				max = mid
			}
			stack = append(stack,
				frame{ot.index[pos+0], min, max},
				frame{ot.index[pos+1], min, max},
				frame{ot.index[pos+2], min, max},
				frame{ot.index[pos+3], min, max})
		} else { // leaf
			pos = -pos - 1
			n := ot.index[pos]
			//log.Printf("%d %d %d\n", pos, n, len(ot.index))
			indices := ot.index[pos+1 : pos+1+n]

			for _, idx := range indices {
				if _, found := dupes[idx]; found {
					continue
				}
				tri := ot.triangles[idx]
				t := Triangle{
					ot.vertices[tri[0]],
					ot.vertices[tri[1]],
					ot.vertices[tri[2]],
				}

				if !(math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z)) > h ||
					math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) < h) {
					dupes[idx] = struct{}{}
					fn(&t)
				}
			}
		}
	}
}