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]
lassign $s1 x1 y1
lassign $s2 x2 y2
lassign $s3 x3 y3
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]
lassign $s1 x1 y1
lassign $s2 x2 y2
lassign $s3 x3 y3
lassign $s4 x4 y4
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}]
lassign [RecDelaunay $sl $sm] ldo ldi
lassign [RecDelaunay $sm $sh] rdi rdo
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} {
lassign $a x1 y1
lassign $b x2 y2
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]
lassign [RecDelaunay 0 [llength $sites]] le re
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]
lassign $s x y
set minx $x; set maxx $x
set miny $y; set maxy $y
foreach s $sites {
lassign $s x y
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
lassign $p x y
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]
lassign $cs x y
$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]
lassign $cs1 x1 y1
lassign $cs2 x2 y2
$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]
lassign $sa xa ya
lassign $sb xb yb
lassign $sc xc yc
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
lassign [COORD $ap] xa ya
lassign [COORD $bp] xb yb
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]
lassign $cs1 x1 y1
lassign $cs2 x2 y2
$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