Mercurial > gemma

comparison pkg/octree/vertex.go @ 788:9f3a4a60dc04

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleaned up interpolation mess in vertical triangle intersection.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 27 Sep 2018 08:15:28 +0200
parents 3d927e06b92c
children d43e61044ad8
comparison
equal deleted inserted replaced
787:3d927e06b92c 788:9f3a4a60dc04
426 426
427 func dist(x1, y1, x2, y2 float64) float64 { 427 func dist(x1, y1, x2, y2 float64) float64 {
428 return math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)) 428 return math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))
429 } 429 }
430 430
431 func interpolate(v1, v2 Vertex) func(x, y float64) float64 { 431 func relative(v1, v2 Vertex) func(x, y float64) float64 {
432 432
433 length := dist(v1.X, v1.Y, v2.X, v2.Y) 433 length := dist(v1.X, v1.Y, v2.X, v2.Y)
434 434
435 // f(0) = v1.Z 435 // f(0) = v1.Z
436 // f(length) = v2.Z 436 // f(length) = v2.Z
441 m := (v2.Z - v1.Z) / length 441 m := (v2.Z - v1.Z) / length
442 442
443 return func(x, y float64) float64 { 443 return func(x, y float64) float64 {
444 l := dist(v1.X, v1.Y, x, y) 444 l := dist(v1.X, v1.Y, x, y)
445 return m*l + b 445 return m*l + b
446 }
447 }
448
449 func interpolate(a, b float64) func(float64) float64 {
450 return func(x float64) float64 {
451 return (b-a)*x + a
446 } 452 }
447 } 453 }
448 454
449 func linear(v1, v2 Vertex) func(t float64) Vertex { 455 func linear(v1, v2 Vertex) func(t float64) Vertex {
450 return func(t float64) Vertex { 456 return func(t float64) Vertex {
454 (v2.Z-v1.Z)*t + v1.Z, 460 (v2.Z-v1.Z)*t + v1.Z,
455 } 461 }
456 } 462 }
457 } 463 }
458 464
459 func inRange(a, b, c float64) bool { 465 func inRange(a float64) bool { return 0 <= a && a <= 1 }
460 if a > b {
461 a, b = b, a
462 }
463 return c >= a && c <= b
464 }
465 466
466 func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ { 467 func (vl *VerticalLine) Intersection(t *Triangle) LineStringZ {
467 468
468 var out LineStringZ 469 var out LineStringZ
469
470 length := dist(vl.x1, vl.y1, vl.x2, vl.y2)
471 470
472 edges: 471 edges:
473 for i := 0; i < 3 && len(out) < 2; i++ { 472 for i := 0; i < 3 && len(out) < 2; i++ {
474 j := (i + 1) % 3 473 j := (i + 1) % 3
475 edge := NewPlane2D(t[i].X, t[i].Y, t[j].X, t[j].Y) 474 edge := NewPlane2D(t[i].X, t[i].Y, t[j].X, t[j].Y)
480 o1 := onPlane(s1) 479 o1 := onPlane(s1)
481 o2 := onPlane(s2) 480 o2 := onPlane(s2)
482 481
483 switch { 482 switch {
484 case o1 && o2: 483 case o1 && o2:
485 // XXX: Broken! 484 pos := relative(t[i], t[j])
486 pos := interpolate(t[i], t[j])
487 t1 := pos(vl.x1, vl.y1) 485 t1 := pos(vl.x1, vl.y1)
488 t2 := pos(vl.x2, vl.y2) 486 t2 := pos(vl.x2, vl.y2)
489 487
490 r1 := inRange(0, length, t1) 488 r1 := inRange(t1)
491 r2 := inRange(0, length, t2) 489 r2 := inRange(t2)
492 490
493 switch { 491 switch {
494 case r1 && r2: 492 case r1 && r2:
495 l := linear(t[i], t[j]) 493 lin := linear(t[i], t[j])
496 // XXX: Broken! 494 out = append(out, lin(t1), lin(t2))
497 out = append(out, l(t1/length), l(t2/length))
498 return out 495 return out
499 496
500 case !r1 && !r2: // the hole edge 497 case !r1 && !r2: // the hole edge
501 out = append(out, t[i], t[j]) 498 out = append(out, t[i], t[j])
502 return out 499 return out
505 case !r2: 502 case !r2:
506 // TODO 503 // TODO
507 } 504 }
508 505
509 case o1: 506 case o1:
510 t1 := interpolate(t[i], t[j])(vl.x1, vl.y1) 507 t1 := relative(t[i], t[j])(vl.x1, vl.y1)
511 if !inRange(0, length, t1) { 508 if !inRange(t1) {
512 continue edges 509 continue edges
513 } 510 }
514 out = append(out, linear(t[i], t[j])(t1)) 511 out = append(out, linear(t[i], t[j])(t1))
515 512
516 case o2: 513 case o2:
517 t2 := interpolate(t[i], t[j])(vl.x2, vl.y2) 514 t2 := relative(t[i], t[j])(vl.x2, vl.y2)
518 if !inRange(0, length, t2) { 515 if !inRange(t2) {
519 continue edges 516 continue edges
520 } 517 }
521 out = append(out, linear(t[i], t[j])(t2)) 518 out = append(out, linear(t[i], t[j])(t2))
522 519
523 default: 520 default:
524 x, y, intersects := vl.line.Intersection(edge) 521 x, y, intersects := vl.line.Intersection(edge)
525 if !intersects { 522 if !intersects {
526 continue edges 523 continue edges
527 } 524 }
528 t1 := interpolate(t[i], t[j])(x, y) 525 t1 := relative(t[i], t[j])(x, y)
529 if !inRange(0, length, t1) { 526 if !inRange(t1) {
530 continue edges 527 continue edges
531 } 528 }
532 529
533 // XXX: Ugly! 530 z := interpolate(t[j].Z, t[i].Z)(t1)
534 z := (t[j].Z-t[i].Z)*(t1/length) + t[1].Z
535 n := Vertex{x, y, z} 531 n := Vertex{x, y, z}
536 532
537 if math.Signbit(s1) != math.Signbit(s2) { 533 if math.Signbit(s1) != math.Signbit(s2) {
538 // different sides -> vertex on edge. 534 // different sides -> vertex on edge.
539 out = append(out) 535 out = append(out, n)
540 } else { // same side -> inside. Find peer 536 } else { // same side -> inside. Find peer
541 if len(out) > 0 { // we already have our peer 537 if len(out) > 0 { // we already have our peer
542 out = append(out, n) 538 out = append(out, n)
543 return out 539 return out
544 } 540 }
551 other := NewPlane2D(t[k].X, t[k].Y, t[l].X, t[l].Y) 547 other := NewPlane2D(t[k].X, t[k].Y, t[l].X, t[l].Y)
552 xo, yo, intersects := vl.line.Intersection(edge) 548 xo, yo, intersects := vl.line.Intersection(edge)
553 if !intersects { 549 if !intersects {
554 continue 550 continue
555 } 551 }
556 // XXX: Broken!!! 552 t2 := relative(t[k], t[j])(xo, yo)
557 zo := interpolate(t[k], t[j])(xo, yo) 553 if !inRange(t2) {
558 if !inRange(t[k].Z, t[k].Z, zo) {
559 continue 554 continue
560 } 555 }
556 zo := interpolate(t[k].Z, t[j].Z)(t2)
557
561 m := Vertex{xo, yo, zo} 558 m := Vertex{xo, yo, zo}
562 559
563 var xn, yn, xf, yf float64 560 var xn, yn, xf, yf float64
564 561
565 // Which is nearer to current edge? 562 // Which is nearer to current edge?
572 } 569 }
573 570
574 if onPlane(other.Eval(xn, yn)) { 571 if onPlane(other.Eval(xn, yn)) {
575 // triangle intersect in only point 572 // triangle intersect in only point
576 // treat as no intersection 573 // treat as no intersection
577 return out 574 break edges
578 } 575 }
579 576
580 // XXX: Broken!!! 577 pos := relative(n, m)
581 zint := interpolate(n, m) 578
582 579 tzn := pos(xn, yn)
583 zn := zint(xn, yn) 580 tzf := pos(xf, yf)
584 zf := zint(xf, yf) 581
585 582 if !inRange(tzn) {
586 if !inRange(n.Z, m.Z, zn) {
587 // if near is out of range far is, too. 583 // if near is out of range far is, too.
588 return out 584 return out
589 } 585 }
590 586
591 if inRange(n.Z, m.Z, zf) { 587 lin := interpolate(n.Z, m.Z)
588
589 if inRange(tzf) {
592 m.X = xf 590 m.X = xf
593 m.Y = yf 591 m.Y = yf
594 m.Z = zf 592 m.Z = lin(tzf)
595 } // else m is clipping 593 } // else m is clipping
596 594
597 n.X = xn 595 n.X = xn
598 n.Y = yn 596 n.Y = yn
599 n.Z = zn 597 n.Z = lin(tzn)
600 598
601 out = append(out, n, m) 599 out = append(out, n, m)
602 return out 600 return out
603 } 601 }
604 } 602 }