Tcl Library Source Code

Check-in [c76dc6a3fc]
Login
Bounty program for improvements to Tcl and certain Tcl packages.
Tcl 2019 Conference, Houston/TX, US, Nov 4-8
Send your abstracts to [email protected]
or submit via the online form by Sep 9.

Many hyperlinks are disabled.
Use anonymous login to enable hyperlinks.

Overview
Comment:Merged fix for ticket [cb043ecc70e0e90bf].
Timelines: family | ancestors | descendants | both | trunk
Files: files | file ages | folders
SHA1: c76dc6a3fca3b337b144fc3447410fa60cc68e0a
User & Date: aku 2017-08-07 23:41:56
References
2017-08-07
23:42 Closed ticket [cb043ecc70]: areaPolygon miscalculates area plus 6 other changes artifact: a6be40a4e2 user: aku
Context
2017-08-07
23:54
Merged fix for ticket [85a0ec8f50]. check-in: 98a32c3965 user: aku tags: trunk
23:53
Fix infinite loop when splitx is invoked with a regexp matching the empty string. Bumped version. Extended testsuite. Closed-Leaf check-in: 31bfe3ad8e user: aku tags: tkt-85a0ec8f50-ak
23:41
Merged fix for ticket [cb043ecc70e0e90bf]. check-in: c76dc6a3fc user: aku tags: trunk
23:40
Fixed missing version bump for bugfix. Fixed warning in statistics documentation. Closed-Leaf check-in: bf6666afc5 user: aku tags: avl-fix-areaPolygon-and-more
2017-06-05
21:00
Merged fixes for tkt [71deadcf96]. check-in: 78d7722e1e user: aku tags: trunk
Changes
Hide Diffs Unified Diffs Ignore Whitespace Patch

Changes to modules/math/geometry.tcl.

10
11
12
13
14
15
16

17
18
19
20
21
22
23
..
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
...
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
...
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
...
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
...
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
...
263
264
265
266
267
268
269
270

271
272
273
274
275
276
277
278
...
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
...
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
...
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
...
418
419
420
421
422
423
424
425
426
427
428
429
430

431
432
433
434
435
436
437
...
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481

482
483
484
485
486
487
488
...
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
...
849
850
851
852
853
854
855
856
857
858
859
860
861
862

863

864
865
866
867
868
869
870
....
1059
1060
1061
1062
1063
1064
1065

1066
1067
1068
1069

1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
....
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106

1107
1108
1109
1110
1111
1112
1113
....
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141

1142
1143
1144
1145
1146
1147
1148
....
1317
1318
1319
1320
1321
1322
1323
1324

1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
....
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
....
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
....
1550
1551
1552
1553
1554
1555
1556
1557
# See the file "license.terms" for information on usage and redistribution
# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
# RCS: @(#) $Id: geometry.tcl,v 1.12 2010/05/24 21:44:16 andreas_kupries Exp $

namespace eval ::math::geometry {}


package require math

###
#
# POINTS
#
#    A point P consists of an x-coordinate, Px, and a y-coordinate, Py,
................................................................................
# Point constructor
proc ::math::geometry::p {x y} {
    return [list $x $y]
}

# Vector addition
proc ::math::geometry::+ {pa pb} {
    foreach {ax ay} $pa break
    foreach {bx by} $pb break
    return [list [expr {$ax + $bx}] [expr {$ay + $by}]]
}

# Vector difference
proc ::math::geometry::- {pa pb} {
    foreach {ax ay} $pa break
    foreach {bx by} $pb break
    return [list [expr {$ax - $bx}] [expr {$ay - $by}]]
}

# Distance between 2 points
proc ::math::geometry::distance {pa pb} {
    foreach {ax ay} $pa break
    foreach {bx by} $pb break
    return [expr {hypot($bx-$ax,$by-$ay)}]
}

# Length of a vector
proc ::math::geometry::length {v} {
    foreach {x y} $v break
    return [expr {hypot($x,$y)}]
}

# Scaling a vector by a factor
proc ::math::geometry::s* {factor p} {
    foreach {x y} $p break
    return [list [expr {$x * $factor}] [expr {$y * $factor}]]
}

# Unit vector into specific direction given by angle (degrees)
proc ::math::geometry::direction {angle} {
    variable torad
    set x [expr {  cos($angle * $torad)}]
................................................................................
proc ::math::geometry::between {pa pb s} {
    return [+ $pa [s* $s [- $pb $pa]]]
}

# Find direction octant the point (vector) lies in.
proc ::math::geometry::octant {p} {
    variable todeg
    foreach {x y} $p break

    set a [expr {(atan2(-$y,$x)*$todeg)}]
    while {$a >  360} {set a [expr {$a - 360}]}
    while {$a < -360} {set a [expr {$a + 360}]}
    if {$a < 0} {set a [expr {360 + $a}]}

    #puts "p ($x, $y) @ angle $a | [expr {atan2($y,$x)}] | [expr {atan2($y,$x)*$todeg}]"
................................................................................
	}
	return east ; # a <= 360.0
    }
}

