Mercurial > gemma
view pkg/octree/polygon.go @ 4488:bff6c5c1db4f
client: pdf-gen: improve adding bottleneck info to pdf
* Check if the bottleneck is in the current view to add its info to the exported pdf and the pdf filename, this avoid wrong filename and wrong info in pdf in case view has been changed to another location.
* Set the bottleneck to print after moving to it in map.
author | Fadi Abbud <fadi.abbud@intevation.de> |
---|---|
date | Fri, 27 Sep 2019 11:15:02 +0200 |
parents | 1c3df921361d |
children | a38d846d9fd5 |
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) 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.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 } 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} }