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
 Reset

The 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.)