Mercurial > gemma
view pkg/octree/polygon.go @ 3488:aa264021e2bf zpg-ldc
merged default into zpg-ldc branch.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 27 May 2019 16:38:07 +0200 |
parents | 857bb070b9f1 |
children | 01ce3ba9b0d0 |
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" "fmt" "log" "math" "time" "github.com/tidwall/rtree" "gemma.intevation.de/gemma/pkg/wkb" ) type ( ring []float64 Polygon struct { rings []ring indices []*rtree.RTree } IntersectionType byte lineSegment []float64 ) const ( clippingPolygonSQL = ` WITH joined AS ( SELECT sr.area AS area, sr.date_info AS date_info FROM waterway.sounding_results sr WHERE sr.bottleneck_id = $1 ) SELECT ST_AsBinary( ST_intersection( (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date), (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date) ) ) AS area ` ) const ( IntersectionInside IntersectionType = iota IntersectionOutSide IntersectionOverlaps ) func LoadClippingPolygon( ctx context.Context, conn *sql.Conn, epsg uint32, bottleneck string, first, second time.Time, ) (*Polygon, error) { var clip []byte if err := conn.QueryRowContext( ctx, clippingPolygonSQL, bottleneck, epsg, first, second, ).Scan(&clip); err != nil { return nil, err } var polygon Polygon if err := polygon.FromWKB(clip); err != nil { return nil, err } polygon.Indexify() return &polygon, nil } func (ls lineSegment) Rect(interface{}) ([]float64, []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 := rtree.New(nil) 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]} } index.Insert(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 } for _, index := range p.indices { var intersects bool index.Search(box, func(item rtree.Item) 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. point := []float64{box.X1, box.Y1} // Check holes first: inside a hole means outside. if len(p.rings) > 1 { for _, hole := range p.rings[1:] { if hole.contains(point) { return IntersectionOutSide } } } // Check shell if p.rings[0].contains(point) { return IntersectionInside } return IntersectionOutSide } func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType { box := t.BBox() for _, index := range p.indices { var intersects bool index.Search(box, func(item rtree.Item) 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. point := []float64{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 hole.contains(point) { return IntersectionOutSide } } } // Check shell if p.rings[0].contains(point) { return IntersectionInside } return IntersectionOutSide } func (rng ring) isClosed() bool { return (len(rng) / 2) >= 3 } func (rng ring) contains(point []float64) bool { if !rng.isClosed() { return false } end := len(rng)/2 - 1 contains := intersectsWithRaycast(point, rng[:2], rng[end*2:end*2+2]) for i := 2; i < len(rng); i += 2 { if intersectsWithRaycast(point, rng[i-2:i], rng[i:i+2]) { contains = !contains } } return contains } // 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(point, start, end []float64) bool { // Always ensure that the the first point // has a y coordinate that is less than the second point if start[1] > end[1] { // Switch the points if otherwise. start, end = end, start } // Move the point's y coordinate // outside of the bounds of the testing region // so we can start drawing a ray for point[1] == start[1] || point[1] == end[1] { y := math.Nextafter(point[1], math.Inf(1)) point = []float64{point[0], y} } // If we are outside of the polygon, indicate so. if point[1] < start[1] || point[1] > end[1] { return false } if start[0] > end[0] { if point[0] > start[0] { return false } if point[0] < end[0] { return true } } else { if point[0] > end[0] { return false } if point[0] < start[0] { return true } } raySlope := (point[1] - start[1]) / (point[0] - start[0]) diagSlope := (end[1] - start[1]) / (end[0] - start[0]) return raySlope >= diagSlope } 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.Printf("info: 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.Printf("info: 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 }