# Return the NW and SE corners of the rectangle.
proc ::math::geometry::nwse {rect} {
    foreach {xnw ynw xse yse} $rect break
    return [list [p $xnw $ynw] [p $xse $yse]]
}

# Construct rectangle from NW and SE corners.
proc ::math::geometry::rect {pa pb} {
    foreach {ax ay} $pa break
    foreach {bx by} $pb break
    return [list $ax $ay $bx $by]
}

proc ::math::geometry::conjx {p} {
    foreach {x y} $p break
    return [list [expr {- $x}] $y]
}

proc ::math::geometry::conjy {p} {
    foreach {x y} $p break
    return [list $x [expr {- $y}]]
}

proc ::math::geometry::x {p} {
    foreach {x y} $p break
    return $x
}

proc ::math::geometry::y {p} {
    foreach {x y} $p break
    return $y
}

# ::math::geometry::calculateDistanceToLine
#
#       Calculate the distance between a point and a line.
#
# Arguments:
................................................................................
#     - calculateDistanceToLine {5 10} {0 0 10 10}
#       Result: 3.53553390593
#     - calculateDistanceToLine {-10 0} {0 0 10 10}
#       Result: 7.07106781187
#
proc ::math::geometry::calculateDistanceToLine {P line} {
    # solution based on FAQ 1.02 on comp.graphics.algorithms
    # L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
    #     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
    # s = -----------------------------
    #                 L^2
    # dist = |s|*L
    #
    # =>
    #
................................................................................
    set Bx [lindex $line 2]
    set By [lindex $line 3]
    set Cx [lindex $P 0]
    set Cy [lindex $P 1]
    if {$Ax==$Bx && $Ay==$By} {
	return [lengthOfPolyline [concat $P [lrange $line 0 1]]]
    } else {
	set L [expr {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}]
	return [expr {abs(($Ay-$Cy)*($Bx-$Ax)-($Ax-$Cx)*($By-$Ay)) / $L}]
    }
}

# ::math::geometry::findClosestPointOnLine
#
#       Return the point on a line which is closest to a given point.
................................................................................
#                        r=0      P = A
#                        r=1      P = B
#                        r<0      P is on the backward extension of AB
#                        r>1      P is on the forward extension of AB
#                        0<r<1    P is interior to AB
#
proc ::math::geometry::findClosestPointOnLineImpl {P line} {
    # solution based on FAQ 1.02 on comp.graphics.algorithms

    #   L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
    #        (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
    #   r = -------------------------------
    #                     L^2
    #   Px = Ax + r(Bx-Ax)
    #   Py = Ay + r(By-Ay)
    set Ax [lindex $line 0]
    set Ay [lindex $line 1]
................................................................................
    set Bx [lindex $line 2]
    set By [lindex $line 3]
    set Cx [lindex $P 0]
    set Cy [lindex $P 1]
    if {$Ax==$Bx && $Ay==$By} {
	return [list [list $Ax $Ay] 0]
    } else {
	set L [expr {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}]
	set r [expr {(($Cx-$Ax)*($Bx-$Ax) + ($Cy-$Ay)*($By-$Ay))/pow($L,2)}]
	set Px [expr {$Ax + $r*($Bx-$Ax)}]
	set Py [expr {$Ay + $r*($By-$Ay)}]
	return [list [list $Px $Py] $r]
    }
}

# ::math::geometry::calculateDistanceToLineSegment
................................................................................
#                        r=1      P = B
#                        r<0      P is on the backward extension of AB
#                        r>1      P is on the forward extension of AB
#                        0<r<1    P is interior to AB
#
proc ::math::geometry::calculateDistanceToLineSegmentImpl {P linesegment} {
    # solution based on FAQ 1.02 on comp.graphics.algorithms
    # L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
    #     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
    # s = -----------------------------
    #                 L^2
    #      (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
    # r = -------------------------------
    #                   L^2
    # dist = |s|*L
................................................................................
    set Bx [lindex $linesegment 2]
    set By [lindex $linesegment 3]
    set Cx [lindex $P 0]
    set Cy [lindex $P 1]
    if {$Ax==$Bx && $Ay==$By} {
	return [list [lengthOfPolyline [concat $P [lrange $linesegment 0 1]]] 0]
    } else {
	set L [expr {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}]
	set r [expr {(($Cx-$Ax)*($Bx-$Ax) + ($Cy-$Ay)*($By-$Ay))/pow($L,2)}]
	return [list [expr {abs(($Ay-$Cy)*($Bx-$Ax)-($Ax-$Cx)*($By-$Ay)) / $L}] $r]
    }
}

