Mercurial > gemma
view pkg/octree/tin.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 | f4dcbe8941a1 |
children | 63475c8e710e |
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" "gemma.intevation.de/gemma/pkg/models" ) 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 } // 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.Printf("info: 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("info: 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 } 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 ) // GenerateTinByID generated a TIN by triangulating a point cloud // from the database. 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 } // 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 } for _, v := range t.Vertices { if err := v.Write(w); err != nil { return err } } log.Printf("info: 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("info: compressed tin indices in bytes: %d (%d)\n", written, 3*4*len(t.Triangles)) return nil }