Convex hull
From a posting in comp.lang.tcl:
> 1. DOn't suppose anyone has a convex hull procedure in TCL that they could > share? > > 2. Don't suppose anyone wants three versions of a convex hull procedure that > doesn't work?
Arjen Markus replies:
> I do not have any procedures to share (working or none), but have you > checked a repository like <http://www.netlib.org>? They may have a > working version for you in, say, C or Fortran: > - Steal^H^H^H^H^HExamine the code and re-implement in Tcl > - Make a small extension that directly uses the code > And above all: if you do the first, put it on the Wiki :)
KBK (5 November 2002) - Arjen, your wish is my command!
The convex hull of a set of points in two dimensions is the smallest convex polygon that encloses the points. In spaces of higher dimension, it is the smallest convex polytope that encloses the points.
Finding the convex hull in two dimensions is fairly easy using Graham's algorithm (any good undergrad text on algorithms or computational geometry will give it). Essentially, what one does is:
- (Optional) Pre-scan the points to find a set that can be trivially excluded from the hull. This is usually done by finding the points that give the minimum and maximum values of x+y and x-y, and excluding either the interior of that quadrilateral or (faster and almost as effective) the interior of the largest rectangle with sides parallel to the axes that is enclosed by the quadrilateral.
- Locate the point p with minimum y. Sort the remaining points q by the slope of the line joining p and q so that they are in counterclockwise sequence.
- Walk the counterclockwise sequence of points. When adding a point to the hull, walk backward through the points added earlier and remove any that are now interior points.
Finding the convex hull in spaces of dimension greater than two is more of a challenge. What I've done on the occasions that I've needed such a calculation is simply to write the point list to a file, execute the qconvex program from the qhull library [L1 ] on it, and read in the results. Since qhull is also available as a C-callable function library, it would be possible to wrap it into Tcl directly, but I've never had an application that was performance-critical enough to bother.
Without further ado, here is the code for the Graham scan in two dimensions. Note the calls to the [illustrate] function; in Tcl 8.4 and beyond, the function as supplied generates no bytecode and has zero cost of evaluation. The [illustrate] function can be redefined (see below) to show the progress of the algorithm.
# Procedures for developing the convex hull of a set of points in
# two dimensions.
# Copyright (c) 2002 by Kevin B. Kenny. All rights reserved.
# See the file,
# 'http://cvs.sourceforge.net/cgi-bin/viewcvs.cgi/tcllib/tcllib/license.terms'
# for terms and conditions of redistribution
package require Tcl 8.4; # This code uses [lset]
# illustrate --
#
# Illustrate the behavior of the convex hull algorithm.
#
# The default illustration does nothing.
proc illustrate args {}
# ccw --
#
# Determines the direction of turn when going from a first point
# through a second to a third.
#
# Parameters:
# p0, p1, p2 -- Three points, in order on the path. Each point is
# a two-element list comprising abscissa and ordinate.
# Additional list elements are ignored.
#
# Results:
# Returns 1 if the path p0-p1-p2 proceeds in a counterclockwise
# direction, -1 if it is clockwise, 0 if it is neither.
proc ccw { p0 p1 p2 } {
lassign $p0 x0 y0
lassign $p1 x1 y1
lassign $p2 x2 y2
set dx1 [expr { $x1 - $x0 }]
set dy1 [expr { $y1 - $y0 }]
set dx2 [expr { $x2 - $x0 }]
set dy2 [expr { $y2 - $y0 }]
set prod1 [expr { $dx1 * $dy2 }]
set prod2 [expr { $dy1 * $dx2 }]
if { $prod1 > $prod2 } {
return 1
}
if { $prod1 < $prod2 } {
return -1
}
if { $dx1*$dx2 < 0 || $dy1*$dy2 < 0 } {
return -1
}
if { $dx1*$dx1 + $dy1*$dy1 < $dx2*$dx2 + $dy2*$dy2 } {
return 1
}
return 0
}
# hull2d --
#
# Computes the convex hull of a set of points in two dimensions.
#
# Parameters:
# points - List of points for which the hull is to be determined.
# Each point is a 2-element list comprising abscissa and
# ordinate.
#
# Results:
# Returns the convex hull.
proc hull2d { points } {
# Rule out the trivial cases
if { [llength $points] < 3 } { return $points }
illustrate { plotPoints $points }
# Find the points having the largest and smallest values of x+y,
# largest and smallest values of x-y, and smallest value of y.
# Break ties on "value of y" by choosing the point with the smallest x.
set p0 [lindex $points 0]
lassign $p0 x0 y0
set minxpy [expr { $x0 + $y0 }]; set xminxpy $x0; set yminxpy $y0
set maxxpy [expr { $x0 + $y0 }]; set xmaxxpy $x0; set ymaxxpy $y0
set minxmy [expr { $x0 - $y0 }]; set xminxmy $x0; set yminxmy $y0
set maxxmy [expr { $x0 - $y0 }]; set xmaxxmy $x0; set ymaxxmy $y0
set miny $y0; set iminy 0; set xminy $x0
set n 0
foreach p $points {
lassign $p x y
if { $y < $miny || ( $y == $miny && $x < $xminy ) } {
set miny $y; set iminy $n; set xminy $x
}
set xpy [expr { $x + $y }]
if { $xpy < $minxpy } {
set minxpy $xpy; set xminxpy $x; set yminxpy $y
}
if { $xpy > $maxxpy } {
set maxxpy $xpy; set xmaxxpy $x; set ymaxxpy $y
}
set xmy [expr { $x - $y }]
if { $xmy < $minxmy } {
set minxmy $xmy; set xminxmy $x; set yminxmy $y
}
if { $xmy > $maxxmy } {
set maxxmy $xmy; set xmaxxmy $x; set ymaxxmy $y
}
incr n
}
illustrate { corner $xminxpy $yminxpy min(x+y) }
illustrate { corner $xmaxxpy $ymaxxpy max(x+y) }
illustrate { corner $xminxmy $yminxmy min(x-y) }
illustrate { corner $xmaxxmy $ymaxxmy max(x-y) }
illustrate { corner $xminy $miny min(y) }
# Determine bounds on x and y that define a rectangle in whose interior
# no point can be on the hull.
if { $xmaxxpy > $xmaxxmy } {
set xmaxint $xmaxxmy
} else {
set xmaxint $xmaxxpy
}
if { $ymaxxpy > $yminxmy } {
set ymaxint $yminxmy
} else {
set ymaxint $ymaxxpy
}
if { $xminxpy > $xminxmy } {
set xminint $xminxpy
} else {
set xminint $xminxmy
}
if { $yminxpy > $ymaxxmy } {
set yminint $yminxpy
} else {
set yminint $ymaxxmy
}
illustrate { rectangle $xminint $yminint $xmaxint $ymaxint }
# For all points outside the "interior rectangle", compute the angle of
# the ray joining $xminy,$miny to $x,$y. List the points, excluding
# $xminy, $miny itself, sorted by this angle. (This list is a simple
# closed curve comprising all the points on the hull plus a few interior
# points).
set n 0
foreach p $points {
lassign $p x y
if { $n != $iminy
&& ( $x <= $xminint || $x >= $xmaxint
|| $y <= $yminint || $y >= $ymaxint ) } {
lappend pts2 [linsert $p 0 [expr { atan2( $y-$miny, $x-$xminy ) }]]
}
incr n
}
set pts3 [list [list $xminy $miny]]
foreach p [lsort -real -index 0 $pts2] {
lassign $p - x y
lappend pts3 [list $x $y]
}
lappend pts3 [list $xminy $miny]
# Eliminate interior points from the hull under construction
set i 0
set M 2 ; set Mm1 1
foreach pi $pts3 {
if { $i < 3 } {
incr i ; continue
}
while { [ccw [lindex $pts3 $M] [lindex $pts3 $Mm1] $pi] >= 0 } {
incr M -1
incr Mm1 -1
}
incr M
incr Mm1
set temp [lindex $pts3 $M]
lset pts3 $i $temp
lset pts3 $M $pi
set pu $temp
set partial [lrange $pts3 0 $M]
lappend partial [list $xminy $miny]
illustrate { polyline partial $partial -fill blue }
incr i
}
return [lrange $pts3 1 $M]
}# To see the algorithm in operation, include the following code, and the convex hull will be found for a typical set of 32 points. The progress of the algorithm will be displayed in a Tk canvas.
#----------------------------------------------------------------------
#
# Procedures for illustrating the progress of the algorithm
package require Tk
proc wait {ms} {
if {$ms == 0} return
after $ms set ::_wait_sem 1
vwait ::_wait_sem
}
proc illustrate { command } {
uplevel 1 $command
wait $::delay
}
# Plot the raw data
proc plotPoints { points } {
set n 0
foreach p $points {
lassign $p x y
set cx [expr { 500 * $x + 6 }]
set cy [expr { 506 - 500 * $y }]
.c create oval [expr { $cx - 3 }] [expr { $cy - 3 }] \
[expr { $cx + 3 }] [expr { $cy + 3 }] \
-fill black -tags point
}
return
}
# Plot the corner points used for interior point elimination
proc corner { x y label } {
set cx [expr { 500 * $x + 6 }]
set cy [expr { 506 - 500 * $y }]
.c create rectangle [expr { $cx - 4 }] [expr { $cy - 4 }] \
[expr { $cx+4 }] [expr { $cy + 4 }] \
-fill {} -outline red -tags corner
.c create text [expr { $cx + 5 }] $cy -text $label -anchor w
}
# Plot the rectangle used for interior point elimination
proc rectangle { x0 y0 x1 y1 } {
set cx0 [expr { 500 * $x0 + 6 }]
set cy0 [expr { 506 - 500 * $y0 }]
set cx1 [expr { 500 * $x1 + 6 }]
set cy1 [expr { 506 - 500 * $y1 }]
.c create rectangle $cx0 $cy0 $cx1 $cy1 -fill \#cccccc -outline black \
-tags interior
.c lower interior point
return
}
# Plot a partial hull under construction
proc polyline { tag pts args } {
set cmd [list .c create line]
foreach p $pts {
lassign $p x y
lappend cmd [expr { 500 * $x + 6 }] [expr { 506 - 500 * $y}]
}
lappend cmd -tags $tag
foreach arg $args {
lappend cmd $arg
}
eval $cmd
.c lower $tag point
return
}
proc doHull {} {
.c delete all
for { set i 0 } { $i < 32 } { incr i } {
lappend points [list [expr { rand() }] [expr { rand() }]]
}
hull2d $points
}
set delay 400
expr { srand( 4 ) }
grid [canvas .c -width 586 -height 512 -background white]
doHull
bind . <Button-1> doHull
bind . <q> exitIDG Jan08/10
- What are the license terms? Above url gives 404.
- It's not robust with respect to duplicate points. Removing these from the pts3 list right after its formation seems to work.
Kieran Sorry, it does not seem to be working for me on this set of points ...
hull2d {{516063 176648} {520512 178524} {516852 179878} {518523 181191} {489189 412073}}Lars H: You're probably running into the fact that Tcl 8.4 integers on most platforms are 32 bits wide, which means products of pretty much any pair of numbers from the above list silently overflow, producing complete garbage results. Proper integers for Tcl 8.5 or 9.0 contains some grumbling about related issues. Things that you can do to overcome this problem are:
- Switch to (a sufficiently recent) Tcl 8.5. This may be a bit difficult, as the change in integer handling happened after the most recent official release 8.5a3, but many Batteries Included distributions have Tcl8.5 builds from more recent sources.
- Change the numbers above to reals by appending a ".0" to each of them. This prevents overflows by may lead to some loss of precision for very large numbers. (You should however not experience any such thing as long as your numbers are below 16 millions, I think.)
- Change the numbers above to wide integers:
proc wide_int_list {L} {
set res {}
foreach int $L {lappend res [expr {wide($int)}]}
return $res
}
proc wide_int_list_list {L} {
set res {}
foreach l $L {lappend res [wide_int_list $l]}
return $res
}
hull2d [wide_int_list_list {{516063 176648} {520512 178524} {516852 179878} {518523 181191} {489189 412073}}](That this last thing works in Tcl8.4 is really unholy black magic as far as The Tcl Way is concerned, becuase the thing going into wide_int_list_list has the same string representation as the thing coming out, but they don't behave the same when you start doing arithmetic on them!)
2007-09-24 UKo: added a little bit of sugar to make a real animation of the algorithm. Just change the delay or set it to 0. Mouseclick redraws with new random points and 'q' exits.
AMG: Here is an alternate implementation which avoids floating point, division, and trigonometric functions. It is interactive--- just drag the points around. The numbers show the sort order, the orange polygon is the convex hull, the small red circle indicates the pivot point, the dashed blue lines are angles to the pivot, and the red triangles are the rejected clockwise turns. All points are considered as potential hull vertices; the interior rectangle optimization is not applied.
#!/bin/sh
# The next line restarts with tclsh.\
exec tclsh "$0" ${1+"$@"}
package require Tk
proc compare {pivot points a b} {
lassign [lindex $points $pivot] xp yp
lassign [lindex $points $a] xa ya
lassign [lindex $points $b] xb yb
if {$a == $pivot} {
return -1
} elseif {$b == $pivot} {
return 1
} else {
expr {($xb - $xp) * ($ya - $yp) - ($xa - $xp) * ($yb - $yp)}
}
}
proc convhull {} {
global count
.c delete hull pivot label angle error
set h [winfo height .c]
set pivot 0
for {set i 0} {$i < $count} {incr i} {
set coords [.c coords point$i]
set x [expr {int(([lindex $coords 0] + [lindex $coords 2]) / 2)}]
set y [expr {int(([lindex $coords 1] + [lindex $coords 3]) / 2)}]
set y [expr {$h - $y}]
lappend indexes $i
lappend points [list $x $y]
if {$y < [lindex $points $pivot 1]} {
set pivot $i
}
}
set indexes [lsort -command [list compare $pivot $points] $indexes]
lassign [lindex $points $pivot] xp yp
.c create oval [expr {$xp - 5}] [expr {$h + 5 - $yp}]\
[expr {$xp + 5}] [expr {$h - 5 - $yp}] -outline red -tags pivot
for {set i 0} {$i < $count} {incr i} {
lassign [lindex $points [lindex $indexes $i]] x y
.c create text [expr {$x + 5}] [expr {$h - $y}]\
-anchor w -tags label -text $i
.c create line $xp [expr {$h - $yp}] $x [expr {$h - $y}]\
-fill blue -dash - -tags angle
}
set curIndex 4
while {1} {
lassign [lindex $points [lindex $indexes $curIndex-3]] x0 y0
lassign [lindex $points [lindex $indexes $curIndex-2]] x1 y1
lassign [lindex $points [lindex $indexes $curIndex-1]] x2 y2
if {($x1 - $x0) * ($y2 - $y0) < ($y1 - $y0) * ($x2 - $x0)} {
.c create line $x0 [expr {$h - $y0}] $x1 [expr {$h - $y1}]\
$x2 [expr {$h - $y2}] $x0 [expr {$h - $y0}]\
-fill red -tags error
set indexes [lreplace $indexes $curIndex-2 $curIndex-2]
incr curIndex -1
} elseif {$curIndex == [llength $indexes]} {
break
} else {
incr curIndex
}
}
set coords {}
foreach index $indexes {
lappend coords [lindex $points $index 0]
lappend coords [expr {$h - [lindex $points $index 1]}]
}
lappend coords [lindex $points $pivot 0]
lappend coords [expr {$h - [lindex $points $pivot 1]}]
.c create line {*}$coords -fill orange -tags hull
.c lower hull
.c lower angle
.c lower error
}
proc down {x y} {
set ::click [list $x $y]
}
proc drag {id x y} {
global click
set dx [expr {$x - [lindex $click 0]}]
set dy [expr {$y - [lindex $click 1]}]
set click [list $x $y]
.c move point$id $dx $dy
convhull
}
set count 10
canvas .c -highlightthickness 0 -background white -width 500 -height 500
pack .c -fill both -expand true
for {set i 0} {$i < $count} {incr i} {
set y [expr {int(rand() * 480 + 10)}]
set x [expr {int(rand() * 480 + 10)}]
.c create oval [expr {$x - 3}] [expr {$y - 3}]\
[expr {$x + 3}] [expr {$y + 3}]\
-fill black -tags [list point point$i]
.c bind point$i <B1-Motion> [list drag $i %x %y]
}
.c bind point <ButtonPress-1> [list down %x %y]
bind .c <Configure> {bind .c <Configure> {}; convhull}
# vim: set sts=4 sw=4 tw=80 et ft=tcl:TV Hm, reminds me of my graduation work, where the concept of the convex hull applied to (rational, qubic) Bezier surfaces in 3D (more dimensions would be possible).
In that case, apart from how to compute the exact bounding 3-point polygons which together form the outer edges of a surrounding polyhedron, the control points of such surfaces, can be used to form a polyhedron which is a closer bounding volume than a rectangular bounding box and which can be proven to always enclose the Bezier surface spanned by the control points. Like the scaffolding for the hull of a boat, plus a roof encloses the hull being built (which is then the 3D surface and the control polygons are the scaffolding).
See Bwise deCasteljau algorithm example where the principle is used to approximate the bezier surfaces by subdividing them and using the theoretical rule that the control polygons (and thus in almost all cases the convext hull spanned by them, too) converge to the actual curve. No black magic at all, just 3D vector/surface processing.
tclfreak - 2017-04-12 05:15:58
I know this code is old, but I found some bugs in the code:
proc ccw{} does not handle the scenario where there are several co-linear points with identical y axis values, and these points reside along what would be the bottom-left edge of the convex hull. With my data, I see that the 4th 'if' clause evaluates to true and therefore incorrectly returns a '1' for these points. It should return '0'. This leads to a catastrophic failure for the subsequent set of co-linear points. :
Error: can't read "x1": no such variable
Changing the code to this helps prevent the Error above:
if { !($dy1 && == 0 && $dy2 == 0 } { return 0 }
if { $prod1 > $prod2 } {
return 1
}
...Even with this fix, the code does not grab all the correct points...at least one of those co-linear points will be 'eliminated'. To fix that, the following code
while { [ccw [lindex $pts3 $M] [lindex $pts3 $Mm1] $pi] >= 0 } {
can be changed to:
while { [ccw [lindex $pts3 $M] [lindex $pts3 $Mm1] $pi] > 0 } {
With *my* data, at least, this solves the problem.
arjen - 2017-04-12 07:26:41
The code serves as an illustration mainly, but if you want to use it for serious purposes, maybe we should add the algorithm to the math::geometry module in Tcllib. Of course, this is one of the better known problems in computational geometry, which has its very own interesting numerical issues (it is for instance quite common to see references to exact arithmetics).
As a side note: I intend to incorporate a package called Triangle, by Jonathan Shewchuk, in my mathemaTcl metapackage. It may take a while before that is done though, as progress is slow and the package deals with triangulations, Voronoi diagrams and so on.