# ::math::geometry::findClosestPointOnLineSegment
#
................................................................................
# Examples:
#     - calculateDistanceToPolyline {10 10} {0 0 10 5 20 0}
#       Result: 5.0
#     - calculateDistanceToPolyline {5 10} {0 0 10 5 20 0}
#       Result: 6.7082039325
#
proc ::math::geometry::calculateDistanceToPolyline {P polyline} {
    set minDist "none"
    foreach {Ax Ay} [lrange $polyline 0 end-2] {Bx By} [lrange $polyline 2 end] {
	set dist [calculateDistanceToLineSegment $P [list $Ax $Ay $Bx $By]]
	if {$minDist=="none" || $dist < $minDist} {
	    set minDist $dist
	}

    }
    return $minDist
}

# ::math::geometry::calculateDistanceToPolygon
#
#       Calculate the distance between a point and a polygon.
................................................................................
# Examples:
#     - findClosestPointOnPolyline {10 10} {0 0 10 5 20 0}
#       Result: 10 5
#     - findClosestPointOnPolyline {5 10} {0 0 10 5 20 0}
#       Result: 8.0 4.0
#
proc ::math::geometry::findClosestPointOnPolyline {P polyline} {
    set closestPoint "none"
    foreach {Ax Ay} [lrange $polyline 0 end-2] {Bx By} [lrange $polyline 2 end] {
	set Q [findClosestPointOnLineSegment $P [list $Ax $Ay $Bx $By]]
	set dist [lengthOfPolyline [concat $P $Q]]
	if {$closestPoint=="none" || $dist<$closestDistance} {
	    set closestPoint $Q
	    set closestDistance $dist
	}

    }
    return $closestPoint
}




................................................................................
#
# Examples:
#     - lengthOfPolyline {0 0 5 0 5 10}
#       Result: 15.0
#
proc ::math::geometry::lengthOfPolyline {polyline} {
    set length 0
    foreach {x1 y1} [lrange $polyline 0 end-2] {x2 y2} [lrange $polyline 2 end] {
	set length [expr {$length + sqrt(pow($x1-$x2,2) + pow($y1-$y2,2))}]
	#set length [expr {$length + sqrt(($x1-$x2)*($x1-$x2) + ($y1-$y2)*($y1-$y2))}]
    }
    return $length
}




................................................................................
	foreach part2 $polyline2parts {
	    set part1bbox [bbox $part1]
	    set part2bbox [bbox $part2]
	    if {[rectanglesOverlap [lrange $part1bbox 0 1] [lrange $part1bbox 2 3] \
		    [lrange $part2bbox 0 1] [lrange $part2bbox 2 3] 0]} {
		# the lines' bounding boxes intersect
		if {$perfectmatch} {
		    foreach {l1x1 l1y1} [lrange $part1 0 end-2] {l1x2 l1y2} [lrange $part1 2 end] {
			foreach {l2x1 l2y1} [lrange $part2 0 end-2] {l2x2 l2y2} [lrange $part2 2 end] {
			    if {[lineSegmentsIntersect [list $l1x1 $l1y1 $l1x2 $l1y2] \
				    [list $l2x1 $l2y1 $l2x2 $l2y2]]} {
				# two line segments overlap
				return 1
			    }

			}

		    }
		    return 0
		} else {
		    return 1
		}
	    }
	}
................................................................................
#
# Results:
#       closedpolygon a polygon whose first and last vertices
#                     coincide
#
proc ::math::geometry::ClosedPolygon {polygon} {


    if { [lindex $polygon 0] != [lindex $polygon end-1] ||
         [lindex $polygon 1] != [lindex $polygon end]     } {

        return [concat $polygon [lrange $polygon 0 1]]


    } else {

        return $polygon
    }
}


