view pkg/mesh/tin.go @ 5520:05db984d3db1

Improve performance of bottleneck area calculation Avoid buffer calculations by replacing them with simple distance comparisons and calculate the boundary of the result geometry only once per iteration. In some edge cases with very large numbers of iterations, this reduced the runtime of a bottleneck import by a factor of more than twenty.
author Tom Gottfried <tom@intevation.de>
date Thu, 21 Oct 2021 19:50:39 +0200
parents 5f47eeea988d
children 1222b777f51f
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 mesh

import (
	"bytes"
	"encoding/binary"
	"errors"
	"fmt"
	"io"
	"math"

	"gemma.intevation.de/gemma/pkg/log"
	"gemma.intevation.de/gemma/pkg/models"
	"gemma.intevation.de/gemma/pkg/wkb"
)

var (
	errNoByteSlice   = errors.New("not a byte slice")
	errTooLessPoints = errors.New("too less points")
)

// Tin stores a mesh of triangles with common vertices.
type Tin struct {
	// EPSG holds the projection.
	EPSG uint32
	// Vertices are the shared vertices.
	Vertices []Vertex
	// Triangles are the triangles.
	Triangles [][]int32

	// Min is the lower left corner of the bbox.
	Min Vertex
	// Max is the upper right corner of the bbox.
	Max Vertex
}

func (t *Tin) Clip(polygon *Polygon) map[int32]struct{} {
	var tree STRTree
	tree.Build(t)
	return tree.Clip(polygon)
}

// FromWKB constructs the TIN from a WKB representation.
// Shared vertices are identified and referenced by the
// same index.
func (t *Tin) FromWKB(data []byte) error {
	log.Infof("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 == wkb.NDR:
		order = binary.LittleEndian
	case endian == wkb.XDR:
		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 != wkb.TinZ:
		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 == wkb.NDR:
			order = binary.LittleEndian
		case endian == wkb.XDR:
			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 != wkb.TriangleZ:
			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.Infof("bbox: [[%f, %f], [%f, %f]]\n",
		min.X, min.Y, max.X, max.Y)

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

	return nil
}

// Scan implements the sql.Scanner interface.
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
	}

	var err error
	vwrite := func(v float64) {
		if err == nil {
			err = binary.Write(w, binary.LittleEndian, math.Float64bits(v))
		}
	}

	for _, v := range t.Vertices {
		vwrite(v.X)
		vwrite(v.Y)
		vwrite(v.Z)
	}

	if err != nil {
		return err
	}
	log.Infof("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.Infof("compressed tin indices in bytes: %d (%d)\n",
		written, 3*4*len(t.Triangles))

	return nil
}