view pkg/octree/tin.go @ 1328:d753ce6cf588

To make golint happier made context.Context to be the first argument in all calls.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sun, 25 Nov 2018 16:26:41 +0100
parents a244b18cb916
children 84e78d2e2d95
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"
)

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

	loadTinByIDSQL = tinSQLPrefix + `WHERE id = $2` + tinSQLSuffix
)

func GenerateTinByID(
	ctx context.Context,
	conn *sql.Conn,
	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
}