# ::math::geometry::pointInsidePolygon
#
#       Determine if a point is completely inside a polygon. If the point
#       touches the polygon, then the point is not complete inside the
................................................................................
#       Result: 0
#
proc ::math::geometry::pointInsidePolygon {P polygon} {
    # check if P is on one of the polygon's sides (if so, P is not
    # inside the polygon)
    set closedPolygon [ClosedPolygon $polygon]

    foreach {x1 y1} [lrange $closedPolygon 0 end-2] {x2 y2} [lrange $closedPolygon 2 end] {
	if {[calculateDistanceToLineSegment $P [list $x1 $y1 $x2 $y2]]<0.0000001} {
	    return 0
	}

    }

    # Algorithm
    #
    # Consider a straight line going from P to a point far away from both
    # P and the polygon (in particular outside the polygon).
    #   - If the line intersects with 0 of the polygon's sides, then
................................................................................
        [expr {[lindex $polygonBbox 1]-0.1*[lindex $polygonBbox 3]}]]

    set infinityLine [concat $pointFarAway $P]

    # calculate number of intersections
    set noOfIntersections 0
    #   1. count intersections between the line and the polygon's sides
    foreach {x1 y1} [lrange $closedPolygon 0 end-2] {x2 y2} [lrange $closedPolygon 2 end] {
	if {[lineSegmentsIntersect $infinityLine [list $x1 $y1 $x2 $y2]]} {
	    incr noOfIntersections
	}

    }
    #   2. count intersections between the line and the polygon's points
    foreach {x1 y1} $closedPolygon {
	if {[calculateDistanceToLineSegment [list $x1 $y1] $infinityLine]<0.0000001} {
	    incr noOfIntersections
	}
    }
................................................................................
#
# Examples:
#     - areaPolygon {-10 -10 10 -10 10 10 -10 10}
#       Result: 400
#
proc ::math::geometry::areaPolygon {polygon} {

    foreach {a1 a2 b1 b2} $polygon {break}


    set area 0.0
    foreach {c1 c2} [lrange $polygon 4 end] {
        set area [expr {$area + $b1*$c2 - $b2*$c1}]
        set b1   $c1
        set b2   $c2
    }
    expr {0.5*abs($area)}
}

