Mercurial > gemma
comparison pkg/mesh/polygon.go @ 4827:f4abfd0ee8ad remove-octree-debris
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 } |