Triangulation

Richard Suchenwirth 2003-08-09 - Given a set of points in a plane, one can construct a planar graph so that all points (vertices) are connected, directly or indirectly, with edges, such that no two edges cross. A maximum planar graph is called a triangulation of the vertex set.

My playing with triangulations starts with an experimental UI, where you can seed vertices by left-clicking (right-click on one to remove it). Click on "Complete" to see the complete graph, each vertex connected with every other (gray edges); "Triangulate" to test the code below and see the triangulation in red.

Vertices in the graph are just named $x/$y by their coordinates; edges are named $from/$to where from and to are vertices.

 proc demo {} {
    canvas .c -bg white
    frame .f
    button .f.c  -text Clear       -command {.c delete all}
    button .f.co -text Complete    -command {showComplete .c}
    button .f.tr -text Triangulate -command {showTriangulate .c}
    eval pack [winfo children .f] -side left
    pack .c .f -fill x -expand 1
    bind .c <1> {addVertex %W %x %y}
    .c bind vertex <3> {%W delete current}
 }
 proc showComplete w {
    $w delete edge
    foreach edge [completeGraph [get vertex $w]] {
        showEdge $w $edge
    }
 }
 proc showEdge {w edge {fill gray}} {
    regexp {(.+)/(.+),(.+)/(.+)} $edge -> x0 y0 x1 y1
    set ::length($edge) [expr {hypot($x1-$x0,$y1-$y0)}]
    $w create line $x0 $y0 $x1 $y1 -tags "edge $edge" -fill $fill
 }
 proc get {tag w} {
    set res {}
    foreach v [$w find withtag $tag] {
        lappend res [lindex [$w gettags $v] 1]
    }
    set res
 }
 proc completeGraph vertices {
    set graph {}
    foreach i $vertices {
        foreach j $vertices {
            if {$i<$j} {lappend graph $i,$j}
        }
    }
    set graph
 }
 proc showTriangulate w {
    $w delete edge
    showComplete $w
    wm title . Wait...
    set t0 [clock clicks -milliseconds]
    foreach edge [triangulate [get edge $w]] {
        showEdge $w $edge red
    }
    wm title . [expr {[clock clicks -milliseconds] - $t0}]
 }

My idea of triangulation is to iterate over all pairs of edges and see if they cross - then remove the longer of them; until no more crossing edges are found. However, sometimes it deletes one edge too many... at least the convex hull is always kept right.

 proc triangulate graph {
    while 1 {
        set found 0
        foreach i $graph {
            foreach j $graph {
                if {$i!=$j && [crossing $i $j]} {
                    lremove graph [longer $i $j]
                    set found 1
                    break
                }
            }
            if $found break
        }
        if {!$found} break
    }
    set graph
 }

To detect whether two line segments cross each other took me a bit of heavy math:

 proc crossing {edge1 edge2} {
    regexp {(.+)/(.+),(.+)/(.+)} $edge1 -> x0 y0 x1 y1
    regexp {(.+)/(.+),(.+)/(.+)} $edge2 -> x2 y2 x3 y3
    if [adjacent $x0/$y0 $x1/$y1 $x2/$y2 $x3/$y3] {return 0}
    set m1 [slope $x0 $y0 $x1 $y1]
    set b1 [expr {$y0-$m1*$x0}]
    set m2 [slope $x2 $y2 $x3 $y3]    
    set b2 [expr {$y2-$m2*$x2}]
    set x [slope $m2 $b1 $m1 $b2]
    expr {[between $x0 $x $x1] && [between $x2 $x $x3]}
 }
 proc adjacent args {
    expr {[llength [lsort -unique $args]]<[llength $args]}
 }
 proc slope {x0 y0 x1 y1} {
    # slightly "bend" a vertical line, to avoid division by zero
    if {$x1==$x0} {set x1 [expr {$x1+0.00000001}]}
    expr {double($y1-$y0)/($x1-$x0)}
 }
 proc between {a b c} {
    expr {$b==[lindex [lsort -real [list $a $b $c]] 1]}
 }
 proc longer {edge1 edge2} {
    global length
    expr {$length($edge1) > $length($edge2)? $edge1: $edge2}
 }
 proc addVertex {w x y} {
    $w create rect [expr $x-2] [expr $y-2] [expr $x+2] [expr $y+2] \
        -tags "vertex $x/$y" -fill blue
 }
 proc lremove {varName element} {
    upvar 1 $varName var
    set pos [lsearch $var $element]
    set var [lreplace $var $pos $pos]
 }
 
 demo
 bind . <Escape> {exec wish $argv0 &; exit}
 bind . <F1> {console show}

