view cmd/octree2contour/loader.go @ 687:be90ab542aa7 octree

octree: contouring: Do the math to calculate the intersection points of the triangles and the planes.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 20 Sep 2018 11:32:03 +0200
parents a8672ba9ebad
children 9db4ae29ded9
line wrap: on
line source

package main

import (
	"bufio"
	"encoding/binary"
	"io"
	"log"
	"math"
	"os"

	"github.com/golang/snappy"
)

type octree struct {
	epsg uint32

	vertices  []vertex
	triangles [][]int32
	index     []int32

	min vertex
	max vertex
}

func (v *vertex) read(r io.Reader) error {
	var buf [8]byte
	b := buf[:]
	if _, err := io.ReadFull(r, b); err != nil {
		return nil
	}
	v.x = math.Float64frombits(binary.LittleEndian.Uint64(b))
	if _, err := io.ReadFull(r, b); err != nil {
		return nil
	}
	v.y = math.Float64frombits(binary.LittleEndian.Uint64(b))
	if _, err := io.ReadFull(r, b); err != nil {
		return nil
	}
	v.z = math.Float64frombits(binary.LittleEndian.Uint64(b))
	return nil
}

func loadOctreeReader(r *bufio.Reader) (*octree, error) {
	tree := new(octree)

	if err := binary.Read(r, binary.LittleEndian, &tree.epsg); err != nil {
		return nil, err
	}

	log.Printf("EPSG: %d\n", tree.epsg)

	if err := tree.min.read(r); err != nil {
		return nil, err
	}

	if err := tree.max.read(r); err != nil {
		return nil, err
	}

	log.Printf("BBOX: [[%f, %f, %f], [%f, %f, %f]]\n",
		tree.min.x, tree.min.y, tree.min.z,
		tree.max.x, tree.max.y, tree.max.z)

	var numVertices uint32
	if err := binary.Read(r, binary.LittleEndian, &numVertices); err != nil {
		return nil, err
	}

	log.Printf("vertices: %d\n", numVertices)

	vertices := make([]vertex, numVertices)
	tree.vertices = vertices

	for i := range vertices {
		if err := vertices[i].read(r); err != nil {
			return nil, err
		}
	}

	var numTriangles uint32
	if err := binary.Read(r, binary.LittleEndian, &numTriangles); err != nil {
		return nil, err
	}

	log.Printf("triangles: %d\n", numTriangles)

	indices := make([]int32, 3*numTriangles)
	triangles := make([][]int32, numTriangles)
	tree.triangles = triangles

	var last int32

	for i := range triangles {
		tri := indices[:3]
		indices = indices[3:]
		triangles[i] = tri
		for j := range tri {
			v, err := binary.ReadVarint(r)
			if err != nil {
				return nil, err
			}
			value := int32(v) + last
			tri[j] = value
			last = value
		}
	}

	var numNodes uint32
	if err := binary.Read(r, binary.LittleEndian, &numNodes); err != nil {
		return nil, err
	}

	log.Printf("num nodes: %d\n", numNodes)

	tree.index = make([]int32, numNodes)
	entries := tree.index[1:]

	last = 0
	for i := range entries {
		v, err := binary.ReadVarint(r)
		if err != nil {
			return nil, err
		}
		value := int32(v) + last
		entries[i] = value
		last = value
	}

	return tree, nil
}

func (ot *octree) horizontal(h float64, fn func(*triangle)) {

	type frame struct {
		pos int32
		min float64
		max float64
	}

	if h < ot.min.z || ot.max.z < h {
		return
	}

	stack := []frame{{1, ot.min.z, ot.max.z}}

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

	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)
				}
			}
		}
	}
}

func loadOctree(fname string) (*octree, error) {

	f, err := os.Open(fname)
	if err != nil {
		return nil, err
	}
	defer f.Close()
	return loadOctreeReader(bufio.NewReader(snappy.NewReader(f)))
}