Mercurial > gemma
view pkg/mesh/polygon.go @ 5591:0011f50cf216 surveysperbottleneckid
Removed no longer used alternative api for surveys/ endpoint.
As bottlenecks in the summary for SR imports are now identified by
their id and no longer by the (not guarantied to be unique!) name,
there is no longer the need to request survey data by the name+date
tuple (which isn't reliable anyway). So the workaround was now
reversed.
author | Sascha Wilde <wilde@sha-bang.de> |
---|---|
date | Wed, 06 Apr 2022 13:30:29 +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" "fmt" "math" "github.com/tidwall/rtree" "gemma.intevation.de/gemma/pkg/log" "gemma.intevation.de/gemma/pkg/wkb" ) type ( ring []float64 Polygon struct { rings []ring indices []*rtree.RTree } IntersectionType byte lineSegment []float64 ) const ( IntersectionInside IntersectionType = iota IntersectionOutSide IntersectionOverlaps ) func (ls lineSegment) Rect() ([2]float64, [2]float64) { var min, max [2]float64 if ls[0] < ls[2] { min[0] = ls[0] max[0] = ls[2] } else { min[0] = ls[2] max[0] = ls[0] } if ls[1] < ls[3] { min[1] = ls[1] max[1] = ls[3] } else { min[1] = ls[3] max[1] = ls[1] } return min, max } func (p *Polygon) Indexify() { indices := make([]*rtree.RTree, len(p.rings)) for i := range indices { index := new(rtree.RTree) indices[i] = index rng := p.rings[i] for i := 0; i < len(rng); i += 2 { var ls lineSegment if i+4 <= len(rng) { ls = lineSegment(rng[i : i+4]) } else { ls = []float64{rng[i], rng[i+1], rng[0], rng[1]} } min, max := ls.Rect() index.Insert(min, max, ls) } } p.indices = indices } func (ls lineSegment) intersects(a Box2D) bool { p1x := ls[0] p1y := ls[1] p2x := ls[2] p2y := ls[3] left := a.X1 right := a.X2 top := a.Y1 bottom := a.Y2 // The direction of the ray dx := p2x - p1x dy := p2y - p1y min, max := 0.0, 1.0 var t0, t1 float64 // Left and right sides. // - If the line is parallel to the y axis. if dx == 0 { if p1x < left || p1x > right { return false } } else { // - Make sure t0 holds the smaller value by checking the direction of the line. if dx > 0 { t0 = (left - p1x) / dx t1 = (right - p1x) / dx } else { t1 = (left - p1x) / dx t0 = (right - p1x) / dx } if t0 > min { min = t0 } if t1 < max { max = t1 } if min > max || max < 0 { return false } } // The top and bottom side. // - If the line is parallel to the x axis. if dy == 0 { if p1y < top || p1y > bottom { return false } } else { // - Make sure t0 holds the smaller value by checking the direction of the line. if dy > 0 { t0 = (top - p1y) / dy t1 = (bottom - p1y) / dy } else { t1 = (top - p1y) / dy t0 = (bottom - p1y) / dy } if t0 > min { min = t0 } if t1 < max { max = t1 } if min > max || max < 0 { return false } } // The point of intersection // ix = p1x + dx*min // iy = p1y + dy*min return true } func (ls lineSegment) intersectsLineSegment(o lineSegment) bool { p0 := ls[:2] p1 := ls[2:4] p2 := o[:2] p3 := o[2:4] s10x := p1[0] - p0[0] s10y := p1[1] - p0[1] s32x := p3[0] - p2[0] s32y := p3[1] - p2[1] den := s10x*s32y - s32x*s10y if den == 0 { return false } denPos := den > 0 s02x := p0[0] - p2[0] s02y := p0[1] - p2[1] sNum := s10x*s02y - s10y*s02x if sNum < 0 == denPos { return false } tNum := s32x*s02y - s32y*s02x if tNum < 0 == denPos { return false } if sNum > den == denPos || tNum > den == denPos { return false } // t := tNum / den // intersection at( p0[0] + (t * s10x), p0[1] + (t * s10y) ) return true } func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType { if len(p.rings) == 0 { return IntersectionOutSide } min, max := box.Rect() for _, index := range p.indices { var intersects bool index.Search(min, max, func(_, _ [2]float64, item interface{}) bool { if item.(lineSegment).intersects(box) { intersects = true return false } return true }) if intersects { return IntersectionOverlaps } } // No intersection -> check inside or outside // if an abritrary point is inside or not. // Check holes first: inside a hole means outside. if len(p.rings) > 1 { for _, hole := range p.rings[1:] { if contains(hole, box.X1, box.Y1) { return IntersectionOutSide } } } // Check shell if contains(p.rings[0], box.X1, box.Y1) { return IntersectionInside } return IntersectionOutSide } func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType { box := t.BBox() min, max := box.Rect() for _, index := range p.indices { var intersects bool index.Search(min, max, func(_, _ [2]float64, item interface{}) bool { ls := item.(lineSegment) other := make(lineSegment, 4) for i := range t { other[0] = t[i].X other[1] = t[i].Y other[2] = t[(i+1)%len(t)].X other[3] = t[(i+1)%len(t)].Y if ls.intersectsLineSegment(other) { intersects = true return false } } return true }) if intersects { return IntersectionOverlaps } } // No intersection -> check inside or outside // if an abritrary point is inside or not. pX, pY := t[0].X, t[0].Y // Check holes first: inside a hole means outside. if len(p.rings) > 1 { for _, hole := range p.rings[1:] { if contains(hole, pX, pY) { return IntersectionOutSide } } } // Check shell if contains(p.rings[0], pX, pY) { return IntersectionInside } return IntersectionOutSide } func (rng ring) length() int { return len(rng) / 2 } func (rng ring) point(i int) (float64, float64) { i *= 2 return rng[i], rng[i+1] } type segments interface { length() int point(int) (float64, float64) } func contains(s segments, pX, pY float64) bool { n := s.length() if n < 3 { return false } sX, sY := s.point(0) eX, eY := s.point(n - 1) const eps = 0.0000001 if math.Abs(sX-eX) > eps || math.Abs(sY-eY) > eps { // It's not closed! return false } var inside bool for i := 1; i < n; i++ { eX, eY := s.point(i) if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) { inside = !inside } sX, sY = eX, eY } return inside } // Using the raycast algorithm, this returns whether or not the passed in point // Intersects with the edge drawn by the passed in start and end points. // Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go func intersectsWithRaycast(pX, pY, sX, sY, eX, eY float64) bool { // Always ensure that the the first point // has a y coordinate that is less than the second point if sY > eY { // Switch the points if otherwise. sX, sY, eX, eY = eX, eY, sX, sY } // Move the point's y coordinate // outside of the bounds of the testing region // so we can start drawing a ray for pY == sY || pY == eY { pY = math.Nextafter(pY, math.Inf(1)) } // If we are outside of the polygon, indicate so. if pY < sY || pY > eY { return false } if sX > eX { if pX > sX { return false } if pX < eX { return true } } else { if pX > eX { return false } if pX < sX { return true } } raySlope := (pY - sY) / (pX - sX) diagSlope := (eY - sY) / (eX - sX) return raySlope >= diagSlope } func (p *Polygon) NumVertices(ring int) int { if ring < 0 || ring >= len(p.rings) { return 0 } return len(p.rings[ring]) / 2 } func (p *Polygon) Vertices(ring int, fn func(float64, float64)) { if ring < 0 || ring >= len(p.rings) { return } rng := p.rings[ring] for i := 0; i < len(rng); i += 2 { fn(rng[i+0], rng[i+1]) } } func (p *Polygon) AsWKB() []byte { size := 1 + 4 + 4 for _, r := range p.rings { size += 4 + len(r)*8 } buf := bytes.NewBuffer(make([]byte, 0, size)) binary.Write(buf, binary.LittleEndian, wkb.NDR) binary.Write(buf, binary.LittleEndian, wkb.Polygon) binary.Write(buf, binary.LittleEndian, uint32(len(p.rings))) for _, r := range p.rings { binary.Write(buf, binary.LittleEndian, uint32(len(r)/2)) for i := 0; i < len(r); i += 2 { binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0])) binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1])) } } return buf.Bytes() } func (p *Polygon) FromWKB(data []byte) error { 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.Polygon: return fmt.Errorf("unknown geometry type %x", geomType) } var numRings uint32 if err = binary.Read(r, order, &numRings); err != nil { return err } rngs := make([]ring, numRings) log.Infof("number of rings: %d\n", len(rngs)) for rng := uint32(0); rng < numRings; rng++ { var numVertices uint32 if err = binary.Read(r, order, &numVertices); err != nil { return err } log.Infof("number of vertices in ring %d: %d\n", rng, numVertices) numVertices *= 2 vertices := make([]float64, numVertices) for v := uint32(0); v < numVertices; v += 2 { var lat, lon uint64 if err = binary.Read(r, order, &lat); err != nil { return err } if err = binary.Read(r, order, &lon); err != nil { return err } vertices[v] = math.Float64frombits(lat) vertices[v+1] = math.Float64frombits(lon) } rngs[rng] = vertices } p.rings = rngs return nil } func (p *Polygon) FromLineStringZ(ls LineStringZ) { r := make([]float64, 2*len(ls)) var pos int for i := range ls { r[pos+0] = ls[i].X r[pos+1] = ls[i].Y pos += 2 } p.rings = []ring{r} }