Please comment if you know better (more reliable) ways of triangulation! I originally wanted to play with the Travelling Salesman problem on the triangulated graph, but then it was so hot, and soon the weekend was over...

DGP Use qhull, http://www.qhull.org/

RS Hm, yes, thanks, but I wanted to do it in pure-Tcl to find out how it can be done ;)

VH (29 Nov 2011) The processing time required for this intersection method grows exponentially as the number of points increases. The benchmark algorithm seem to be Fortune's method .

AM See the Digital Elevation Model page for an application of triangulation

AM (4 august 2005) I have thought of a relatively simple algorithm to automatically determine the triangulation. It will not guarantee anything about the form of the triangles, but it can be done in pure Tcl with little effort - AM see Simple triangulation algorithm}

EKB I played around with this a bit. I had an idea for providing an alternative routine for determining whether line segments cross. The motivation for creating an alternative algorithm is the step in the slope proc for taking care of a vertical line segment. The choice of coordinate system here is completely arbitrary, so "up" shouldn't have any meaning. There should be some way of deciding whether lines cross that doesn't make assumptions about the coordinate system.

After I worked this out I saw AM's Digital Elevation Model, which seems to do some similar manipulations.

Here's the approach:

1) Take two line segments

 x0, y0 to x0, y1
 u0, v0 to u1, v1

2) Introduce variables s and t that run from 0 to 1. To get any point on the first line segment, set

 x(s) = x0 + s * (x1 - x0)
 y(s) = y0 + s * (y1 - y0)

and on the second line segment:

 u(t) = u0 + t * (u1 - u0)
 v(t) = v0 + t * (v1 - v0)

Now the line segments cross whenever

 x(s) = u(t)
 y(s) = v(t)

for some s and t between 0 and 1.

3) Setting x(s) = u(t) and y(s) = v(t) gives a set of equations for s and t. After some rearrangement, it gives:

 (u1 - u0) * t - (x1 - x0) * s = x0 - u0
 (v1 - v0) * t - (y1 - y0) * s = y0 - v0

It's handy to replace (u1 - u0), etc. with new variables and write this as

 du * t - dx * s = xu
 dv * t - dy * s = yv

where du = u1 - u0, etc. This is a matrix equation that can be solved for s and t.

4) Solving for s and t (using matrix algebra tricks or brute force) gives

 -dy * xu + dx * yv = s * (dx * dv - du * dy)
 -dv * xu + du * yv = t * (dx * dv - du * dy)

For convenience, the left-hand-side can be called S and T, and the expression in parenthesis can be called D (because it's the determinant of a matrix). So the equation is

 S = s * D
 T = t * D

5) Since s and t can go from 0 to 1, then if the lines intersect, S and T must be between 0 and D. That gives the test

 [between 0 $S $D] && [between 0 $T $D]

Here's the new proc:

 proc crossing {edge1 edge2} {
    regexp {(.+)/(.+),(.+)/(.+)} $edge1 -> x0 y0 x1 y1
    regexp {(.+)/(.+),(.+)/(.+)} $edge2 -> u0 v0 u1 v1
    if [adjacent $x0/$y0 $x1/$y1 $u0/$v0 $u1/$v1] {return 0}
    set dx [expr {$x1 - $x0}]
    set dy [expr {$y1 - $y0}]
    set du [expr {$u1 - $u0}]
    set dv [expr {$v1 - $v0}]
    set xu [expr {$x0 - $u0}]
    set yv [expr {$y0 - $v0}]
    set D [expr {$dx * $dv - $du * $dy}]
    set S [expr {-$dy * $xu + $dx * $yv}]
    set T [expr {-$dv * $xu + $du * $yv}]
    expr {[between 0 $S $D] && [between 0 $T $D]}
 }

PD See my attempt at Delaunay Triangulation in pure-Tcl

HJG See also math or rather [L1 ] for ::math::geometry::lineSegmentsIntersect and ::math::geometry::findLineSegmentIntersection