Mercurial > gemma

comparison pkg/mesh/polygon.go @ 4827:f4abfd0ee8ad remove-octree-debris

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Renamed octree package to mesh as there is no octree any more.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 05 Nov 2019 14:30:22 +0100
parents pkg/octree/polygon.go@ea570f43d7a9
children 866eae1bd888
comparison
equal deleted inserted replaced
4826:ec5afada70ec 4827:f4abfd0ee8ad
1 // This is Free Software under GNU Affero General Public License v >= 3.0
2 // without warranty, see README.md and license for details.
3 //
4 // SPDX-License-Identifier: AGPL-3.0-or-later
5 // License-Filename: LICENSES/AGPL-3.0.txt
6 //
7 // Copyright (C) 2018 by via donau
8 // – Österreichische Wasserstraßen-Gesellschaft mbH
9 // Software engineering by Intevation GmbH
10 //
11 // Author(s):
12 // * Sascha L. Teichmann <sascha.teichmann@intevation.de>
13
14 package mesh
15
16 import (
17 "bytes"
18 "encoding/binary"
19 "fmt"
20 "log"
21 "math"
22
23 "github.com/tidwall/rtree"
24
25 "gemma.intevation.de/gemma/pkg/wkb"
26 )
27
28 type (
29 ring []float64
30
31 Polygon struct {
32 rings []ring
33 indices []*rtree.RTree
34 }
35
36 IntersectionType byte
37
38 lineSegment []float64
39 )
40
41 const (
42 IntersectionInside IntersectionType = iota
43 IntersectionOutSide
44 IntersectionOverlaps
45 )
46
47 func (ls lineSegment) Rect(interface{}) ([]float64, []float64) {
48
49 var min, max [2]float64
50
51 if ls[0] < ls[2] {
52 min[0] = ls[0]
53 max[0] = ls[2]
54 } else {
55 min[0] = ls[2]
56 max[0] = ls[0]
57 }
58
59 if ls[1] < ls[3] {
60 min[1] = ls[1]
61 max[1] = ls[3]
62 } else {
63 min[1] = ls[3]
64 max[1] = ls[1]
65 }
66
67 return min[:], max[:]
68 }
69
70 func (p *Polygon) Indexify() {
71 indices := make([]*rtree.RTree, len(p.rings))
72
73 for i := range indices {
74 index := rtree.New(nil)
75 indices[i] = index
76
77 rng := p.rings[i]
78
79 for i := 0; i < len(rng); i += 2 {
80 var ls lineSegment
81 if i+4 <= len(rng) {
82 ls = lineSegment(rng[i : i+4])
83 } else {
84 ls = []float64{rng[i], rng[i+1], rng[0], rng[1]}
85 }
86 index.Insert(ls)
87 }
88 }
89
90 p.indices = indices
91 }
92
93 func (ls lineSegment) intersects(a Box2D) bool {
94 p1x := ls[0]
95 p1y := ls[1]
96 p2x := ls[2]
97 p2y := ls[3]
98
99 left := a.X1
100 right := a.X2
101 top := a.Y1
102 bottom := a.Y2
103
104 // The direction of the ray
105 dx := p2x - p1x
106 dy := p2y - p1y
107
108 min, max := 0.0, 1.0
109
110 var t0, t1 float64
111
112 // Left and right sides.
113 // - If the line is parallel to the y axis.
114 if dx == 0 {
115 if p1x < left || p1x > right {
116 return false
117 }
118 } else {
119 // - Make sure t0 holds the smaller value by checking the direction of the line.
120 if dx > 0 {
121 t0 = (left - p1x) / dx
122 t1 = (right - p1x) / dx
123 } else {
124 t1 = (left - p1x) / dx
125 t0 = (right - p1x) / dx
126 }
127
128 if t0 > min {
129 min = t0
130 }
131 if t1 < max {
132 max = t1
133 }
134 if min > max || max < 0 {
135 return false
136 }
137 }
138
139 // The top and bottom side.
140 // - If the line is parallel to the x axis.
141 if dy == 0 {
142 if p1y < top || p1y > bottom {
143 return false
144 }
145 } else {
146 // - Make sure t0 holds the smaller value by checking the direction of the line.
147 if dy > 0 {
148 t0 = (top - p1y) / dy
149 t1 = (bottom - p1y) / dy
150 } else {
151 t1 = (top - p1y) / dy
152 t0 = (bottom - p1y) / dy
153 }
154
155 if t0 > min {
156 min = t0
157 }
158 if t1 < max {
159 max = t1
160 }
161 if min > max || max < 0 {
162 return false
163 }
164 }
165
166 // The point of intersection
167 // ix = p1x + dx*min
168 // iy = p1y + dy*min
169 return true
170 }
171
172 func (ls lineSegment) intersectsLineSegment(o lineSegment) bool {
173
174 p0 := ls[:2]
175 p1 := ls[2:4]
176 p2 := o[:2]
177 p3 := o[2:4]
178
179 s10x := p1[0] - p0[0]
180 s10y := p1[1] - p0[1]
181 s32x := p3[0] - p2[0]
182 s32y := p3[1] - p2[1]
183
184 den := s10x*s32y - s32x*s10y
185
186 if den == 0 {
187 return false
188 }
189
190 denPos := den > 0
191
192 s02x := p0[0] - p2[0]
193 s02y := p0[1] - p2[1]
194
195 sNum := s10x*s02y - s10y*s02x
196 if sNum < 0 == denPos {
197 return false
198 }
199
200 tNum := s32x*s02y - s32y*s02x
201 if tNum < 0 == denPos {
202 return false
203 }
204
205 if sNum > den == denPos || tNum > den == denPos {
206 return false
207 }
208
209 // t := tNum / den
210 // intersection at( p0[0] + (t * s10x), p0[1] + (t * s10y) )
211 return true
212 }
213
214 func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType {
215
216 if len(p.rings) == 0 {
217 return IntersectionOutSide
218 }
219
220 for _, index := range p.indices {
221 var intersects bool
222 index.Search(box, func(item rtree.Item) bool {
223 if item.(lineSegment).intersects(box) {
224 intersects = true
225 return false
226 }
227 return true
228 })
229 if intersects {
230 return IntersectionOverlaps
231 }
232 }
233
234 // No intersection -> check inside or outside
235 // if an abritrary point is inside or not.
236
237 // Check holes first: inside a hole means outside.
238 if len(p.rings) > 1 {
239 for _, hole := range p.rings[1:] {
240 if contains(hole, box.X1, box.Y1) {
241 return IntersectionOutSide
242 }
243 }
244 }
245
246 // Check shell
247 if contains(p.rings[0], box.X1, box.Y1) {
248 return IntersectionInside
249 }
250 return IntersectionOutSide
251 }
252
253 func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType {
254 box := t.BBox()
255 for _, index := range p.indices {
256 var intersects bool
257 index.Search(box, func(item rtree.Item) bool {
258 ls := item.(lineSegment)
259 other := make(lineSegment, 4)
260 for i := range t {
261 other[0] = t[i].X
262 other[1] = t[i].Y
263 other[2] = t[(i+1)%len(t)].X
264 other[3] = t[(i+1)%len(t)].Y
265 if ls.intersectsLineSegment(other) {
266 intersects = true
267 return false
268 }
269 }
270 return true
271 })
272 if intersects {
273 return IntersectionOverlaps
274 }
275 }
276 // No intersection -> check inside or outside
277 // if an abritrary point is inside or not.
278 pX, pY := t[0].X, t[0].Y
279
280 // Check holes first: inside a hole means outside.
281 if len(p.rings) > 1 {
282 for _, hole := range p.rings[1:] {
283 if contains(hole, pX, pY) {
284 return IntersectionOutSide
285 }
286 }
287 }
288
289 // Check shell
290 if contains(p.rings[0], pX, pY) {
291 return IntersectionInside
292 }
293 return IntersectionOutSide
294 }
295
296 func (rng ring) length() int {
297 return len(rng) / 2
298 }
299
300 func (rng ring) point(i int) (float64, float64) {
301 i *= 2
302 return rng[i], rng[i+1]
303 }
304
305 type segments interface {
306 length() int
307 point(int) (float64, float64)
308 }
309
310 func contains(s segments, pX, pY float64) bool {
311
312 n := s.length()
313 if n < 3 {
314 return false
315 }
316
317 sX, sY := s.point(0)
318 eX, eY := s.point(n - 1)
319
320 const eps = 0.0000001
321
322 if math.Abs(sX-eX) > eps || math.Abs(sY-eY) > eps {
323 // It's not closed!
324 return false
325 }
326
327 var inside bool
328
329 for i := 1; i < n; i++ {
330 eX, eY := s.point(i)
331 if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) {
332 inside = !inside
333 }
334 sX, sY = eX, eY
335 }
336
337 return inside
338 }
339
340 // Using the raycast algorithm, this returns whether or not the passed in point
341 // Intersects with the edge drawn by the passed in start and end points.
342 // Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go
343 func intersectsWithRaycast(pX, pY, sX, sY, eX, eY float64) bool {
344
345 // Always ensure that the the first point
346 // has a y coordinate that is less than the second point
347 if sY > eY {
348 // Switch the points if otherwise.
349 sX, sY, eX, eY = eX, eY, sX, sY
350 }
351
352 // Move the point's y coordinate
353 // outside of the bounds of the testing region
354 // so we can start drawing a ray
355 for pY == sY || pY == eY {
356 pY = math.Nextafter(pY, math.Inf(1))
357 }
358
359 // If we are outside of the polygon, indicate so.
360 if pY < sY || pY > eY {
361 return false
362 }
363
364 if sX > eX {
365 if pX > sX {
366 return false
367 }
368 if pX < eX {
369 return true
370 }
371 } else {
372 if pX > eX {
373 return false
374 }
375 if pX < sX {
376 return true
377 }
378 }
379
380 raySlope := (pY - sY) / (pX - sX)
381 diagSlope := (eY - sY) / (eX - sX)
382
383 return raySlope >= diagSlope
384 }
385
386 func (p *Polygon) NumVertices(ring int) int {
387 if ring < 0 || ring >= len(p.rings) {
388 return 0
389 }
390 return len(p.rings[ring]) / 2
391 }
392
393 func (p *Polygon) Vertices(ring int, fn func(float64, float64)) {
394 if ring < 0 || ring >= len(p.rings) {
395 return
396 }
397 rng := p.rings[ring]
398 for i := 0; i < len(rng); i += 2 {
399 fn(rng[i+0], rng[i+1])
400 }
401 }
402
403 func (p *Polygon) AsWKB() []byte {
404
405 size := 1 + 4 + 4
406 for _, r := range p.rings {
407 size += 4 + len(r)*8
408 }
409
410 buf := bytes.NewBuffer(make([]byte, 0, size))
411
412 binary.Write(buf, binary.LittleEndian, wkb.NDR)
413 binary.Write(buf, binary.LittleEndian, wkb.Polygon)
414 binary.Write(buf, binary.LittleEndian, uint32(len(p.rings)))
415
416 for _, r := range p.rings {
417 binary.Write(buf, binary.LittleEndian, uint32(len(r)/2))
418 for i := 0; i < len(r); i += 2 {
419 binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0]))
420 binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1]))
421 }
422 }
423
424 return buf.Bytes()
425 }
426
427 func (p *Polygon) FromWKB(data []byte) error {
428
429 r := bytes.NewReader(data)
430
431 endian, err := r.ReadByte()
432
433 var order binary.ByteOrder
434
435 switch {
436 case err != nil:
437 return err
438 case endian == wkb.NDR:
439 order = binary.LittleEndian
440 case endian == wkb.XDR:
441 order = binary.BigEndian
442 default:
443 return fmt.Errorf("unknown byte order %x", endian)
444 }
445
446 var geomType uint32
447 err = binary.Read(r, order, &geomType)
448
449 switch {
450 case err != nil:
451 return err
452 case geomType != wkb.Polygon:
453 return fmt.Errorf("unknown geometry type %x", geomType)
454 }
455
456 var numRings uint32
457 if err = binary.Read(r, order, &numRings); err != nil {
458 return err
459 }
460
461 rngs := make([]ring, numRings)
462
463 log.Printf("info: Number of rings: %d\n", len(rngs))
464
465 for rng := uint32(0); rng < numRings; rng++ {
466 var numVertices uint32
467 if err = binary.Read(r, order, &numVertices); err != nil {
468 return err
469 }
470
471 log.Printf("info: Number of vertices in ring %d: %d\n", rng, numVertices)
472
473 numVertices *= 2
474 vertices := make([]float64, numVertices)
475
476 for v := uint32(0); v < numVertices; v += 2 {
477 var lat, lon uint64
478 if err = binary.Read(r, order, &lat); err != nil {
479 return err
480 }
481 if err = binary.Read(r, order, &lon); err != nil {
482 return err
483 }
484 vertices[v] = math.Float64frombits(lat)
485 vertices[v+1] = math.Float64frombits(lon)
486 }
487
488 rngs[rng] = vertices
489 }
490
491 p.rings = rngs
492
493 return nil
494 }
495
496 func (p *Polygon) FromLineStringZ(ls LineStringZ) {
497 r := make([]float64, 2*len(ls))
498 var pos int
499 for i := range ls {
500 r[pos+0] = ls[i].X
501 r[pos+1] = ls[i].Y
502 pos += 2
503 }
504 p.rings = []ring{r}
505 }