Tcl Library Source Code

# 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} {
    foreach {ax ay} $pa break
    lassign $pa ax ay; lassign $pb bx by
    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
    lassign $pa ax ay; lassign $pb bx by
    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
    lassign $pa ax ay; lassign $pb bx by
    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
    lassign $v x y
    return [expr {hypot($x,$y)}]
}

# Scaling a vector by a factor
proc ::math::geometry::s* {factor p} {
    foreach {x y} $p break
    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
    foreach {x y} $p break
    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} {
    foreach {xnw ynw xse yse} $rect break
    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} {
    foreach {ax ay} $pa break
    lassign $pa ax ay; lassign $pb bx by
    foreach {bx by} $pb break
    return [list $ax $ay $bx $by]
}

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

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

proc ::math::geometry::x {p} {
    foreach {x y} $p break
    return $x
    return [lindex $p 0]
}

proc ::math::geometry::y {p} {
    foreach {x y} $p break
    return $y
    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 = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
    # L = hypot( Bx-Ax, By-Ay )
    #     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
    # s = -----------------------------
    #                 L^2
    # dist = |s|*L
    #
    # =>
    #
    #        | (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) |
    # dist = ---------------------------------
    #                       L
    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 [lengthOfPolyline [concat $P [lrange $line 0 1]]]
    } else {
	set L [expr {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}]
	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
    #   L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
    # 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 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 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 = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
    # 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 {sqrt(pow($Bx-$Ax,2) + pow($By-$Ay,2))}]
	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 "none"
    foreach {Ax Ay} [lrange $polyline 0 end-2] {Bx By} [lrange $polyline 2 end] {
    set minDist "Inf"
    foreach {Bx By} [lassign $polyline Ax Ay] {
	set dist [calculateDistanceToLineSegment $P [list $Ax $Ay $Bx $By]]
	if {$minDist=="none" || $dist < $minDist} {
	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"
    foreach {Ax Ay} [lrange $polyline 0 end-2] {Bx By} [lrange $polyline 2 end] {
    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 {$closestPoint=="none" || $dist<$closestDistance} {
	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 {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))}]
    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 {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] {
		    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 { [lindex $polygon 0] != [lindex $polygon end-1] ||
         [lindex $polygon 1] != [lindex $polygon end]     } {
    if { $x != [lindex $polygon end-1] ||
         $y != [lindex $polygon end]     } {

        return [concat $polygon [lrange $polygon 0 1]]
        lappend polygon $x $y

    } else {

        return $polygon
    }
    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] {
    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 {x1 y1} [lrange $closedPolygon 0 end-2] {x2 y2} [lrange $closedPolygon 2 end] {
    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:
    foreach {a1 a2 b1 b2} $polygon {break}
    set b1 [lindex $polygon end-1]; set b2 [lindex $polygon end]

    set area 0.0
    foreach {c1 c2} [lrange $polygon 4 end] {
        set area [expr {$area + $b1*$c2 - $b2*$c1}]
    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} {

    foreach {x1 y1} $vector1 {break}
    lassign $vector1 x1 y1; lassign $vector2 x2 y2
    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}
    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.2
package provide math::geometry 1.2.3

    ::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

[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]]
[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.

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]]
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]]

[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]
[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]