Mercurial > gemma
view pkg/mesh/polygon.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" "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} }