# ::math::geometry::inproduct
................................................................................
#       vector2       second vector
#
# Results:
#       area          the area of the parallellogram
#
proc ::math::geometry::areaParallellogram {vector1 vector2} {

    foreach {x1 y1} $vector1 {break}
    foreach {x2 y2} $vector2 {break}

    set area [expr {abs($x2 * $y1 - $x1 * $y2}]

    return $area
}

# ::math::geometry::translate
................................................................................
# Results:
#       newPolyline   Translated poyline
#
proc ::math::geometry::translate {vector polyline} {

    set newPolyline $polyline

    foreach {xt yt} $vector {break}

    set idx 0
    foreach {x y} $polyline {
        lset newPolyline $idx [expr {$x + $xt}]
        incr idx
        lset newPolyline $idx [expr {$y + $yt}]
        incr idx
................................................................................
	calculateDistanceToLineSegment findClosestPointOnLineSegment \
	calculateDistanceToPolylineSegment findClosestPointOnPolyline lengthOfPolyline \
	movePointInDirection lineSegmentsIntersect findLineSegmentIntersection findLineIntersection \
	polylinesIntersect polylinesBoundingIntersect intervalsOverlap rectanglesOverlap pointInsidePolygon pointInsidePolygonAlt \
	rectangleInsidePolygon areaPolygon translate rotate reflect degToRad radToDeg
}

package provide math::geometry 1.2.2






>







 







|
<





|
<





|
<





|





|







 







|







 







|





|
<




|




|




<
|



<
|







 







|







 







|







 







|
>
|







 







|
|







 







|







 







|







 







|
|

|


>







 







|
|


|



>







 







|
|
|







 







|
|





>

>







 







>
|
|

<
>

<
|
|
<







 







|



>







 







|



>







 







|
>


|
|







 







|
<







 







|







 







|
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
..
57
58
59
60
61
62
63
64

65
66
67
68
69
70

71
72
73
74
75
76

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
...
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
...
144
145
146
147
148
149
150
151
152
153
154
155
156
157

158
159
160
161
162
163
164
165
166
167
168
169
170
171

172
173
174
175

176
177
178
179
180
181
182
183
...
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
...
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
...
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
...
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
...
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
...
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
...
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
...
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
...
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
...
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
....
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069

1070
1071

1072
1073

1074
1075
1076
1077
1078
1079
1080
....
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
....
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
....
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
....
1391
1392
1393
1394
1395
1396
1397
1398

1399
1400
1401
1402
1403
1404
1405
....
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
....
1551
1552
1553
1554
1555
1556
1557
1558
# See the file "license.terms" for information on usage and redistribution
# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
# RCS: @(#) $Id: geometry.tcl,v 1.12 2010/05/24 21:44:16 andreas_kupries Exp $

namespace eval ::math::geometry {}

package require Tcl 8.5
package require math

###
#
# POINTS
#
#    A point P consists of an x-coordinate, Px, and a y-coordinate, Py,
................................................................................
# Point constructor
proc ::math::geometry::p {x y} {
    return [list $x $y]
}

# Vector addition
proc ::math::geometry::+ {pa pb} {
    lassign $pa ax ay; lassign $pb bx by

    return [list [expr {$ax + $bx}] [expr {$ay + $by}]]
}

# Vector difference
proc ::math::geometry::- {pa pb} {
    lassign $pa ax ay; lassign $pb bx by

    return [list [expr {$ax - $bx}] [expr {$ay - $by}]]
}

# Distance between 2 points
proc ::math::geometry::distance {pa pb} {
    lassign $pa ax ay; lassign $pb bx by

    return [expr {hypot($bx-$ax,$by-$ay)}]
}

# Length of a vector
proc ::math::geometry::length {v} {
    lassign $v x y
    return [expr {hypot($x,$y)}]
}

# Scaling a vector by a factor
proc ::math::geometry::s* {factor p} {
    lassign $p x y
    return [list [expr {$x * $factor}] [expr {$y * $factor}]]
}

# Unit vector into specific direction given by angle (degrees)
proc ::math::geometry::direction {angle} {
    variable torad
    set x [expr {  cos($angle * $torad)}]
................................................................................
proc ::math::geometry::between {pa pb s} {
    return [+ $pa [s* $s [- $pb $pa]]]
}

# Find direction octant the point (vector) lies in.
proc ::math::geometry::octant {p} {
    variable todeg
    lassign $p x y

    set a [expr {(atan2(-$y,$x)*$todeg)}]
    while {$a >  360} {set a [expr {$a - 360}]}
    while {$a < -360} {set a [expr {$a + 360}]}
    if {$a < 0} {set a [expr {360 + $a}]}

    #puts "p ($x, $y) @ angle $a | [expr {atan2($y,$x)}] | [expr {atan2($y,$x)*$todeg}]"
................................................................................
	}
	return east ; # a <= 360.0
    }
}

# Return the NW and SE corners of the rectangle.
proc ::math::geometry::nwse {rect} {
    lassign $rect xnw ynw xse yse
    return [list [p $xnw $ynw] [p $xse $yse]]
}

# Construct rectangle from NW and SE corners.
proc ::math::geometry::rect {pa pb} {
    lassign $pa ax ay; lassign $pb bx by

    return [list $ax $ay $bx $by]
}

proc ::math::geometry::conjx {p} {
    lassign $p x y
    return [list [expr {- $x}] $y]
}

proc ::math::geometry::conjy {p} {
    lassign $p x y
    return [list $x [expr {- $y}]]
}

proc ::math::geometry::x {p} {

    return [lindex $p 0]
}

proc ::math::geometry::y {p} {

    return [lindex $p 1]
}

# ::math::geometry::calculateDistanceToLine
#
#       Calculate the distance between a point and a line.
#
# Arguments:
................................................................................
#     - calculateDistanceToLine {5 10} {0 0 10 10}
#       Result: 3.53553390593
#     - calculateDistanceToLine {-10 0} {0 0 10 10}
#       Result: 7.07106781187
#
proc ::math::geometry::calculateDistanceToLine {P line} {
    # solution based on FAQ 1.02 on comp.graphics.algorithms
    # L = hypot( Bx-Ax, By-Ay )
    #     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
    # s = -----------------------------
    #                 L^2
    # dist = |s|*L
    #
    # =>
    #
................................................................................
    set Bx [lindex $line 2]
    set By [lindex $line 3]
    set Cx [lindex $P 0]
    set Cy [lindex $P 1]
    if {$Ax==$Bx && $Ay==$By} {
	return [lengthOfPolyline [concat $P [lrange $line 0 1]]]
    } else {
	set L [expr {hypot($Bx-$Ax,$By-$Ay)}]
	return [expr {abs(($Ay-$Cy)*($Bx-$Ax)-($Ax-$Cx)*($By-$Ay)) / $L}]
    }
}

# ::math::geometry::findClosestPointOnLine
#
#       Return the point on a line which is closest to a given point.
................................................................................
#                        r=0      P = A
#                        r=1      P = B
#                        r<0      P is on the backward extension of AB
#                        r>1      P is on the forward extension of AB
#                        0<r<1    P is interior to AB
#
proc ::math::geometry::findClosestPointOnLineImpl {P line} {
    # solution based on FAQ 1.02 on comp.graphics.algorithms - but avoid the
    # chain of pow( sqrt(...) ,2) for better precision (& performance).
    #   L^2 = (Bx-Ax)^2 + (By-Ay)^2
    #        (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
    #   r = -------------------------------
    #                     L^2
    #   Px = Ax + r(Bx-Ax)
    #   Py = Ay + r(By-Ay)
    set Ax [lindex $line 0]
    set Ay [lindex $line 1]
................................................................................
    set Bx [lindex $line 2]
    set By [lindex $line 3]
    set Cx [lindex $P 0]
    set Cy [lindex $P 1]
    if {$Ax==$Bx && $Ay==$By} {
	return [list [list $Ax $Ay] 0]
    } else {
	set Lsquared [expr {pow($Bx-$Ax,2) + pow($By-$Ay,2)}]
	set r [expr {(($Cx-$Ax)*($Bx-$Ax) + ($Cy-$Ay)*($By-$Ay))/$Lsquared}]
	set Px [expr {$Ax + $r*($Bx-$Ax)}]
	set Py [expr {$Ay + $r*($By-$Ay)}]
	return [list [list $Px $Py] $r]
    }
}

# ::math::geometry::calculateDistanceToLineSegment
................................................................................
#                        r=1      P = B
#                        r<0      P is on the backward extension of AB
#                        r>1      P is on the forward extension of AB
#                        0<r<1    P is interior to AB
#
proc ::math::geometry::calculateDistanceToLineSegmentImpl {P linesegment} {
    # solution based on FAQ 1.02 on comp.graphics.algorithms
    # L = hypot( Bx-Ax , By-Ay )
    #     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
    # s = -----------------------------
    #                 L^2
    #      (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
    # r = -------------------------------
    #                   L^2
    # dist = |s|*L
................................................................................
    set Bx [lindex $linesegment 2]
    set By [lindex $linesegment 3]
    set Cx [lindex $P 0]
    set Cy [lindex $P 1]
    if {$Ax==$Bx && $Ay==$By} {
	return [list [lengthOfPolyline [concat $P [lrange $linesegment 0 1]]] 0]
    } else {
	set L [expr {hypot($Bx-$Ax,$By-$Ay)}]
	set r [expr {(($Cx-$Ax)*($Bx-$Ax) + ($Cy-$Ay)*($By-$Ay))/pow($L,2)}]
	return [list [expr {abs(($Ay-$Cy)*($Bx-$Ax)-($Ax-$Cx)*($By-$Ay)) / $L}] $r]
    }
}

# ::math::geometry::findClosestPointOnLineSegment
#
................................................................................
# Examples:
#     - calculateDistanceToPolyline {10 10} {0 0 10 5 20 0}
#       Result: 5.0
#     - calculateDistanceToPolyline {5 10} {0 0 10 5 20 0}
#       Result: 6.7082039325
#
proc ::math::geometry::calculateDistanceToPolyline {P polyline} {
    set minDist "Inf"
    foreach {Bx By} [lassign $polyline Ax Ay] {
	set dist [calculateDistanceToLineSegment $P [list $Ax $Ay $Bx $By]]
	if {$dist < $minDist} {
	    set minDist $dist
	}
	set Ax $Bx; set Ay $By
    }
    return $minDist
}

# ::math::geometry::calculateDistanceToPolygon
#
#       Calculate the distance between a point and a polygon.
................................................................................
# Examples:
#     - findClosestPointOnPolyline {10 10} {0 0 10 5 20 0}
#       Result: 10 5
#     - findClosestPointOnPolyline {5 10} {0 0 10 5 20 0}
#       Result: 8.0 4.0
#
proc ::math::geometry::findClosestPointOnPolyline {P polyline} {
    set closestPoint "none"; set closestDistance "Inf"
    foreach {Bx By} [lassign $polyline Ax Ay] {
	set Q [findClosestPointOnLineSegment $P [list $Ax $Ay $Bx $By]]
	set dist [lengthOfPolyline [concat $P $Q]]
	if {$dist<$closestDistance} {
	    set closestPoint $Q
	    set closestDistance $dist
	}
	set Ax $Bx; set Ay $By
    }
    return $closestPoint
}




................................................................................
#
# Examples:
#     - lengthOfPolyline {0 0 5 0 5 10}
#       Result: 15.0
#
proc ::math::geometry::lengthOfPolyline {polyline} {
    set length 0
    foreach {x2 y2} [lassign $polyline x1 y1] {
	set length [expr {$length + hypot($x1-$x2,$y1-$y2)}]
	set x1 $x2; set y1 $y2
    }
    return $length
}




................................................................................
	foreach part2 $polyline2parts {
	    set part1bbox [bbox $part1]
	    set part2bbox [bbox $part2]
	    if {[rectanglesOverlap [lrange $part1bbox 0 1] [lrange $part1bbox 2 3] \
		    [lrange $part2bbox 0 1] [lrange $part2bbox 2 3] 0]} {
		# the lines' bounding boxes intersect
		if {$perfectmatch} {
		    foreach {l1x2 l1y2} [lassign $part1 l1x1 l1y1] {
			foreach {l2x2 l2y2} [lassign $part2 l2x1 l2y1] {
			    if {[lineSegmentsIntersect [list $l1x1 $l1y1 $l1x2 $l1y2] \
				    [list $l2x1 $l2y1 $l2x2 $l2y2]]} {
				# two line segments overlap
				return 1
			    }
			    set l2x1 $l2x2; set l2y1 $l2y2
			}
			set l1x1 $l1x2; set l1y1 $l1y2
		    }
		    return 0
		} else {
		    return 1
		}
	    }
	}
................................................................................
#
# Results:
#       closedpolygon a polygon whose first and last vertices
#                     coincide
#
proc ::math::geometry::ClosedPolygon {polygon} {

    lassign $polygon x y
    if { $x != [lindex $polygon end-1] ||
         $y != [lindex $polygon end]     } {


        lappend polygon $x $y


    }
    return $polygon

}


# ::math::geometry::pointInsidePolygon
#
#       Determine if a point is completely inside a polygon. If the point
#       touches the polygon, then the point is not complete inside the
................................................................................
#       Result: 0
#
proc ::math::geometry::pointInsidePolygon {P polygon} {
    # check if P is on one of the polygon's sides (if so, P is not
    # inside the polygon)
    set closedPolygon [ClosedPolygon $polygon]

    foreach {x2 y2} [lassign $closedPolygon x1 y1] {
	if {[calculateDistanceToLineSegment $P [list $x1 $y1 $x2 $y2]]<0.0000001} {
	    return 0
	}
	set x1 $x2; set y1 $y2
    }

    # Algorithm
    #
    # Consider a straight line going from P to a point far away from both
    # P and the polygon (in particular outside the polygon).
    #   - If the line intersects with 0 of the polygon's sides, then
................................................................................
        [expr {[lindex $polygonBbox 1]-0.1*[lindex $polygonBbox 3]}]]

    set infinityLine [concat $pointFarAway $P]

    # calculate number of intersections
    set noOfIntersections 0
    #   1. count intersections between the line and the polygon's sides
    foreach {x2 y2} [lassign $closedPolygon x1 y1] {
	if {[lineSegmentsIntersect $infinityLine [list $x1 $y1 $x2 $y2]]} {
	    incr noOfIntersections
	}
	set x1 $x2; set y1 $y2
    }
    #   2. count intersections between the line and the polygon's points
    foreach {x1 y1} $closedPolygon {
	if {[calculateDistanceToLineSegment [list $x1 $y1] $infinityLine]<0.0000001} {
	    incr noOfIntersections
	}
    }
................................................................................
#
# Examples:
#     - areaPolygon {-10 -10 10 -10 10 10 -10 10}
#       Result: 400
#
proc ::math::geometry::areaPolygon {polygon} {

    # get last pair of the polygon for start:
    set b1 [lindex $polygon end-1]; set b2 [lindex $polygon end]

    set area 0.0
    foreach {c1 c2} $polygon {
        set area [expr {$area + ($b1*$c2 - $b2*$c1)}]
        set b1   $c1
        set b2   $c2
    }
    expr {0.5*abs($area)}
}

# ::math::geometry::inproduct
................................................................................
#       vector2       second vector
#
# Results:
#       area          the area of the parallellogram
#
proc ::math::geometry::areaParallellogram {vector1 vector2} {

    lassign $vector1 x1 y1; lassign $vector2 x2 y2


    set area [expr {abs($x2 * $y1 - $x1 * $y2}]

    return $area
}

# ::math::geometry::translate
................................................................................
# Results:
#       newPolyline   Translated poyline
#
proc ::math::geometry::translate {vector polyline} {

    set newPolyline $polyline

    lassign $vector xt yt

    set idx 0
    foreach {x y} $polyline {
        lset newPolyline $idx [expr {$x + $xt}]
        incr idx
        lset newPolyline $idx [expr {$y + $yt}]
        incr idx
................................................................................
	calculateDistanceToLineSegment findClosestPointOnLineSegment \
	calculateDistanceToPolylineSegment findClosestPointOnPolyline lengthOfPolyline \
	movePointInDirection lineSegmentsIntersect findLineSegmentIntersection findLineIntersection \
	polylinesIntersect polylinesBoundingIntersect intervalsOverlap rectanglesOverlap pointInsidePolygon pointInsidePolygonAlt \
	rectangleInsidePolygon areaPolygon translate rotate reflect degToRad radToDeg
}

package provide math::geometry 1.2.3

Changes to modules/math/geometry.test.

612
613
614
615
616
617
618







619
620
621
    ::math::geometry::radToDeg [expr {acos(-1.0)}]
} 180.0

test geometry-19.6 {geometry::degToRad} {
    ::math::geometry::degToRad 180.0
} [expr {acos(-1.0)}]









###
testsuiteCleanup






>
>
>
>
>
>
>



612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
    ::math::geometry::radToDeg [expr {acos(-1.0)}]
} 180.0

test geometry-19.6 {geometry::degToRad} {
    ::math::geometry::degToRad 180.0
} [expr {acos(-1.0)}]

##
# areaPolygon
##

test geometry-20.0 {geometry::areaPolygon} {
    math::geometry::areaPolygon {-10 -10 10 -10 10 10 -10 10}
} 400.0

###
testsuiteCleanup

Changes to modules/math/math_geometry.man.

10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
[copyright {2004 by Arjen Markus}]
[copyright {2010 by Andreas Kupries}]
[copyright {2010 by Kevin Kenny}]
[moddesc   {Tcl Math Library}]
[titledesc {Geometrical computations}]
[category  Mathematics]
[require Tcl [opt 8.5]]
[require math::geometry [opt 1.2.2]]

[description]
[para]
The [package math::geometry] package is a collection of functions for
computations and manipulations on two-dimensional geometrical objects,
such as points, lines and polygons.







|







10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
[copyright {2004 by Arjen Markus}]
[copyright {2010 by Andreas Kupries}]
[copyright {2010 by Kevin Kenny}]
[moddesc   {Tcl Math Library}]
[titledesc {Geometrical computations}]
[category  Mathematics]
[require Tcl [opt 8.5]]
[require math::geometry [opt 1.2.3]]

[description]
[para]
The [package math::geometry] package is a collection of functions for
computations and manipulations on two-dimensional geometrical objects,
such as points, lines and polygons.

Changes to modules/math/pkgIndex.tcl.

23
24
25
26
27
28
29
30
31
32
package ifneeded math::machineparameters 0.1   [list source [file join $dir machineparameters.tcl]]

if {![package vsatisfies [package provide Tcl] 8.5]} {return}
package ifneeded math::calculus::symdiff 1.0.1 [list source [file join $dir symdiff.tcl]]
package ifneeded math::bigfloat          2.0.2 [list source [file join $dir bigfloat2.tcl]]
package ifneeded math::numtheory         1.1   [list source [file join $dir numtheory.tcl]]
package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]
package ifneeded math::geometry          1.2.2 [list source [file join $dir geometry.tcl]]
if {![package vsatisfies [package require Tcl] 8.6]} {return}
package ifneeded math::exact             1.0   [list source [file join $dir exact.tcl]]






