view pkg/octree/tin.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 d753ce6cf588
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 (
	"bytes"
	"context"
	"database/sql"
	"encoding/binary"
	"errors"
	"fmt"
	"io"
	"log"
	"math"
	"time"
)

var (
	errNoByteSlice   = errors.New("Not a byte slice")
	errTooLessPoints = errors.New("Too less points")
)

const wgs84 = 4326

type Tin struct {
	EPSG      uint32
	Vertices  []Vertex
	Triangles [][]int32

	Min Vertex
	Max Vertex
}

func (t *Tin) FromWKB(data []byte) error {
	log.Printf("data length %d\n", len(data))

	r := bytes.NewReader(data)

	endian, err := r.ReadByte()

	var order binary.ByteOrder

	switch {
	case err != nil:
		return err
	case endian == wkbNDR:
		order = binary.LittleEndian
	case endian == wkbXDR:
		order = binary.BigEndian
	default:
		return fmt.Errorf("unknown byte order %x", endian)
	}

	var geomType uint32
	err = binary.Read(r, order, &geomType)

	switch {
	case err != nil:
		return err
	case geomType != wkbTinZ:
		return fmt.Errorf("unknown geometry type %x", geomType)
	}

	var num uint32
	if err = binary.Read(r, order, &num); err != nil {
		return err
	}

	vertices := make([]Vertex, 0, 100000)

	var v Vertex

	v2i := make(map[Vertex]int32, 100000)

	var indexPool []int32

	allocIndices := func() []int32 {
		if len(indexPool) == 0 {
			indexPool = make([]int32, 3*8*1024)
		}
		ids := indexPool[:3]
		indexPool = indexPool[3:]
		return ids
	}

	var triangles [][]int32

	min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64}
	max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64}

	for i := uint32(0); i < num; i++ {

		endian, err = r.ReadByte()
		switch {
		case err != nil:
			return err
		case endian == wkbNDR:
			order = binary.LittleEndian
		case endian == wkbXDR:
			order = binary.BigEndian
		default:
			return fmt.Errorf("unknown byte order %x", endian)
		}

		err = binary.Read(r, order, &geomType)
		switch {
		case err != nil:
			return err
		case geomType != wkbTriangleZ:
			return fmt.Errorf("unknown geometry type %d", geomType)
		}

		var rings uint32
		if err = binary.Read(r, order, &rings); err != nil {
			return err
		}
		triangle := allocIndices()

		for ring := uint32(0); ring < rings; ring++ {
			var npoints uint32
			if err = binary.Read(r, order, &npoints); err != nil {
				return err
			}

			if npoints < 3 {
				return errTooLessPoints
			}

			for p := uint32(0); p < npoints; p++ {
				var x, y, z uint64
				for _, addr := range []*uint64{&x, &y, &z} {
					if err = binary.Read(r, order, addr); err != nil {
						return err
					}
				}
				if p >= 3 || ring >= 1 {
					// Don't store the forth point.
					continue
				}
				// Do this conversion later to spare reflect calls
				// and allocs in binary.Read.
				v.X = math.Float64frombits(x)
				v.Y = math.Float64frombits(y)
				v.Z = math.Float64frombits(z)
				idx, found := v2i[v]
				if !found {
					idx = int32(len(vertices))
					v2i[v] = idx
					vertices = append(vertices, v)
					min.Minimize(v)
					max.Maximize(v)
				}
				triangle[p] = idx
			}
		}
		triangles = append(triangles, triangle)
	}

	log.Printf("bbox: [[%f, %f], [%f, %f]]\n",
		min.X, min.Y, max.X, max.Y)

	*t = Tin{
		EPSG:      wgs84,
		Vertices:  vertices,
		Triangles: triangles,
		Min:       min,
		Max:       max,
	}

	return nil
}

const (
	tinSQLPrefix = `WITH trans AS (
  SELECT
	ST_Buffer(ST_Transform(area::geometry, $1::int), 0.001) AS area,
	ST_Transform(point_cloud::geometry, $1::int) AS point_cloud
  FROM waterway.sounding_results
`
	tinSQLSuffix = `
),
triangles AS (
  SELECT t.geom AS geom, ST_MakePolygon(ST_ExteriorRing(t.geom)) AS poly FROM (
    SELECT (ST_Dump(
      ST_DelaunayTriangles(point_cloud, 0, 2))).geom
    FROM trans) t
)
SELECT ST_AsBinary(ST_Collect(triangles.geom)) FROM triangles, trans
WHERE ST_Covers(trans.area, triangles.poly)`

	loadTinSQL     = tinSQLPrefix + `WHERE bottleneck_id = $2 AND date_info = $3` + tinSQLSuffix
	loadTinByIDSQL = tinSQLPrefix + `WHERE id = $2` + tinSQLSuffix
)

func GenerateTin(
	conn *sql.Conn,
	ctx context.Context,
	bottleneck string,
	date time.Time,
	epsg uint32,
) (*Tin, error) {
	var tin Tin
	err := conn.QueryRowContext(ctx, loadTinSQL, epsg, bottleneck, date).Scan(&tin)
	switch {
	case err == sql.ErrNoRows:
		return nil, nil
	case err != nil:
		return nil, err
	}
	tin.EPSG = epsg
	return &tin, nil
}

func GenerateTinByID(
	conn *sql.Conn,
	ctx context.Context,
	id int64,
	epsg uint32,
) (*Tin, error) {
	var tin Tin
	err := conn.QueryRowContext(ctx, loadTinByIDSQL, epsg, id).Scan(&tin)
	switch {
	case err == sql.ErrNoRows:
		return nil, nil
	case err != nil:
		return nil, err
	}
	tin.EPSG = epsg
	return &tin, nil
}

func (t *Tin) Scan(raw interface{}) error {
	if raw == nil {
		return nil
	}
	data, ok := raw.([]byte)
	if !ok {
		return errNoByteSlice
	}
	return t.FromWKB(data)
}

func (t *Tin) Serialize(w io.Writer) error {

	if err := binary.Write(w, binary.LittleEndian, t.EPSG); err != nil {
		return err
	}

	if err := t.Min.Write(w); err != nil {
		return err
	}
	if err := t.Max.Write(w); err != nil {
		return err
	}

	if err := binary.Write(
		w, binary.LittleEndian, uint32(len(t.Vertices))); err != nil {
		return err
	}

	for _, v := range t.Vertices {
		if err := v.Write(w); err != nil {
			return err
		}
	}
	log.Printf("vertices %d (%d)\n", len(t.Vertices), len(t.Vertices)*3*8)

	if err := binary.Write(
		w, binary.LittleEndian, uint32(len(t.Triangles))); err != nil {
		return err
	}

	var buf [binary.MaxVarintLen32]byte
	var written int
	var last int32
	for _, triangle := range t.Triangles {
		for _, idx := range triangle {
			value := idx - last
			n := binary.PutVarint(buf[:], int64(value))
			for p := buf[:n]; len(p) > 0; p = p[n:] {
				var err error
				if n, err = w.Write(p); err != nil {
					return err
				}
				written += n
			}
			last = idx
		}
	}
	log.Printf("compressed tin indices in bytes: %d (%d)\n",
		written, 3*4*len(t.Triangles))

	return nil
}