[PD] I was just curious whether this would be too slow to be useable. It's actually quite quick.  But the data structure needs to be "tclified".
----
    #!/bin/sh
    # the next line restarts using tclsh \
    exec wish8.5 "$0" "$@"

    if {0} {
    For details, see 
    "Primitives for the Manipulation of General Subdivisions and the Computation of Voronoi Diagrams"
    L. Guibas, J. Stolfi, ACM Transactions on Graphics, April 1985
    http://portal.acm.org/citation.cfm?doid=282918.282923 ($10)
    and
    Jorge Stolfi's C libraries
    http://www.ic.unicamp.br/~stolfi/EXPORT/software/c/Index.html
    (this tcl code is a direct translation, right down to the bit twiddling of pointers)
    also
    http://www.cs.cmu.edu/afs/andrew/scs/cs/15-463/2001/pub/src/a2/lischinski/114.ps
    describes an incremental algorithm using the quadedge data structure in c++
    }

    package require Tcl 8.5

    namespace eval quad {}

    proc quad::reset {} {
     variable next; array unset next 
     variable org;  array unset org 
     variable mark; array unset mark 
     variable nextedge 0
     variable nextmark 1
     variable sites {}
     variable index {}
    }

    proc quad::EDGE e {expr {$e&~3}}
    proc quad::ROT  e {expr {($e&~3)+(($e+1)&3)}}
    proc quad::SYM  e {expr {($e&~3)+(($e+2)&3)}}
    proc quad::TOR  e {expr {($e&~3)+(($e+3)&3)}}

    proc quad::ONEXT e {
      variable next
      return $next($e) 
    }
    proc quad::SETNEXT {e n} {
      variable next
      set next($e) $n
    }

    proc quad::RNEXT e {TOR [ONEXT [ROT $e]]}
    proc quad::DNEXT e {SYM [ONEXT [SYM $e]]}
    proc quad::LNEXT e {ROT [ONEXT [TOR $e]]}
    proc quad::OPREV e {ROT [ONEXT [ROT $e]]}
    proc quad::DPREV e {TOR [ONEXT [TOR $e]]}
    proc quad::RPREV e {ONEXT [SYM $e]}
    proc quad::LPREV e {SYM [ONEXT $e]}

    proc quad::MARK e {
     variable mark
     return $mark($e)
    }
    proc quad::SETMARK {e m} {
     variable mark
     set mark($e) $m
    }

    proc quad::ORG e {
     variable org
     return $org($e)
    }
    proc quad::SETORG {e p} {
      variable org
      set org($e) $p
    }
    proc quad::DEST e {ORG [SYM $e]}
    proc quad::SETDEST {e p} {SETORG [SYM $e] $p}

    proc quad::K {a b} {set a}

    proc quad::MALLOC {} {
        variable nextedge
        K [set nextedge] [incr nextedge 4]
    }

    proc quad::FREE e {
     variable next
     variable org
     variable mark

     foreach i [list $e [ROT $e] [SYM $e] [TOR $e]] {
      array unset next $i
      array unset org $i
      array unset mark $i
     }
    }

    proc quad::MakeEdge {} {
        set e [EDGE [MALLOC]]
        SETNEXT $e $e
        SETNEXT [SYM $e] [SYM $e]
        SETNEXT [ROT $e] [TOR $e]
        SETNEXT [TOR $e] [ROT $e]
        foreach i [list $e [ROT $e] [SYM $e] [TOR $e]] {
           SETMARK $i 0
        }
        return $e
    }

    proc quad::Splice {a b} {
        set alpha [ROT [ONEXT $a]]
        set beta  [ROT [ONEXT $b]]
        set ta [ONEXT $a]
        set tb [ONEXT $b]
        SETNEXT $a $tb
        SETNEXT $b $ta
        set ta [ONEXT $alpha]
        set tb [ONEXT $beta]
        SETNEXT $alpha $tb
        SETNEXT $beta $ta
    }

    proc quad::DestroyEdge e {
        set f [SYM $e]
        if {[ONEXT $e] != $e} {Splice $e [OPREV $e]}
        if {[ONEXT $f] != $f} {Splice $f [OPREV $f]}
        FREE [EDGE $e]
    }

    # Recursive DFS
    # can overflow stack!

    proc quad::DoEnum {e visit closure mark} {
        while {[MARK [EDGE $e]] != $mark} {
          SETMARK [EDGE $e] $mark
          $visit $e $closure
            DoEnum [RPREV $e] $visit $closure $mark
            set e [ONEXT $e]
        }
    }

    proc quad::enum {a visit closure} {
        variable nextmark
        set mark $nextmark
        incr nextmark
        if {$nextmark==0} {set nextmark 1}
        DoEnum $a $visit $closure $mark
    }

    # non-recursive bfs
    # but is it right?

    proc quad::bfs {e visit closure} {
      variable nextmark
      set mark $nextmark
      incr nextmark
      if {$nextmark==0} { set nextmark 1}
      set q [list $e]
      while {[llength $q]} {
        set n [lindex $q 0]
        set q [lreplace $q 0 0]
        if {[MARK [EDGE $n]] != $mark} {
          SETMARK [EDGE $n] $mark
          $visit $n $closure
          set t [SYM $n]
          while {1} {
            set t [ONEXT $t]
            if {$t == [SYM $n]} break
            if {[MARK [EDGE $t]] != $mark} {   
              lappend q $t
            }
          }
        }
      }
    }
      

    ########################
    # Delaunay Triangulation       
    ########################

    proc quad::MAP p {
      variable index
      lindex $index $p
    }

    proc quad::COORD p {
      variable sites
      lindex $sites $p
    }

    proc quad::CCW {a b c} {
        set s1 [COORD $a]
        set s2 [COORD $b]
        set s3 [COORD $c]
        foreach {x1 y1} $s1 {x2 y2} $s2 {x3 y3} $s3 break
        expr {($x2*$y3-$y2*$x3)-($x1*$y3-$y1*$x3)+($x1*$y2-$y1*$x2)}
    }

    proc quad::RightOf {s e} {
        expr {[CCW $s [DEST $e] [ORG $e]] > 0.0}
    }

    proc quad::LeftOf {s e} {
        expr {[CCW $s [ORG $e] [DEST $e]] > 0.0}
    }

    proc quad::InCircle {a b c d} {
        set s1 [COORD $a]
        set s2 [COORD $b]
        set s3 [COORD $c]
        set s4 [COORD $d]
        foreach {x1 y1} $s1 {x2 y2} $s2 {x3 y3} $s3 {x4 y4} $s4 break
        expr {(($y4-$y1)*($x2-$x3)+($x4-$x1)*($y2-$y3)) * \
              (($x4-$x3)*($x2-$x1)-($y4-$y3)*($y2-$y1)) > \
              (($y4-$y3)*($x2-$x1)+($x4-$x3)*($y2-$y1)) * \
              (($x4-$x1)*($x2-$x3)-($y4-$y1)*($y2-$y3))}
    }

    proc quad::Connect {a b} {
        set e [MakeEdge]
        SETORG $e [DEST $a] 
        SETDEST $e [ORG $b]
        Splice $e [LNEXT $a]
        Splice [SYM $e] $b
        return $e
    } 

    # Divide and Conquer

    proc quad::RecDelaunay {sl sh} {

        if {$sh == $sl+2} {
            set a [MakeEdge]
            SETORG $a [MAP $sl]
            SETDEST $a [MAP [incr sl]]
            return [list $a [SYM $a]]
        } elseif {$sh == $sl+3} {
            set s1 [MAP $sl]
            set s2 [MAP [incr sl]]
            set s3 [MAP [incr sl]]
            set a [MakeEdge]
            set b [MakeEdge]
            set ct [CCW $s1 $s2 $s3]
            Splice [SYM $a] $b
            SETORG $a $s1
            SETDEST $a $s2
            SETORG $b $s2
            SETDEST $b $s3
            if {$ct == 0.0} {
                return [list $a [SYM $b]]
            } else {
                set c [Connect $b $a]
                if {$ct > 0.0} {
                    return [list $a [SYM $b]]
                } else {
                    return [list [SYM $c] $c]
                }
            }
        } else {
            set sm [expr {($sl+$sh)/2}]
            foreach {ldo ldi} [RecDelaunay $sl $sm] break
            foreach {rdi rdo} [RecDelaunay $sm $sh] break

            while {1} {
                if {[LeftOf [ORG $rdi] $ldi]} {
                    set ldi [LNEXT $ldi]
                } elseif {[RightOf [ORG $ldi] $rdi]} {
                    set rdi [RPREV $rdi]
                } else break
            }

            set base1 [Connect [SYM $rdi] $ldi]
            if {[ORG $ldi] == [ORG $ldo]} { set ldo [SYM $base1]}
            if {[ORG $rdi] == [ORG $rdo]} { set rdo $base1}

            while {1} {

                set lcand [ONEXT [SYM $base1]]
                if {[RightOf [DEST $lcand] $base1]} {
                   while {[InCircle \
                       [DEST $base1] [ORG $base1] [DEST $lcand] \
                     [DEST [ONEXT $lcand]]]} {
                           set t [ONEXT $lcand]
                           DestroyEdge $lcand
                           set lcand $t
                   }
                }

                set rcand [OPREV $base1]
                if {[RightOf [DEST $rcand] $base1]} {
                   while {[InCircle \
                       [DEST $base1] [ORG $base1] [DEST $rcand] \
                         [DEST [OPREV $rcand]]]} {
                           set t [OPREV $rcand]
                           DestroyEdge $rcand
                           set rcand $t
                   }
                }

                if {![RightOf [DEST $lcand] $base1] && \
                    ![RightOf [DEST $rcand] $base1]} {break}

                if {![RightOf [DEST $lcand] $base1] || \
                    ([RightOf [DEST $rcand] $base1] && \
                     [InCircle \
                       [DEST $lcand] [ORG $lcand] [ORG $rcand] \
                       [DEST $rcand]])} {
                    set base1 [Connect $rcand [SYM $base1]]
                } else {
                    set base1 [Connect [SYM $base1] [SYM $lcand]]
                }
            }
            return [list $ldo $rdo]
        }
    }

    proc quad::Ptcmp {a b} {
      foreach {x1 y1} $a {x2 y2} $b break
      if {$x1==$x2} {
         set f [expr {$y2-$y1}]
      } else {
         set f [expr {$x2-$x1}]
      }
      if {$f==0.0} {
          return 0
      } elseif {$f>0.0} {
          return -1 
      } else {
          return 1
      }
    }

    proc quad::delaunay points {
        variable sites
        variable index
        variable le
        variable re
        reset
        set sites $points
        set index [lsort -indices -command quad::Ptcmp $points]
        foreach {le re} [RecDelaunay 0 [llength $sites]] break
        return $le
    }

    ######################
    # some plotting stuff
    ######################
    package require Tk

    proc quad::calctransform c {
        variable sites
        variable scale
        variable fx
        variable fy
        set vy [$c cget -height]
        set vx [$c cget -width]
        set s [lindex $sites 0]
        foreach {x y} $s break
        set minx $x; set maxx $x
        set miny $y; set maxy $y
        foreach s $sites {
            foreach {x y} $s break
            if {$x<$minx} {set minx $x}
            if {$x>$maxx} {set maxx $x}
            if {$y<$miny} {set miny $y}
            if {$y>$maxy} {set maxy $y}
        }
        set sx [expr {($vx-20)/($maxx-$minx)}]
        set sy [expr {($vy-20)/($maxy-$miny)}]
        set scale [expr {$sx>$sy?$sy:$sx}]
        set fx [expr {$vx/2-$scale*($minx+$maxx)/2}]
        set fy [expr {$vy/2+$scale*($miny+$maxy)/2}]
    }

    proc quad::trans p {
        variable scale
        variable fx
        variable fy
        foreach {x y} $p break
        set x [expr {$scale*$x+$fx}]
        set y [expr {$fy-$scale*$y}]
        return [list $x $y]
    }

    proc quad::drawsites c {
        variable sites
        set n 0
        foreach s $sites {
            set cs [trans $s]
            foreach {x y} $cs break 
            $c create oval [expr {$x-3}] [expr {$y-3}] \
                [expr {$x+3}] [expr {$y+3}] \
                -fill black 
            $c create text [expr {$x+9}] [expr {$y-0}] -text $n -fill red
            incr n
       } 
    }

    proc quad::drawedge {e c} {
        set s1 [COORD [ORG $e]]
        set s2 [COORD [DEST $e]]
        set cs1 [trans $s1]
        set cs2 [trans $s2]
        foreach {x1 y1} $cs1 {x2 y2} $cs2 break
        $c create line $x1 $y1 $x2 $y2 -fill blue
    }

    proc quad::drawedges {e c} {
        bfs $e quad::drawedge $c
    }

    proc quad::circumcenter {a b c} { 
      set sa [COORD $a]
      set sb [COORD $b]
      set sc [COORD $c]
      foreach {xa ya} $sa {xb yb} $sb {xc yc} $sc break

      set A00 [expr {2*($xa-$xb)}]
      set A01 [expr {2*($ya-$yb)}]
      set B0 [expr {($xa-$xb)*($xa+$xb) + ($ya-$yb)*($ya+$yb)}]
      set A10 [expr {2*($xa-$xc)}]
      set A11 [expr {2*($ya-$yc)}]
      set B1 [expr { ($xa-$xc)*($xa+$xc) + ($ya-$yc)*($ya+$yc)}]
      set det [expr {$A00*$A11 - $A01*$A10}]
      set vx [expr {($B0*$A11 - $B1*$A01)/$det}]
      set vy [expr {($A00*$B1 - $A10*$B0)/$det}]
        return [list $vx $vy]
    }

    proc quad::far_left {ap bp} {
      set RBIG 100
      foreach {xa ya} [COORD $ap] {xb yb} [COORD $bp] break
      set xm [expr {($xa + $xb)/2}]
      set ym [expr {($ya + $yb)/2}]
      set xd [expr {($xb - $xa)}]
      set yd [expr {($yb - $ya)}]
      set d [expr {sqrt($xd*$xd+$yd*$yd)}]
      set vx [expr {$xm-$RBIG*$yd/$d}]
      set vy [expr {$ym+$RBIG*$xd/$d}]
      return [list $vx $vy]
    }

    proc quad::left_voronoi_vertex e { 
      set ap [ORG $e]
      set bp [DEST $e]
      set cp [DEST [LNEXT $e]]
      if {[CCW $ap $bp $cp] > 0.0} {     
         return [circumcenter $ap $bp $cp]
      } else { 
        return [far_left $ap $bp]
      }
    }

    proc quad::drawdualedge {e c} {
        set s1 [left_voronoi_vertex $e]
        set s2 [left_voronoi_vertex [SYM $e]]
        set cs1 [trans $s1]
        set cs2 [trans $s2]
        foreach {x1 y1} $cs1 {x2 y2} $cs2 break
        $c create line $x1 $y1 $x2 $y2 -fill green
    }

    proc quad::drawdualedges {e c} {
        bfs $e quad::drawdualedge $c
    }


    ########
    # lists 
    ########

    proc quad::Edge {e c} {
        variable counter
        lappend counter [list [ORG $e] [DEST $e]]
    }

    proc quad::edges e {
      variable counter [list]
      bfs $e quad::Edge {} 
      return $counter
    }

    proc quad::Faces {e f} {
      variable counter
      set mark [MARK [EDGE $e]]
      foreach i [list $e [SYM $e]] {
        if {[MARK [TOR $i]]==$mark} continue
        set a [ORG $i]
        set b [DEST $i]
        set c [DEST [LNEXT $i]]
        if {[CCW $a $b $c] > 0} {
          lappend counter [list $a $b $c]
          SETMARK [TOR $i] $mark
          SETMARK [TOR [LNEXT $i]] $mark
          SETMARK [ROT [ONEXT $i]] $mark
        }
      }
    }
            
    proc quad::faces e {
      variable counter [list] 
      bfs $e quad::Faces {}
      return $counter
    }

    #########
    # Demo
    #########

    proc makesites {n} {
      set list {}
      for {set i 0} {$i<$n} {incr i} {
         set x [expr {rand()}]
         set y [expr {rand()}]
        lappend list [list $x $y]
      }
      return $list
    }

    proc lists e {
     set lf [quad::faces $e]
     puts "[llength $lf] faces=$lf"
     set le [quad::edges $e]
     puts "[llength $le] edges=$le" 
    }

    proc plot {e c} {
        $c delete all
        quad::calctransform $c
        quad::drawsites $c
        quad::drawedges $e $c
        quad::drawdualedges $e $c
    }

    proc demo {c} {
      set l [makesites 10]
      set e [quad::delaunay $l]
      plot $e $c
      lists $e
    }

    frame .x
    pack .x
    canvas .x.c -bg white
    frame .x.f
    button .x.f.q -text QUIT -command exit
    button .x.f.d -text DEMO -command [list demo .x.c]
    pack .x.c .x.f -fill x -expand 1
    pack .x.f.q .x.f.d -side left
    raise .x
    expr {srand(5)}
    demo .x.c


----
!!!!!!
%| [Category Mathematics] |%
!!!!!!