|


23
24
25
26
27
28
29
30
31
32
package ifneeded math::machineparameters 0.1   [list source [file join $dir machineparameters.tcl]]

if {![package vsatisfies [package provide Tcl] 8.5]} {return}
package ifneeded math::calculus::symdiff 1.0.1 [list source [file join $dir symdiff.tcl]]
package ifneeded math::bigfloat          2.0.2 [list source [file join $dir bigfloat2.tcl]]
package ifneeded math::numtheory         1.1   [list source [file join $dir numtheory.tcl]]
package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]
package ifneeded math::geometry          1.2.3 [list source [file join $dir geometry.tcl]]
if {![package vsatisfies [package require Tcl] 8.6]} {return}
package ifneeded math::exact             1.0   [list source [file join $dir exact.tcl]]

Changes to modules/math/statistics.man.

308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
[call [cmd ::math::statistics::test-Dunnett] [arg alpha] [arg control] [arg args]]
Determine if one or more groups with normally distributed data have the same means as
the group of control data, using Dunnett's test. It is complementary to the ANOVA test.
The procedure returns a list of the comparison results for each group with the control group. Each
element of this list contains: whether the means are likely to be different (1) or not (0)
and the confidence interval the conclusion is based on. The groups may also be stored in a
nested list, just as with the ANOVA test.
[nl]
Note: some care is required if there is only one group to compare the control with:
[example {
    test-Dunnett-F 0.05 $control [list $A]
}]
Otherwise the group A is split up into groups of one element - this is due to an ambiguity.

[list_begin arguments]






|







308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
[call [cmd ::math::statistics::test-Dunnett] [arg alpha] [arg control] [arg args]]
Determine if one or more groups with normally distributed data have the same means as
the group of control data, using Dunnett's test. It is complementary to the ANOVA test.
The procedure returns a list of the comparison results for each group with the control group. Each
element of this list contains: whether the means are likely to be different (1) or not (0)
and the confidence interval the conclusion is based on. The groups may also be stored in a
nested list, just as with the ANOVA test.
[para]
Note: some care is required if there is only one group to compare the control with:
[example {
    test-Dunnett-F 0.05 $control [list $A]
}]
Otherwise the group A is split up into groups of one element - this is due to an ambiguity.

[list_begin arguments]