Mercurial > gemma

comparison pkg/octree/polygon.go @ 4729:1137c5a18242 stack-polygons

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Virtualized point in polygon test with an interface to be usable for contours, too.
author Sascha L. Teichmann <teichmann@intevation.de>
date Thu, 17 Oct 2019 23:33:20 +0200
parents a38d846d9fd5
children 5c80a33edd44
comparison
equal deleted inserted replaced
4728:bfbdcf67ae55 4729:1137c5a18242
279 } 279 }
280 } 280 }
281 281
282 // No intersection -> check inside or outside 282 // No intersection -> check inside or outside
283 // if an abritrary point is inside or not. 283 // if an abritrary point is inside or not.
284 point := []float64{box.X1, box.Y1}
285 284
286 // Check holes first: inside a hole means outside. 285 // Check holes first: inside a hole means outside.
287 if len(p.rings) > 1 { 286 if len(p.rings) > 1 {
288 for _, hole := range p.rings[1:] { 287 for _, hole := range p.rings[1:] {
289 if hole.contains(point) { 288 if contains(hole, box.X1, box.Y1) {
290 return IntersectionOutSide 289 return IntersectionOutSide
291 } 290 }
292 } 291 }
293 } 292 }
294 293
295 // Check shell 294 // Check shell
296 if p.rings[0].contains(point) { 295 if contains(p.rings[0], box.X1, box.Y1) {
297 return IntersectionInside 296 return IntersectionInside
298 } 297 }
299 return IntersectionOutSide 298 return IntersectionOutSide
300 } 299 }
301 300
322 return IntersectionOverlaps 321 return IntersectionOverlaps
323 } 322 }
324 } 323 }
325 // No intersection -> check inside or outside 324 // No intersection -> check inside or outside
326 // if an abritrary point is inside or not. 325 // if an abritrary point is inside or not.
327 point := []float64{t[0].X, t[0].Y} 326 pX, pY := t[0].X, t[0].Y
328 327
329 // Check holes first: inside a hole means outside. 328 // Check holes first: inside a hole means outside.
330 if len(p.rings) > 1 { 329 if len(p.rings) > 1 {
331 for _, hole := range p.rings[1:] { 330 for _, hole := range p.rings[1:] {
332 if hole.contains(point) { 331 if contains(hole, pX, pY) {
333 return IntersectionOutSide 332 return IntersectionOutSide
334 } 333 }
335 } 334 }
336 } 335 }
337 336
338 // Check shell 337 // Check shell
339 if p.rings[0].contains(point) { 338 if contains(p.rings[0], pX, pY) {
340 return IntersectionInside 339 return IntersectionInside
341 } 340 }
342 return IntersectionOutSide 341 return IntersectionOutSide
343 } 342 }
344 343
345 func (rng ring) isClosed() bool { return (len(rng) / 2) >= 3 } 344 func (rng ring) closed() bool {
346 345 return (len(rng) / 2) >= 3
347 func (rng ring) contains(point []float64) bool { 346 }
348 if !rng.isClosed() { 347
348 func (rng ring) length() int {
349 return len(rng) / 2
350 }
351
352 func (rng ring) point(i int) (float64, float64) {
353 i *= 2
354 return rng[i], rng[i+1]
355 }
356
357 type segments interface {
358 closed() bool
359 length() int
360 point(int) (float64, float64)
361 }
362
363 func contains(s segments, pX, pY float64) bool {
364 if !s.closed() {
349 return false 365 return false
350 } 366 }
351 367
352 end := len(rng)/2 - 1 368 n := s.length()
353 369
354 contains := intersectsWithRaycast(point, rng[:2], rng[end*2:end*2+2]) 370 sX, sY := s.point(0)
355 371 eX, eY := s.point(n - 1)
356 for i := 2; i < len(rng); i += 2 { 372
357 if intersectsWithRaycast(point, rng[i-2:i], rng[i:i+2]) { 373 contains := intersectsWithRaycast(pX, pY, sX, sY, eX, eY)
374
375 for i := 1; i < n; i++ {
376 sX, sY := s.point(i - 1)
377 eX, eY := s.point(i)
378 if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) {
358 contains = !contains 379 contains = !contains
359 } 380 }
360 } 381 }
361 382
362 return contains 383 return contains
363 } 384 }
364 385
365 // Using the raycast algorithm, this returns whether or not the passed in point 386 // Using the raycast algorithm, this returns whether or not the passed in point
366 // Intersects with the edge drawn by the passed in start and end points. 387 // Intersects with the edge drawn by the passed in start and end points.
367 // Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go 388 // Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go
368 func intersectsWithRaycast(point, start, end []float64) bool { 389 func intersectsWithRaycast(pX, pY, sX, sY, eX, eY float64) bool {
369 390
370 // Always ensure that the the first point 391 // Always ensure that the the first point
371 // has a y coordinate that is less than the second point 392 // has a y coordinate that is less than the second point
372 if start[1] > end[1] { 393 if sY > eY {
373 // Switch the points if otherwise. 394 // Switch the points if otherwise.
374 start, end = end, start 395 sX, sY, eX, eY = eX, eY, sX, sY
375 } 396 }
376 397
377 // Move the point's y coordinate 398 // Move the point's y coordinate
378 // outside of the bounds of the testing region 399 // outside of the bounds of the testing region
379 // so we can start drawing a ray 400 // so we can start drawing a ray
380 for point[1] == start[1] || point[1] == end[1] { 401 for pY == sY || pY == eY {
381 y := math.Nextafter(point[1], math.Inf(1)) 402 pY = math.Nextafter(pY, math.Inf(1))
382 point = []float64{point[0], y}
383 } 403 }
384 404
385 // If we are outside of the polygon, indicate so. 405 // If we are outside of the polygon, indicate so.
386 if point[1] < start[1] || point[1] > end[1] { 406 if pY < sY || pY > eY {
387 return false 407 return false
388 } 408 }
389 409
390 if start[0] > end[0] { 410 if sX > eX {
391 if point[0] > start[0] { 411 if pX > sX {
392 return false 412 return false
393 } 413 }
394 if point[0] < end[0] { 414 if pX < eX {
395 return true 415 return true
396 } 416 }
397 } else { 417 } else {
398 if point[0] > end[0] { 418 if pX > eX {
399 return false 419 return false
400 } 420 }
401 if point[0] < start[0] { 421 if pX < sX {
402 return true 422 return true
403 } 423 }
404 } 424 }
405 425
406 raySlope := (point[1] - start[1]) / (point[0] - start[0]) 426 raySlope := (pY - sY) / (pX - sX)
407 diagSlope := (end[1] - start[1]) / (end[0] - start[0]) 427 diagSlope := (eY - sY) / (eX - sX)
408 428
409 return raySlope >= diagSlope 429 return raySlope >= diagSlope
410 } 430 }
411 431
412 func (p *Polygon) NumVertices(ring int) int { 432 func (p *Polygon) NumVertices(ring int) int {