Sutherland-Hodgson polygon clipping, with coroutines

# 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}