Mercurial > gemma
view pkg/mesh/polygon.go @ 5652:9b3779f649c2 v5.7.1
Bumped version for release.
author | Sascha Wilde <wilde@sha-bang.de> |
---|---|
date | Mon, 07 Aug 2023 17:59:18 +0200 |
parents | 1222b777f51f |
children | 6270951dda28 |
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 has a border and holes. // The line segments are spatially indexed. Polygon struct { rings []ring indices []*rtree.RTree } // IntersectionType represents an enum // of the type of intersection. IntersectionType byte lineSegment []float64 ) const ( // IntersectionInside is inside the polygon. IntersectionInside IntersectionType = iota // IntersectionOutSide is outside the polygon. IntersectionOutSide // IntersectionOverlaps overlaps the polygon. 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 } // Indexify creates a spatial index for thw polygon. 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 } // IntersectionBox2D checks the type of intersection of the // given box. 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 } // IntersectionWithTriangle checks the intersection type // for the given triangle. 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 } // NumVertices returns the number of vertices of a given ring. func (p *Polygon) NumVertices(ring int) int { if ring < 0 || ring >= len(p.rings) { return 0 } return len(p.rings[ring]) / 2 } // Vertices passes the vertices of a given ring // to the given fn function. 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]) } } // AsWKB serializes the polygon as WKB. 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() } // FromWKB deserializes a polygon from WKB. 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 } // FromLineStringZ creates a polygon from a given linestring z. 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} }