Smallest Enclosing Disc
Keith Vetter 2005-11-21 : I was perusing a book on Computational Geometry and I came across an interesting algorithm: consider a robot arm whose base is affixed to the work floor. The arm has to pick up various items and place them elsewhere. What would be a good position for the base of the arm? How about "somewhere in the middle" of the points it has to reach. If you locate it in the middle of the smallest enclosing disc then you minimize the maximum distance the arm has to travel.
The algorithm used below works by incrementally adding points to an existing solution, noting that if the new point is not in the current smallest disc, then the new disc will have that point on its perimeter. This produces an O(n) solution; there are faster ones out there but this was fun to program.
##+##########################################################################
#
# SmallestDisc.tcl -- Computes the smallest enclosing disc for a set of points
# by Keith Vetter, Nov 21, 2005
#
package require Tk
if {! [catch {package require tile}]} { ;# Be tile friendly
namespace import -force ::ttk::button
namespace import -force ::ttk::label
}
# Random algorithm for computing smallest enclosing disk
proc MiniDisc {PTS0} {
set PTS [Shuffle $PTS0]
set pre0 [lrange $PTS 0 1]
set D [SegmentCenter [lindex $PTS 0] [lindex $PTS 1]]
foreach pt [lrange $PTS 2 end] {
if {! [Within $D $pt]} {
set D [MiniDiscWithPoint $pre0 $pt]
}
lappend pre0 $pt
}
return $D
}
# Returns smallest enclosing disk for points PTS with q on perimeter
proc MiniDiscWithPoint {PTS q} {
set pre1 [lindex $PTS 0]
set D [SegmentCenter $pre1 $q]
foreach pt [lrange $PTS 1 end] {
if {! [Within $D $pt]} {
set D [MiniDiscWith2Points $pre1 $pt $q]
}
lappend pre1 $pt
}
return $D
}
# Returns smallest enclosing disk for points PTS with q1, q2 on perimeter
proc MiniDiscWith2Points {PTS q1 q2} {
set D [SegmentCenter $q1 $q2]
foreach pt $PTS {
if {! [Within $D $pt]} {
set D [CircumCenter $pt $q1 $q2]
}
}
return $D
}
# Returns circle that has P1-P2 as a diameter
proc SegmentCenter {P1 P2} {
set h12 [VAdd $P1 [VSub $P2 $P1] .5]
foreach {rx ry} [XY [VSub $P1 $h12]] break ;# Radius vector
set radius [expr {hypot($rx,$ry)}] ;# Radius magnitude
return [concat $h12 $radius]
}
# compute circumcenter for three points
proc CircumCenter {P1 P2 P3} {
set h12 [VAdd $P1 [VSub $P2 $P1] .5] ;# Midpoints of each side
set h13 [VAdd $P1 [VSub $P3 $P1] .5]
set n12 [VNormal [VSub $P2 $P1]] ;# Normal to side p1-p2
set n13 [VNormal [VSub $P3 $P1]]
set O [IntersectV $h12 $n12 $h13 $n13] ;# The circumcenter
foreach {rx ry} [XY [VSub $P1 $O]] break ;# Radius vector
set radius [expr {hypot($rx,$ry)}] ;# Radius magnitude
return [concat $O $radius]
}
# returns true if pt is inside circle D
proc Within {D pt} {
foreach {ox oy} [XY [lindex $D 0]] break
set r [lindex $D 1]
foreach {x y} [XY $pt] break
return [expr {(($x-$ox)*($x-$ox) + ($y-$oy)*($y-$oy)) <= $r*$r}]
}
# splits apart our point data structure
proc XY {P} {
return [split $P ","]
}
# adds to vectors, possibly scaling the second one
proc VAdd {v1 v2 {scaling 1}} {
foreach {x1 y1} [XY $v1] {x2 y2} [XY $v2] break
return "[expr {$x1 + $scaling*$x2}],[expr {$y1 + $scaling*$y2}]"
}
proc VSub {v1 v2} { return [VAdd $v1 $v2 -1] }
proc VNormal {v} {foreach {x y} [XY $v] break; return "$y,[expr {-1 * $x}]"}
##+##########################################################################
#
# IntersectV -- find where 2 point/vector intersect
#
# p1+K(v1) = p3+J(v3)
# convert into and solve matrix equation (a b / c d) ( K / J) = ( e / f )
#
proc IntersectV {p1 v1 p3 v3} {
foreach {x1 y1} [XY $p1] {vx1 vy1} [XY $v1] {x3 y3} [XY $p3] {vx3 vy3} [XY $v3] break
set a $vx1
set b [expr {-1 * $vx3}]
set c $vy1
set d [expr {-1 * $vy3}]
set e [expr {$x3 - $x1}]
set f [expr {$y3 - $y1}]
set det [expr {double($a*$d - $b*$c)}]
if {$det == 0} {error "Determinant is 0"}
set k [expr {($d*$e - $b*$f) / $det}]
#set j [expr {($a*$f - $c*$e) / $det}]
return [VAdd $p1 $v1 $k]
}
proc Shuffle { l } {
set len [llength $l]
set len2 $len
for {set i 0} {$i < $len-1} {incr i} {
set n [expr {int($i + $len2 * rand())}]
incr len2 -1
set temp [lindex $l $i] ;# Swap elements at i & n
lset l $i [lindex $l $n]
lset l $n $temp
}
return $l
}
################################################################
#
# GUI stuff below
#
proc DoDisplay {} {
canvas .c -bg white -width 500 -height 500 -highlightthickness 0 -bd 2 -relief ridge
bind .c <1> [list Click %x %y]
label .l -justify l -relief ridge \
-text "Click to create a new point\nRight click to delete point"
button .reset -text "Reset" -command Reset -state disabled
button .about -text "About" -command About
pack .c -fill both -expand 1
pack .l -side left -expand 1 -fill x
pack .reset .about -side left -expand 1
bind all <Key-F2> [list console show]
}
proc Reset {} {
array unset ::P
set ::P(id) 0
.c delete all
.c create text 250 10 -anchor n -font {Times 18 bold underline} \
-text "Smallest Enclosing Disk Demo"
ShowDisc
}
proc Click {x y} {
global P
set x [.c canvasx $x]
set y [.c canvasy $y]
set P(pt,[incr P(id)]) "$x,$y"
set xy [Expand $x $y 5]
.c create oval $xy -tag pt$P(id) -fill magenta -outline black
.c bind pt$P(id) <3> [list DeletePoint $P(id)]
ShowDisc
}
proc DeletePoint {id} {
global P
unset -nocomplain P(pt,$id)
.c delete pt$id
ShowDisc
}
proc Expand {x y dxy} {
set x0 [expr {$x-$dxy}]
set y0 [expr {$y-$dxy}]
set x1 [expr {$x+$dxy}]
set y1 [expr {$y+$dxy}]
return [list $x0 $y0 $x1 $y1]
}
proc ShowDisc {} {
global P
.c delete disc
set PTS {}
foreach a [array name P pt,*] {
lappend PTS $P($a)
}
.reset config -state [expr {$PTS eq {} ? "disabled" : "normal"}]
if {[llength $PTS] < 2} return
set D [MiniDisc $PTS]
foreach {x y} [XY [lindex $D 0]] break
set r [lindex $D 1]
.c delete disc
set xy [Expand $x $y $r]
.c create oval $xy -tag {x disc} -fill turquoise1 -outline turquoise1
set xy [Expand $x $y 3]
.c create oval $xy -tag disc -fill red -outline black
.c create line $x $y [expr {$x+$r}] $y -tag disc -fill red -arrow last
.c lower disc
}
proc About {} {
set txt "Smallest Enclosing Disc\nby Keith Vetter, Nov 2005\n\n"
append txt "Consider a robot arm attached to the work floor. The arm\n"
append txt "has to pick up various items scattered around it and place\n"
append txt "them at other points. What would be the best position for\n"
append txt "the robot arm? Somewhere 'in the middle' of the points it\n"
append txt "has to reach. How about at the center of the smallest disc\n"
append txt "that encloses all the points.\n\n"
append txt "Formally, given a set P of n points in the plane (the points\n"
append txt "the robot arm must be able to reach), find the SMALLEST\n"
append txt "ENCLOSING DISC for P, that is, the smallest disc that\n"
append txt "contains all the points in P\n"
tk_messageBox -message $txt
}
DoDisplay
ResetThe above code needs a "set P(id) 0" in there somewhere...
KPV It does it in Reset which gets called by the last line of code. I bet you cut and pasted the code into a console and the last line didn't get executed.
uniquename 2013jul29
This code could use an image to show what it produces. (I added a frame for the label and 2 buttons and packed them at the top of the GUI, instead of them appearing at the bottom of the GUI.)
