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 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
    #
    # =>
    #
    #        | (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))}]
	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

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

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

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

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