Version 15 of Convex hull

Updated 2004-09-20 00:11:10 by kbk

http://mini.net/files/convexhull.jpg

From a posting in the comp.lang.tcl newsgroup:

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

  1. (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.
  2. 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.
  3. 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 } {

     foreach { x0 y0 } $p0 break
     foreach { x1 y1 } $p1 break
     foreach { x2 y2 } $p2 break

     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]
     foreach { x0 y0 } $p0 break
     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 {
         foreach { x y } $p break
         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 {
         foreach { x y } $p break
         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] {
         foreach { - x y } $p break
         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 illustrate { command } {
     uplevel 1 $command
 }

 # Plot the raw data

 proc plotPoints { points } {
     set n 0
     foreach p $points {
         foreach { x y } $p break
         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 {
         foreach { x y } $p break
         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
 }

 expr { srand( 4 ) }
 grid [canvas .c -width 586 -height 512 -background white]
 for { set i 0 } { $i < 32 } { incr i } {
     lappend points [list [expr { rand() }] [expr { rand() }]]
 }
 hull2d $points

Category Mathematics