1# geometry.tcl -- 2# 3# Collection of geometry functions. 4# 5# Copyright (c) 2001 by Ideogramic ApS and other parties. 6# Copyright (c) 2004 Arjen Markus 7# Copyright (c) 2010 Andreas Kupries 8# Copyright (c) 2010 Kevin Kenny 9# 10# See the file "license.terms" for information on usage and redistribution 11# of this file, and for a DISCLAIMER OF ALL WARRANTIES. 12# 13# RCS: @(#) $Id: geometry.tcl,v 1.12 2010/05/24 21:44:16 andreas_kupries Exp $ 14 15namespace eval ::math::geometry {} 16 17package require math 18 19### 20# 21# POINTS 22# 23# A point P consists of an x-coordinate, Px, and a y-coordinate, Py, 24# and both coordinates are floating point values. 25# 26# Points are usually denoted by A, B, C, P, or Q. 27# 28### 29# 30# LINES 31# 32# There are basically three types of lines: 33# line A line is defined by two points A and B as the 34# _infinite_ line going through these two points. 35# Often a line is given as a list of 4 coordinates 36# instead of 2 points. 37# line segment A line segment is defined by two points A and B 38# as the _finite_ that starts in A and ends in B. 39# Often a line segment is given as a list of 4 40# coordinates instead of 2 points. 41# polyline A polyline is a sequence of connected line segments. 42# 43# Please note that given a point P, the closest point on a line is given 44# by the projection of P onto the line. The closest point on a line segment 45# may be the projection, but it may also be one of the end points of the 46# line segment. 47# 48### 49# 50# DISTANCES 51# 52# The distances in this package are all floating point values. 53# 54### 55 56# Point constructor 57proc ::math::geometry::p {x y} { 58 return [list $x $y] 59} 60 61# Vector addition 62proc ::math::geometry::+ {pa pb} { 63 foreach {ax ay} $pa break 64 foreach {bx by} $pb break 65 return [list [expr {$ax + $bx}] [expr {$ay + $by}]] 66} 67 68# Vector difference 69proc ::math::geometry::- {pa pb} { 70 foreach {ax ay} $pa break 71 foreach {bx by} $pb break 72 return [list [expr {$ax - $bx}] [expr {$ay - $by}]] 73} 74 75# Distance between 2 points 76proc ::math::geometry::distance {pa pb} { 77 foreach {ax ay} $pa break 78 foreach {bx by} $pb break 79 return [expr {hypot($bx-$ax,$by-$ay)}] 80} 81 82# Length of a vector 83proc ::math::geometry::length {v} { 84 foreach {x y} $v break 85 return [expr {hypot($x,$y)}] 86} 87 88# Scaling a vector by a factor 89proc ::math::geometry::s* {factor p} { 90 foreach {x y} $p break 91 return [list [expr {$x * $factor}] [expr {$y * $factor}]] 92} 93 94# Unit vector into specific direction given by angle (degrees) 95proc ::math::geometry::direction {angle} { 96 variable torad 97 set x [expr { cos($angle * $torad)}] 98 set y [expr {- sin($angle * $torad)}] 99 return [list $x $y] 100} 101 102# Vertical vector of specified length. 103proc ::math::geometry::v {h} { 104 return [list 0 $h] 105} 106 107# Horizontal vector of specified length. 108proc ::math::geometry::h {w} { 109 return [list $w 0] 110} 111 112# Find point on a line between 2 points at a distance 113# distance 0 => a, distance 1 => b 114proc ::math::geometry::between {pa pb s} { 115 return [+ $pa [s* $s [- $pb $pa]]] 116} 117 118# Find direction octant the point (vector) lies in. 119proc ::math::geometry::octant {p} { 120 variable todeg 121 foreach {x y} $p break 122 123 set a [expr {(atan2(-$y,$x)*$todeg)}] 124 while {$a > 360} {set a [expr {$a - 360}]} 125 while {$a < -360} {set a [expr {$a + 360}]} 126 if {$a < 0} {set a [expr {360 + $a}]} 127 128 #puts "p ($x, $y) @ angle $a | [expr {atan2($y,$x)}] | [expr {atan2($y,$x)*$todeg}]" 129 # XXX : Add outer conditions to make a log2 tree of checks. 130 131 if {$a <= 157.5} { 132 if {$a <= 67.5} { 133 if {$a <= 22.5} { return east } 134 return northeast 135 } 136 if {$a <= 112.5} { return north } 137 return northwest 138 } else { 139 if {$a <= 247.5} { 140 if {$a <= 202.5} { return west } 141 return southwest 142 } 143 if {$a <= 337.5} { 144 if {$a <= 292.5} { return south } 145 return southeast 146 } 147 return east ; # a <= 360.0 148 } 149} 150 151# Return the NW and SE corners of the rectangle. 152proc ::math::geometry::nwse {rect} { 153 foreach {xnw ynw xse yse} $rect break 154 return [list [p $xnw $ynw] [p $xse $yse]] 155} 156 157# Construct rectangle from NW and SE corners. 158proc ::math::geometry::rect {pa pb} { 159 foreach {ax ay} $pa break 160 foreach {bx by} $pb break 161 return [list $ax $ay $bx $by] 162} 163 164proc ::math::geometry::conjx {p} { 165 foreach {x y} $p break 166 return [list [expr {- $x}] $y] 167} 168 169proc ::math::geometry::conjy {p} { 170 foreach {x y} $p break 171 return [list $x [expr {- $y}]] 172} 173 174proc ::math::geometry::x {p} { 175 foreach {x y} $p break 176 return $x 177} 178 179proc ::math::geometry::y {p} { 180 foreach {x y} $p break 181 return $y 182} 183 184# ::math::geometry::calculateDistanceToLine 185# 186# Calculate the distance between a point and a line. 187# 188# Arguments: 189# P a point 190# line a line 191# 192# Results: 193# dist the smallest distance between P and the line 194# 195# Examples: 196# - calculateDistanceToLine {5 10} {0 0 10 10} 197# Result: 3.53553390593 198# - calculateDistanceToLine {-10 0} {0 0 10 10} 199# Result: 7.07106781187 200# 201proc ::math::geometry::calculateDistanceToLine {P line} { 202 # solution based on FAQ 1.02 on comp.graphics.algorithms 203 # L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 ) 204 # (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) 205 # s = ----------------------------- 206 # L^2 207 # dist = |s|*L 208 # 209 # => 210 # 211 # | (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) | 212 # dist = --------------------------------- 213 # L 214 set Ax [lindex $line 0] 215 set Ay [lindex $line 1] 216 set Bx [lindex $line 2] 217 set By [lindex $line 3] 218 set Cx [lindex $P 0] 219 set Cy [lindex $P 1] 220 if {$Ax==$Bx && $Ay==$By} { 221 return [lengthOfPolyline [concat $P [lrange $line 0 1]]] 222 } else { 223 set L [expr {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}] 224 return [expr {abs(($Ay-$Cy)*($Bx-$Ax)-($Ax-$Cx)*($By-$Ay)) / $L}] 225 } 226} 227 228# ::math::geometry::findClosestPointOnLine 229# 230# Return the point on a line which is closest to a given point. 231# 232# Arguments: 233# P a point 234# line a line 235# 236# Results: 237# Q the point on the line that has the smallest 238# distance to P 239# 240# Examples: 241# - findClosestPointOnLine {5 10} {0 0 10 10} 242# Result: 7.5 7.5 243# - findClosestPointOnLine {-10 0} {0 0 10 10} 244# Result: -5.0 -5.0 245# 246proc ::math::geometry::findClosestPointOnLine {P line} { 247 return [lindex [findClosestPointOnLineImpl $P $line] 0] 248} 249 250# ::math::geometry::findClosestPointOnLineImpl 251# 252# PRIVATE FUNCTION USED BY OTHER FUNCTIONS. 253# Find the point on a line that is closest to a given point. 254# 255# Arguments: 256# P a point 257# line a line defined by points A and B 258# 259# Results: 260# Q the point on the line that has the smallest 261# distance to P 262# r r has the following meaning: 263# r=0 P = A 264# r=1 P = B 265# r<0 P is on the backward extension of AB 266# r>1 P is on the forward extension of AB 267# 0<r<1 P is interior to AB 268# 269proc ::math::geometry::findClosestPointOnLineImpl {P line} { 270 # solution based on FAQ 1.02 on comp.graphics.algorithms 271 # L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 ) 272 # (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay) 273 # r = ------------------------------- 274 # L^2 275 # Px = Ax + r(Bx-Ax) 276 # Py = Ay + r(By-Ay) 277 set Ax [lindex $line 0] 278 set Ay [lindex $line 1] 279 set Bx [lindex $line 2] 280 set By [lindex $line 3] 281 set Cx [lindex $P 0] 282 set Cy [lindex $P 1] 283 if {$Ax==$Bx && $Ay==$By} { 284 return [list [list $Ax $Ay] 0] 285 } else { 286 set L [expr {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}] 287 set r [expr {(($Cx-$Ax)*($Bx-$Ax) + ($Cy-$Ay)*($By-$Ay))/pow($L,2)}] 288 set Px [expr {$Ax + $r*($Bx-$Ax)}] 289 set Py [expr {$Ay + $r*($By-$Ay)}] 290 return [list [list $Px $Py] $r] 291 } 292} 293 294# ::math::geometry::calculateDistanceToLineSegment 295# 296# Calculate the distance between a point and a line segment. 297# 298# Arguments: 299# P a point 300# linesegment a line segment 301# 302# Results: 303# dist the smallest distance between P and any point 304# on the line segment 305# 306# Examples: 307# - calculateDistanceToLineSegment {5 10} {0 0 10 10} 308# Result: 3.53553390593 309# - calculateDistanceToLineSegment {-10 0} {0 0 10 10} 310# Result: 10.0 311# 312proc ::math::geometry::calculateDistanceToLineSegment {P linesegment} { 313 set result [calculateDistanceToLineSegmentImpl $P $linesegment] 314 set distToLine [lindex $result 0] 315 set r [lindex $result 1] 316 if {$r<0} { 317 return [lengthOfPolyline [concat $P [lrange $linesegment 0 1]]] 318 } elseif {$r>1} { 319 return [lengthOfPolyline [concat $P [lrange $linesegment 2 3]]] 320 } else { 321 return $distToLine 322 } 323} 324 325# ::math::geometry::calculateDistanceToLineSegmentImpl 326# 327# PRIVATE FUNCTION USED BY OTHER FUNCTIONS. 328# Find the distance between a point and a line. 329# 330# Arguments: 331# P a point 332# linesegment a line segment A->B 333# 334# Results: 335# dist the smallest distance between P and the line 336# r r has the following meaning: 337# r=0 P = A 338# r=1 P = B 339# r<0 P is on the backward extension of AB 340# r>1 P is on the forward extension of AB 341# 0<r<1 P is interior to AB 342# 343proc ::math::geometry::calculateDistanceToLineSegmentImpl {P linesegment} { 344 # solution based on FAQ 1.02 on comp.graphics.algorithms 345 # L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 ) 346 # (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) 347 # s = ----------------------------- 348 # L^2 349 # (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay) 350 # r = ------------------------------- 351 # L^2 352 # dist = |s|*L 353 # 354 # => 355 # 356 # | (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) | 357 # dist = --------------------------------- 358 # L 359 set Ax [lindex $linesegment 0] 360 set Ay [lindex $linesegment 1] 361 set Bx [lindex $linesegment 2] 362 set By [lindex $linesegment 3] 363 set Cx [lindex $P 0] 364 set Cy [lindex $P 1] 365 if {$Ax==$Bx && $Ay==$By} { 366 return [list [lengthOfPolyline [concat $P [lrange $linesegment 0 1]]] 0] 367 } else { 368 set L [expr {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}] 369 set r [expr {(($Cx-$Ax)*($Bx-$Ax) + ($Cy-$Ay)*($By-$Ay))/pow($L,2)}] 370 return [list [expr {abs(($Ay-$Cy)*($Bx-$Ax)-($Ax-$Cx)*($By-$Ay)) / $L}] $r] 371 } 372} 373 374# ::math::geometry::findClosestPointOnLineSegment 375# 376# Return the point on a line segment which is closest to a given point. 377# 378# Arguments: 379# P a point 380# linesegment a line segment 381# 382# Results: 383# Q the point on the line segment that has the 384# smallest distance to P 385# 386# Examples: 387# - findClosestPointOnLineSegment {5 10} {0 0 10 10} 388# Result: 7.5 7.5 389# - findClosestPointOnLineSegment {-10 0} {0 0 10 10} 390# Result: 0 0 391# 392proc ::math::geometry::findClosestPointOnLineSegment {P linesegment} { 393 set result [findClosestPointOnLineImpl $P $linesegment] 394 set Q [lindex $result 0] 395 set r [lindex $result 1] 396 if {$r<0} { 397 return [lrange $linesegment 0 1] 398 } elseif {$r>1} { 399 return [lrange $linesegment 2 3] 400 } else { 401 return $Q 402 } 403 404} 405 406# ::math::geometry::calculateDistanceToPolyline 407# 408# Calculate the distance between a point and a polyline. 409# 410# Arguments: 411# P a point 412# polyline a polyline 413# 414# Results: 415# dist the smallest distance between P and any point 416# on the polyline 417# 418# Examples: 419# - calculateDistanceToPolyline {10 10} {0 0 10 5 20 0} 420# Result: 5.0 421# - calculateDistanceToPolyline {5 10} {0 0 10 5 20 0} 422# Result: 6.7082039325 423# 424proc ::math::geometry::calculateDistanceToPolyline {P polyline} { 425 set minDist "none" 426 foreach {Ax Ay} [lrange $polyline 0 end-2] {Bx By} [lrange $polyline 2 end] { 427 set dist [calculateDistanceToLineSegment $P [list $Ax $Ay $Bx $By]] 428 if {$minDist=="none" || $dist < $minDist} { 429 set minDist $dist 430 } 431 } 432 return $minDist 433} 434 435# ::math::geometry::findClosestPointOnPolyline 436# 437# Return the point on a polyline which is closest to a given point. 438# 439# Arguments: 440# P a point 441# polyline a polyline 442# 443# Results: 444# Q the point on the polyline that has the smallest 445# distance to P 446# 447# Examples: 448# - findClosestPointOnPolyline {10 10} {0 0 10 5 20 0} 449# Result: 10 5 450# - findClosestPointOnPolyline {5 10} {0 0 10 5 20 0} 451# Result: 8.0 4.0 452# 453proc ::math::geometry::findClosestPointOnPolyline {P polyline} { 454 set closestPoint "none" 455 foreach {Ax Ay} [lrange $polyline 0 end-2] {Bx By} [lrange $polyline 2 end] { 456 set Q [findClosestPointOnLineSegment $P [list $Ax $Ay $Bx $By]] 457 set dist [lengthOfPolyline [concat $P $Q]] 458 if {$closestPoint=="none" || $dist<$closestDistance} { 459 set closestPoint $Q 460 set closestDistance $dist 461 } 462 } 463 return $closestPoint 464} 465 466 467 468 469 470 471# ::math::geometry::lengthOfPolyline 472# 473# Find the length of a polyline, i.e., the sum of the 474# lengths of the individual line segments. 475# 476# Arguments: 477# polyline a polyline 478# 479# Results: 480# length the length of the polyline 481# 482# Examples: 483# - lengthOfPolyline {0 0 5 0 5 10} 484# Result: 15.0 485# 486proc ::math::geometry::lengthOfPolyline {polyline} { 487 set length 0 488 foreach {x1 y1} [lrange $polyline 0 end-2] {x2 y2} [lrange $polyline 2 end] { 489 set length [expr {$length + sqrt(pow($x1-$x2,2) + pow($y1-$y2,2))}] 490 #set length [expr {$length + sqrt(($x1-$x2)*($x1-$x2) + ($y1-$y2)*($y1-$y2))}] 491 } 492 return $length 493} 494 495 496 497 498# ::math::geometry::movePointInDirection 499# 500# Move a point in a given direction. 501# 502# Arguments: 503# P the starting point 504# direction the direction from P 505# The direction is in 360-degrees going counter-clockwise, 506# with "straight right" being 0 degrees 507# dist the distance from P 508# 509# Results: 510# Q the point which is found by starting in P and going 511# in the given direction, until the distance between 512# P and Q is dist 513# 514# Examples: 515# - movePointInDirection {0 0} 45.0 10 516# Result: 7.07106781187 7.07106781187 517# 518proc ::math::geometry::movePointInDirection {P direction dist} { 519 set x [lindex $P 0] 520 set y [lindex $P 1] 521 set pi [expr {4*atan(1)}] 522 set xt [expr {$x + $dist*cos(($direction*$pi)/180)}] 523 set yt [expr {$y + $dist*sin(($direction*$pi)/180)}] 524 return [list $xt $yt] 525} 526 527 528# ::math::geometry::angle 529# 530# Calculates angle from the horizon (0,0)->(1,0) to a line. 531# 532# Arguments: 533# line a line defined by two points A and B 534# 535# Results: 536# angle the angle between the line (0,0)->(1,0) and (Ax,Ay)->(Bx,By). 537# Angle is in 360-degrees going counter-clockwise 538# 539# Examples: 540# - angle {10 10 15 13} 541# Result: 30.9637565321 542# 543proc ::math::geometry::angle {line} { 544 set x1 [lindex $line 0] 545 set y1 [lindex $line 1] 546 set x2 [lindex $line 2] 547 set y2 [lindex $line 3] 548 # - handle vertical lines 549 if {$x1==$x2} {if {$y1<$y2} {return 90} else {return 270}} 550 # - handle other lines 551 set a [expr {atan(abs((1.0*$y1-$y2)/(1.0*$x1-$x2)))}] ; # a is between 0 and pi/2 552 set pi [expr {4*atan(1)}] 553 if {$y1<=$y2} { 554 # line is going upwards 555 if {$x1<$x2} {set b $a} else {set b [expr {$pi-$a}]} 556 } else { 557 # line is going downwards 558 if {$x1<$x2} {set b [expr {2*$pi-$a}]} else {set b [expr {$pi+$a}]} 559 } 560 return [expr {$b/$pi*180}] ; # convert b to degrees 561} 562 563 564 565 566### 567# 568# Intersection procedures 569# 570### 571 572# ::math::geometry::lineSegmentsIntersect 573# 574# Checks whether two line segments intersect. 575# 576# Arguments: 577# linesegment1 the first line segment 578# linesegment2 the second line segment 579# 580# Results: 581# dointersect a boolean saying whether the line segments intersect 582# (i.e., have any points in common) 583# 584# Examples: 585# - lineSegmentsIntersect {0 0 10 10} {0 10 10 0} 586# Result: 1 587# - lineSegmentsIntersect {0 0 10 10} {20 20 20 30} 588# Result: 0 589# - lineSegmentsIntersect {0 0 10 10} {10 10 15 15} 590# Result: 1 591# 592proc ::math::geometry::lineSegmentsIntersect {linesegment1 linesegment2} { 593 # Algorithm based on Sedgewick. 594 set l1x1 [lindex $linesegment1 0] 595 set l1y1 [lindex $linesegment1 1] 596 set l1x2 [lindex $linesegment1 2] 597 set l1y2 [lindex $linesegment1 3] 598 set l2x1 [lindex $linesegment2 0] 599 set l2y1 [lindex $linesegment2 1] 600 set l2x2 [lindex $linesegment2 2] 601 set l2y2 [lindex $linesegment2 3] 602 return [expr {([ccw [list $l1x1 $l1y1] [list $l1x2 $l1y2] [list $l2x1 $l2y1]]\ 603 *[ccw [list $l1x1 $l1y1] [list $l1x2 $l1y2] [list $l2x2 $l2y2]] <= 0) \ 604 && ([ccw [list $l2x1 $l2y1] [list $l2x2 $l2y2] [list $l1x1 $l1y1]]\ 605 *[ccw [list $l2x1 $l2y1] [list $l2x2 $l2y2] [list $l1x2 $l1y2]] <= 0)}] 606} 607 608# ::math::geometry::findLineSegmentIntersection 609# 610# Returns the intersection point of two line segments. 611# Note: may also return "coincident" and "none". 612# 613# Arguments: 614# linesegment1 the first line segment 615# linesegment2 the second line segment 616# 617# Results: 618# P the intersection point of linesegment1 and linesegment2. 619# If linesegment1 and linesegment2 have an infinite number 620# of points in common, the procedure returns "coincident". 621# If there are no intersection points, the procedure 622# returns "none". 623# 624# Examples: 625# - findLineSegmentIntersection {0 0 10 10} {0 10 10 0} 626# Result: 5.0 5.0 627# - findLineSegmentIntersection {0 0 10 10} {20 20 20 30} 628# Result: none 629# - findLineSegmentIntersection {0 0 10 10} {10 10 15 15} 630# Result: 10.0 10.0 631# - findLineSegmentIntersection {0 0 10 10} {5 5 15 15} 632# Result: coincident 633# 634proc ::math::geometry::findLineSegmentIntersection {linesegment1 linesegment2} { 635 if {[lineSegmentsIntersect $linesegment1 $linesegment2]} { 636 set lineintersect [findLineIntersection $linesegment1 $linesegment2] 637 switch -- $lineintersect { 638 639 "coincident" { 640 # lines are coincident 641 set l1x1 [lindex $linesegment1 0] 642 set l1y1 [lindex $linesegment1 1] 643 set l1x2 [lindex $linesegment1 2] 644 set l1y2 [lindex $linesegment1 3] 645 set l2x1 [lindex $linesegment2 0] 646 set l2y1 [lindex $linesegment2 1] 647 set l2x2 [lindex $linesegment2 2] 648 set l2y2 [lindex $linesegment2 3] 649 # check if the line SEGMENTS overlap 650 # (NOT enough to check if the x-intervals overlap (vertical lines!)) 651 set overlapx [intervalsOverlap $l1x1 $l1x2 $l2x1 $l2x2 0] 652 set overlapy [intervalsOverlap $l1y1 $l1y2 $l2y1 $l2y2 0] 653 if {$overlapx && $overlapy} { 654 return "coincident" 655 } else { 656 return "none" 657 } 658 } 659 660 "none" { 661 # should never happen, because we call "lineSegmentsIntersect" first 662 puts stderr "::math::geometry::findLineSegmentIntersection: suddenly no intersection?" 663 return "none" 664 } 665 666 default { 667 # lineintersect = the intersection point 668 return $lineintersect 669 } 670 } 671 } else { 672 return "none" 673 } 674} 675 676# ::math::geometry::findLineIntersection {line1 line2} 677# 678# Returns the intersection point of two lines. 679# Note: may also return "coincident" and "none". 680# 681# Arguments: 682# line1 the first line 683# line2 the second line 684# 685# Results: 686# P the intersection point of line1 and line2. 687# If line1 and line2 have an infinite number of points 688# in common, the procedure returns "coincident". 689# If there are no intersection points, the procedure 690# returns "none". 691# 692# Examples: 693# - findLineIntersection {0 0 10 10} {0 10 10 0} 694# Result: 5.0 5.0 695# - findLineIntersection {0 0 10 10} {20 20 20 30} 696# Result: 20.0 20.0 697# - findLineIntersection {0 0 10 10} {10 10 15 15} 698# Result: coincident 699# - findLineIntersection {0 0 10 10} {5 5 15 15} 700# Result: coincident 701# - findLineIntersection {0 0 10 10} {0 1 10 11} 702# Result: none 703# 704proc ::math::geometry::findLineIntersection {line1 line2} { 705 706 # References: 707 # http://wiki.tcl.tk/12070 (Kevin Kenny) 708 # http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/ 709 710 set l1x1 [lindex $line1 0] 711 set l1y1 [lindex $line1 1] 712 set l1x2 [lindex $line1 2] 713 set l1y2 [lindex $line1 3] 714 715 set l2x1 [lindex $line2 0] 716 set l2y1 [lindex $line2 1] 717 set l2x2 [lindex $line2 2] 718 set l2y2 [lindex $line2 3] 719 720 set d [expr {($l2y2 - $l2y1) * ($l1x2 - $l1x1) - 721 ($l2x2 - $l2x1) * ($l1y2 - $l1y1)}] 722 set na [expr {($l2x2 - $l2x1) * ($l1y1 - $l2y1) - 723 ($l2y2 - $l2y1) * ($l1x1 - $l2x1)}] 724 725 # http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/ 726 if {$d == 0} { 727 if {$na == 0} { 728 return "coincident" 729 } else { 730 return "none" 731 } 732 } 733 set r [list \ 734 [expr {$l1x1 + $na * ($l1x2 - $l1x1) / $d}] \ 735 [expr {$l1y1 + $na * ($l1y2 - $l1y1) / $d}]] 736 return $r 737} 738 739 740# ::math::geometry::polylinesIntersect 741# 742# Checks whether two polylines intersect. 743# 744# Arguments; 745# polyline1 the first polyline 746# polyline2 the second polyline 747# 748# Results: 749# dointersect a boolean saying whether the polylines intersect 750# 751# Examples: 752# - polylinesIntersect {0 0 10 10 10 20} {0 10 10 0} 753# Result: 1 754# - polylinesIntersect {0 0 10 10 10 20} {5 4 10 4} 755# Result: 0 756# 757proc ::math::geometry::polylinesIntersect {polyline1 polyline2} { 758 return [polylinesBoundingIntersect $polyline1 $polyline2 0] 759} 760 761# ::math::geometry::polylinesBoundingIntersect 762# 763# Check whether two polylines intersect, but reduce 764# the correctness of the result to the given granularity. 765# Use this for faster, but weaker, intersection checking. 766# 767# How it works: 768# Each polyline is split into a number of smaller polylines, 769# consisting of granularity points each. If a pair of those smaller 770# lines' bounding boxes intersect, then this procedure returns 1, 771# otherwise it returns 0. 772# 773# Arguments: 774# polyline1 the first polyline 775# polyline2 the second polyline 776# granularity the number of points in each part-polyline 777# granularity<=1 means full correctness 778# 779# Results: 780# dointersect a boolean saying whether the polylines intersect 781# 782# Examples: 783# - polylinesBoundingIntersect {0 0 10 10 10 20} {0 10 10 0} 2 784# Result: 1 785# - polylinesBoundingIntersect {0 0 10 10 10 20} {5 4 10 4} 2 786# Result: 1 787# 788proc ::math::geometry::polylinesBoundingIntersect {polyline1 polyline2 granularity} { 789 if {$granularity<=1} { 790 # Use perfect intersect 791 # => first pin down where an intersection point may be, and then 792 # call MultilinesIntersectPerfect on those parts 793 set granularity 10 ; # optimal search granularity? 794 set perfectmatch 1 795 } else { 796 set perfectmatch 0 797 } 798 799 # split the lines into parts consisting of $granularity points 800 set polyline1parts {} 801 for {set i 0} {$i<[llength $polyline1]} {incr i [expr {2*$granularity-2}]} { 802 lappend polyline1parts [lrange $polyline1 $i [expr {$i+2*$granularity-1}]] 803 } 804 set polyline2parts {} 805 for {set i 0} {$i<[llength $polyline2]} {incr i [expr {2*$granularity-2}]} { 806 lappend polyline2parts [lrange $polyline2 $i [expr {$i+2*$granularity-1}]] 807 } 808 809 # do any of the parts overlap? 810 foreach part1 $polyline1parts { 811 foreach part2 $polyline2parts { 812 set part1bbox [bbox $part1] 813 set part2bbox [bbox $part2] 814 if {[rectanglesOverlap [lrange $part1bbox 0 1] [lrange $part1bbox 2 3] \ 815 [lrange $part2bbox 0 1] [lrange $part2bbox 2 3] 0]} { 816 # the lines' bounding boxes intersect 817 if {$perfectmatch} { 818 foreach {l1x1 l1y1} [lrange $part1 0 end-2] {l1x2 l1y2} [lrange $part1 2 end] { 819 foreach {l2x1 l2y1} [lrange $part2 0 end-2] {l2x2 l2y2} [lrange $part2 2 end] { 820 if {[lineSegmentsIntersect [list $l1x1 $l1y1 $l1x2 $l1y2] \ 821 [list $l2x1 $l2y1 $l2x2 $l2y2]]} { 822 # two line segments overlap 823 return 1 824 } 825 } 826 } 827 return 0 828 } else { 829 return 1 830 } 831 } 832 } 833 } 834 return 0 835} 836 837# ::math::geometry::ccw 838# 839# PRIVATE FUNCTION USED BY OTHER FUNCTIONS. 840# Returns whether traversing from A to B to C is CounterClockWise 841# Algorithm by Sedgewick. 842# 843# Arguments: 844# A first point 845# B second point 846# C third point 847# 848# Reeults: 849# ccw a boolean saying whether traversing from A to B to C 850# is CounterClockWise 851# 852proc ::math::geometry::ccw {A B C} { 853 set Ax [lindex $A 0] 854 set Ay [lindex $A 1] 855 set Bx [lindex $B 0] 856 set By [lindex $B 1] 857 set Cx [lindex $C 0] 858 set Cy [lindex $C 1] 859 set dx1 [expr {$Bx - $Ax}] 860 set dy1 [expr {$By - $Ay}] 861 set dx2 [expr {$Cx - $Ax}] 862 set dy2 [expr {$Cy - $Ay}] 863 if {$dx1*$dy2 > $dy1*$dx2} {return 1} 864 if {$dx1*$dy2 < $dy1*$dx2} {return -1} 865 if {($dx1*$dx2 < 0) || ($dy1*$dy2 < 0)} {return -1} 866 if {($dx1*$dx1 + $dy1*$dy1) < ($dx2*$dx2+$dy2*$dy2)} {return 1} 867 return 0 868} 869 870 871 872 873 874 875 876### 877# 878# Overlap procedures 879# 880### 881 882# ::math::geometry::intervalsOverlap 883# 884# Check whether two intervals overlap. 885# Examples: 886# - (2,4) and (5,3) overlap with strict=0 and strict=1 887# - (2,4) and (1,2) overlap with strict=0 but not with strict=1 888# 889# Arguments: 890# y1,y2 the first interval 891# y3,y4 the second interval 892# strict choosing strict or non-strict interpretation 893# 894# Results: 895# dooverlap a boolean saying whether the intervals overlap 896# 897# Examples: 898# - intervalsOverlap 2 4 4 6 1 899# Result: 0 900# - intervalsOverlap 2 4 4 6 0 901# Result: 1 902# - intervalsOverlap 4 2 3 5 0 903# Result: 1 904# 905proc ::math::geometry::intervalsOverlap {y1 y2 y3 y4 strict} { 906 if {$y1>$y2} { 907 set temp $y1 908 set y1 $y2 909 set y2 $temp 910 } 911 if {$y3>$y4} { 912 set temp $y3 913 set y3 $y4 914 set y4 $temp 915 } 916 if {$strict} { 917 return [expr {$y2>$y3 && $y4>$y1}] 918 } else { 919 return [expr {$y2>=$y3 && $y4>=$y1}] 920 } 921} 922 923# ::math::geometry::rectanglesOverlap 924# 925# Check whether two rectangles overlap (see also intervalsOverlap). 926# 927# Arguments: 928# P1 upper-left corner of the first rectangle 929# P2 lower-right corner of the first rectangle 930# Q1 upper-left corner of the second rectangle 931# Q2 lower-right corner of the second rectangle 932# strict choosing strict or non-strict interpretation 933# 934# Results: 935# dooverlap a boolean saying whether the rectangles overlap 936# 937# Examples: 938# - rectanglesOverlap {0 10} {10 0} {10 10} {20 0} 1 939# Result: 0 940# - rectanglesOverlap {0 10} {10 0} {10 10} {20 0} 0 941# Result: 1 942# 943proc ::math::geometry::rectanglesOverlap {P1 P2 Q1 Q2 strict} { 944 set b1x1 [lindex $P1 0] 945 set b1y1 [lindex $P1 1] 946 set b1x2 [lindex $P2 0] 947 set b1y2 [lindex $P2 1] 948 set b2x1 [lindex $Q1 0] 949 set b2y1 [lindex $Q1 1] 950 set b2x2 [lindex $Q2 0] 951 set b2y2 [lindex $Q2 1] 952 # ensure b1x1<=b1x2 etc. 953 if {$b1x1 > $b1x2} { 954 set temp $b1x1 955 set b1x1 $b1x2 956 set b1x2 $temp 957 } 958 if {$b1y1 > $b1y2} { 959 set temp $b1y1 960 set b1y1 $b1y2 961 set b1y2 $temp 962 } 963 if {$b2x1 > $b2x2} { 964 set temp $b2x1 965 set b2x1 $b2x2 966 set b2x2 $temp 967 } 968 if {$b2y1 > $b2y2} { 969 set temp $b2y1 970 set b2y1 $b2y2 971 set b2y2 $temp 972 } 973 # Check if the boxes intersect 974 # (From: Cormen, Leiserson, and Rivests' "Algorithms", page 889) 975 if {$strict} { 976 return [expr {($b1x2>$b2x1) && ($b2x2>$b1x1) \ 977 && ($b1y2>$b2y1) && ($b2y2>$b1y1)}] 978 } else { 979 return [expr {($b1x2>=$b2x1) && ($b2x2>=$b1x1) \ 980 && ($b1y2>=$b2y1) && ($b2y2>=$b1y1)}] 981 } 982} 983 984 985 986# ::math::geometry::bbox 987# 988# Calculate the bounding box of a polyline. 989# 990# Arguments: 991# polyline a polyline 992# 993# Results: 994# x1,y1,x2,y2 four coordinates where (x1,y1) is the upper-left corner 995# of the bounding box, and (x2,y2) is the lower-right corner 996# 997# Examples: 998# - bbox {0 10 4 1 6 23 -12 5} 999# Result: -12 1 6 23 1000# 1001proc ::math::geometry::bbox {polyline} { 1002 set minX [lindex $polyline 0] 1003 set maxX $minX 1004 set minY [lindex $polyline 1] 1005 set maxY $minY 1006 foreach {x y} $polyline { 1007 if {$x < $minX} {set minX $x} 1008 if {$x > $maxX} {set maxX $x} 1009 if {$y < $minY} {set minY $y} 1010 if {$y > $maxY} {set maxY $y} 1011 } 1012 return [list $minX $minY $maxX $maxY] 1013} 1014 1015# ::math::geometry::ClosedPolygon 1016# 1017# Return a closed polygon - used internally 1018# 1019# Arguments: 1020# polygon a polygon 1021# 1022# Results: 1023# closedpolygon a polygon whose first and last vertices 1024# coincide 1025# 1026proc ::math::geometry::ClosedPolygon {polygon} { 1027 1028 if { [lindex $polygon 0] != [lindex $polygon end-1] || 1029 [lindex $polygon 1] != [lindex $polygon end] } { 1030 1031 return [concat $polygon [lrange $polygon 0 1]] 1032 1033 } else { 1034 1035 return $polygon 1036 } 1037} 1038 1039 1040# ::math::geometry::pointInsidePolygon 1041# 1042# Determine if a point is completely inside a polygon. If the point 1043# touches the polygon, then the point is not complete inside the 1044# polygon. 1045# 1046# Arguments: 1047# P a point 1048# polygon a polygon 1049# 1050# Results: 1051# isinside a boolean saying whether the point is 1052# completely inside the polygon or not 1053# 1054# Examples: 1055# - pointInsidePolygon {5 5} {4 4 4 6 6 6 6 4} 1056# Result: 1 1057# - pointInsidePolygon {5 5} {6 6 6 7 7 7} 1058# Result: 0 1059# 1060proc ::math::geometry::pointInsidePolygon {P polygon} { 1061 # check if P is on one of the polygon's sides (if so, P is not 1062 # inside the polygon) 1063 set closedPolygon [ClosedPolygon $polygon] 1064 1065 foreach {x1 y1} [lrange $closedPolygon 0 end-2] {x2 y2} [lrange $closedPolygon 2 end] { 1066 if {[calculateDistanceToLineSegment $P [list $x1 $y1 $x2 $y2]]<0.0000001} { 1067 return 0 1068 } 1069 } 1070 1071 # Algorithm 1072 # 1073 # Consider a straight line going from P to a point far away from both 1074 # P and the polygon (in particular outside the polygon). 1075 # - If the line intersects with 0 of the polygon's sides, then 1076 # P must be outside the polygon. 1077 # - If the line intersects with 1 of the polygon's sides, then 1078 # P must be inside the polygon (since the other end of the line 1079 # is outside the polygon). 1080 # - If the line intersects with 2 of the polygon's sides, then 1081 # the line must pass into one polygon area and out of it again, 1082 # and hence P is outside the polygon. 1083 # - In general: if the line intersects with the polygon's sides an odd 1084 # number of times, then P is inside the polygon. Note: we also have 1085 # to check whether the line crosses one of the polygon's 1086 # bend points for the same reason. 1087 1088 # get point far away and define the line 1089 set polygonBbox [bbox $polygon] 1090 1091 set pointFarAway [list \ 1092 [expr {[lindex $polygonBbox 0]-[lindex $polygonBbox 2]}] \ 1093 [expr {[lindex $polygonBbox 1]-0.1*[lindex $polygonBbox 3]}]] 1094 1095 set infinityLine [concat $pointFarAway $P] 1096 # calculate number of intersections 1097 set noOfIntersections 0 1098 # 1. count intersections between the line and the polygon's sides 1099 foreach {x1 y1} [lrange $closedPolygon 0 end-2] {x2 y2} [lrange $closedPolygon 2 end] { 1100 if {[lineSegmentsIntersect $infinityLine [list $x1 $y1 $x2 $y2]]} { 1101 incr noOfIntersections 1102 } 1103 } 1104 # 2. count intersections between the line and the polygon's points 1105 foreach {x1 y1} $closedPolygon { 1106 if {[calculateDistanceToLineSegment [list $x1 $y1] $infinityLine]<0.0000001} { 1107 incr noOfIntersections 1108 } 1109 } 1110 return [expr {$noOfIntersections % 2}] 1111} 1112 1113 1114# ::math::geometry::rectangleInsidePolygon 1115# 1116# Determine if a rectangle is completely inside a polygon. If polygon 1117# touches the rectangle, then the rectangle is not complete inside the 1118# polygon. 1119# 1120# Arguments: 1121# P1 upper-left corner of the rectangle 1122# P2 lower-right corner of the rectangle 1123# polygon a polygon 1124# 1125# Results: 1126# isinside a boolean saying whether the rectangle is 1127# completely inside the polygon or not 1128# 1129# Examples: 1130# - rectangleInsidePolygon {0 10} {10 0} {-10 -10 0 11 11 11 11 0} 1131# Result: 1 1132# - rectangleInsidePolygon {0 0} {0 0} {-16 14 5 -16 -16 -25 -21 16 -19 24} 1133# Result: 1 1134# - rectangleInsidePolygon {0 0} {0 0} {2 2 2 4 4 4 4 2} 1135# Result: 0 1136# 1137proc ::math::geometry::rectangleInsidePolygon {P1 P2 polygon} { 1138 # get coordinates of rectangle 1139 set bx1 [lindex $P1 0] 1140 set by1 [lindex $P1 1] 1141 set bx2 [lindex $P2 0] 1142 set by2 [lindex $P2 1] 1143 1144 # if rectangle does not overlap with the bbox of polygon, then the 1145 # rectangle cannot be inside the polygon (this is a quick way to 1146 # get an answer in many cases) 1147 set polygonBbox [bbox $polygon] 1148 set polygonP1x [lindex $polygonBbox 0] 1149 set polygonP1y [lindex $polygonBbox 1] 1150 set polygonP2x [lindex $polygonBbox 2] 1151 set polygonP2y [lindex $polygonBbox 3] 1152 if {![rectanglesOverlap [list $bx1 $by1] [list $bx2 $by2] \ 1153 [list $polygonP1x $polygonP1y] [list $polygonP2x $polygonP2y] 0]} { 1154 return 0 1155 } 1156 1157 # 1. if one of the points of the polygon is inside the rectangle, 1158 # then the rectangle cannot be inside the polygon 1159 foreach {x y} $polygon { 1160 if {$bx1<$x && $x<$bx2 && $by1<$y && $y<$by2} { 1161 return 0 1162 } 1163 } 1164 1165 # 2. if one of the line segments of the polygon intersect with the 1166 # rectangle, then the rectangle cannot be inside the polygon 1167 set rectanglePolyline [list $bx1 $by1 $bx2 $by1 $bx2 $by2 $bx1 $by2 $bx1 $by1] 1168 set closedPolygon [ClosedPolygon $polygon] 1169 if {[polylinesIntersect $closedPolygon $rectanglePolyline]} { 1170 return 0 1171 } 1172 1173 # at this point we know that: 1174 # 1. the polygon has no points inside the rectangle 1175 # 2. the polygon's sides don't intersect with the rectangle 1176 # therefore: 1177 # either the rectangle is (completely) inside the polygon, or 1178 # the rectangle is (completely) outside the polygon 1179 1180 # final test: if one of the points on the rectangle is inside the 1181 # polygon, then the whole rectangle must be inside the rectangle 1182 return [pointInsidePolygon [list $bx1 $by1] $polygon] 1183} 1184 1185 1186# ::math::geometry::areaPolygon 1187# 1188# Determine the area enclosed by a (non-complex) polygon 1189# 1190# Arguments: 1191# polygon a polygon 1192# 1193# Results: 1194# area the area enclosed by the polygon 1195# 1196# Examples: 1197# - areaPolygon {-10 -10 10 -10 10 10 -10 10} 1198# Result: 400 1199# 1200proc ::math::geometry::areaPolygon {polygon} { 1201 1202 foreach {a1 a2 b1 b2} $polygon {break} 1203 1204 set area 0.0 1205 foreach {c1 c2} [lrange $polygon 4 end] { 1206 set area [expr {$area + $b1*$c2 - $b2*$c1}] 1207 set b1 $c1 1208 set b2 $c2 1209 } 1210 expr {0.5*abs($area)} 1211} 1212 1213# # ## ### ##### ############# 1214 1215namespace eval ::math::geometry { 1216 variable pi [expr { 4 * atan(1) }] 1217 variable torad [expr { (4 * atan(1)) / 180.0 }] 1218 variable todeg [expr { 180.0 / (4 * atan(1)) }] 1219 1220 namespace export \ 1221 + - s* direction v h p between distance length \ 1222 nwse rect octant findLineSegmentIntersection \ 1223 findLineIntersection bbox x y conjx conjy 1224} 1225 1226package provide math::geometry 1.1.2 1227