### Posted to tcl by kbk at Wed Mar 10 23:00:23 GMT 2010view pretty

```# Find intersection of an arbitrary polygon with a convex one.

# cw --
#
#	Does the path (x0,y0)->(x1,y1)->(x2,y2) turn clockwise
#	or counterclockwise?
#
# Parameters:
#	x0, y0 - First point
#	x1, y1 - Second point
#	x2, y2 - Third point
#
# Results:
#	Returns 1 if point 2 is to the right of a line joining points 0
#	and 1, and -1 if it is to the left of the line.
#
#	If the three points are collinear, returns 1 if the third point is
#	the middle point, -1 if point 0 is the middle point, 0 if
#	point 1 is the middle point

proc cw {x0 y0 x1 y1 x2 y2} {
set dx1 [expr {\$x1 - \$x0}]; set dy1 [expr {\$y1 - \$y0}]
set dx2 [expr {\$x2 - \$x0}]; set dy2 [expr {\$y2 - \$y0}]
# (0,0,\$dx1*\$dy2 - \$dx2*\$dy1) is the crossproduct of
# (\$x1-\$x0,\$y1-\$y0,0) and (\$x2-\$x0,\$y2-\$y0,0).
# Its z-component is positive if the turn
# is clockwise, negative if the turn is counterclockwise.
set pr1 [expr {\$dx1 * \$dy2}]
set pr2 [expr {\$dx2 * \$dy1}]
if {\$pr1 > \$pr2} {
# Clockwise
return 1
} elseif {\$pr1 < \$pr2} {
# Counter-clockwise
return -1
} elseif {\$dx1*\$dx2 < 0 || \$dy1*\$dy2 < 0} {
# point 0 is the middle point
return 0
} elseif {(\$dx1*\$dx1 + \$dy1*\$dy1) < (\$dx2*\$dx2 + \$dy2+\$dy2)} {
# point 1 is the middle point
return 0
} else {
# point 2 lies on the segment joining points 0 and 1
return 1
}
}

# intersect --
#
#	Calculate the point of intersection of two lines
#	containing the line segments (x1,y1)-(x2,y2) and (x3,y3)-(x4,y4)
#
# Parameters:
#	x1,y1 x2,y2 - Endpoints of the first line segment
#	x3,y3 x4,y4 - Endpoints of the second line segment
#
# Results:
#	Returns a two-element list containing the point of intersection.
#	Returns an empty list if the line segments are parallel
#	(including the case where the segments are concurrent).

proc intersect {x1 y1 x2 y2 x3 y3 x4 y4} {
set d [expr {(\$y4 - \$y3) * (\$x2 - \$x1)
- (\$x4 - \$x3) * (\$y2 - \$y1)}]
set na [expr {(\$x4 - \$x3) * (\$y1 - \$y3)
- (\$y4 - \$y3) * (\$x1 - \$x3)}]
if {\$d == 0} {
return {}
}
set r [list \
[expr {\$x1 + \$na * (\$x2 - \$x1) / \$d}] \
[expr {\$y1 + \$na * (\$y2 - \$y1) / \$d}]]
return \$r
}

# pairs --
#
#	Coroutine that yields the elements of a list in pairs
#
# Parameters:
#	list - List to decompose
#
# Immediate result:
#	Returns the name of the coroutine
#
# Further results:
#	Returns two-element ranges from the given list, one at a time.
#	Returns {} at the end of the iteration.

proc pairs {list} {
yield [info coroutine]
foreach {x y} \$list {
yield [list \$x \$y]
}
return {}
}

# clipsegment --
#
#	Clips one segment of a polygon against a line.
#
# Parameters:
#	inside0 - Flag = 1 if sx0,sy0 is to the right of the clipping line
#	cx0,cy1 cx1,cy1 - Two points determining the clipping line
#	sx0,sy0 sx1,sy1 - Two points determining the subject line
#
# Results:
#	Returns 1 if sx1,sy1 is to the right of the clipping line, 0 otherwise
#
# Yields:
#	Yields, in order:
#
#		The intersection point of the segment and the line, if
#		the segment intersects the line.
#
#		The endpoint of the segment, if the segment ends to the
#		right of the line.

proc clipsegment {inside0 cx0 cy0 cx1 cy1 sx0 sy0 sx1 sy1} {

set inside1 [expr {[cw \$cx0 \$cy0 \$cx1 \$cy1 \$sx1 \$sy1] > 0}]
if {\$inside1} {
if {!\$inside0} {
set int [intersect \$cx0 \$cy0 \$cx1 \$cy1 \
\$sx0 \$sy0 \$sx1 \$sy1]
if {[llength \$int] >= 0} {
yield \$int
}
}
yield [list \$sx1 \$sy1]
} else {
if {\$inside0} {
set int [intersect \$cx0 \$cy0 \$cx1 \$cy1 \
\$sx0 \$sy0 \$sx1 \$sy1]
if {[llength \$int] >= 0} {
yield \$int
}
}
}
return \$inside1
}

# clipstep --
#
#	Coroutine to perform one step of Sutherland-Hodgman polygon clipping
#
# Parameters:
#	source - Name of a coroutine that will return the vertices of a
#		 subject polygon to clip, and return {} at the end of the
#		 iteration.
#	cx0,cy0 cx1,cy1 - Endpoints of an edge of the clip polygon
#
# Immediate result:
#	Returns the name of the coroutine
#
# Further results:
#	Returns the vertices of the clipped polygon

proc clipstep {source cx0 cy0 cx1 cy1} {
yield [info coroutine]
set pt0 [{*}\$source]
if {[llength \$pt0] == 0} {
return
}
lassign \$pt0 sx0 sy0
set inside0 [expr {[cw \$cx0 \$cy0 \$cx1 \$cy1 \$sx0 \$sy0] > 0}]
set finished 0
while {!\$finished} {
set thispt [{*}\$source]
if {[llength \$thispt] == 0} {
set thispt \$pt0
set finished 1
}
lassign \$thispt sx1 sy1
set inside0 [clipsegment \$inside0 \
\$cx0 \$cy0 \$cx1 \$cy1 \$sx0 \$sy0 \$sx1 \$sy1]
set sx0 \$sx1
set sy0 \$sy1
}
return {}
}

# clippoly --
#
#	Perform Sutherland-Hodgman polygon clipping
#
# Parameters:
#	cpoly - Coordinates of the clip polygon, listed in clockwise order.
#	spoly - Coordinates of the subject polygon, listed in clockwise order.
#
# Results:
#	Returns the coordinates of the clipped polygon, listed in clockwise
#	order.
#
# The clip polygon must be convex. The subject polygon may be any polygon,
# including degenerate and self-intersecting ones.

proc clippoly {cpoly spoly} {
variable clipindx
set source [coroutine clipper[incr clipindx] \
pairs \$spoly]
set cx0 [lindex \$cpoly end-1]
set cy0 [lindex \$cpoly end]
foreach {cx1 cy1} \$cpoly {
set source [coroutine clipper[incr clipindx] \
clipstep \$source \$cx0 \$cy0 \$cx1 \$cy1]
set cx0 \$cx1; set cy0 \$cy1
}
set result {}
while {[llength [set pt [{*}\$source]]] > 0} {
lappend result {*}\$pt
}
return \$result
}

if {![info exists ::argv0] || [string compare \$::argv0 [info script]]} {
return
}

package require Tk

grid [canvas .c -width 400 -height 400 -background \#ffffff]
proc demonstrate {cpoly spoly} {
.c create polygon {*}\$cpoly -outline \#ff9999 -fill {} \
-width 5
.c create polygon {*}\$spoly -outline \#9999ff -fill {} \
-width 3
.c create polygon {*}[clippoly \$cpoly \$spoly] \
-fill \#99ff99 -outline black -width 1
}

demonstrate {100 100 300 100 300 300 100 300} \
{50 150 200 50 350 150 350 300 250 300 200 250 150 350 100 250 100 200}```