view pkg/octree/tree.go @ 2130:f3aabc05f9b2

Fix constraints on waterway profiles staging_done in the UNIQUE constraint had no effect, because the exclusion constraint prevented two rows with equal location and validity anyhow. Adding staging_done to the exclusion constraint makes the UNIQUE constraint checking only a corner case of what the exclusion constraint checks. Thus, remove the UNIQUE constraint. Casting staging_done to int is needed because there is no appropriate operator class for booleans. Casting to smallint or even bit would have been better (i.e. should result in smaller index size), but that would have required creating such a CAST, in addition.
author Tom Gottfried <tom@intevation.de>
date Wed, 06 Feb 2019 15:42:32 +0100
parents fe1aa62195c2
children a1e751c08c56
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"
)

// Tree is an Octree holding triangles.
type Tree struct {
	// EPSG is the projection.
	EPSG uint32

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

	// Min is the lower left corner of the bbox.
	Min Vertex
	// Max is the upper right corner of the bbox.
	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},
}

// Vertical does a vertical cross cut from (x1, y1) to (x2, y2).
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
			if index := ot.index[top.pos:]; len(index) > 7 {
				for i := 0; i < 4; i++ {
					a := index[i]
					b := index[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{}{}
			}
		}
	}
}

// Horizontal does a horizontal cross cut.
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
			}
			if index := ot.index[pos:]; len(index) > 3 {
				stack = append(stack,
					frame{index[0], min, max},
					frame{index[1], min, max},
					frame{index[2], min, max},
					frame{index[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)
				}
			}
		}
	}
}