Arjen Markus (15 august 2005) It was surprisingly simple to implement this algorithm to construct a network of non-overlapping triangles from a set of points. But then:
- There is no way that I can see to generalise it to 3D or more dimensions
- There is concern for the form of the triangles (near-collinearity can cause very narrow triangles)
Anyway, enjoy!
# triangulate.tcl --
# Construct a triangle network for a given set of points
# The algorithm is somewhat naive as that it does not
# guarantee any particular constraints on the triangles'
# form or other properties other than that they cover
# the convex hull of the point set, that each vertex
# is a point from the set and that they do not overlap.
# In particular nearly collinear points may cause very
# narrow triangles.
#
namespace eval ::math::triangulate {
variable torad [expr {180.0/(4.0*atan(1.0))}]
}
# triangles --
# Construct the network of triangles for a set of points
# Arguments:
# points List of coordinates (each point is a list of
# the x and y coordinates, e.g.
# {{1 1} {2 3} {1 4} {5 1}}
# Result:
# A list of triangles (a triangle is represented as
# a list of three points)
#
proc ::math::triangulate::triangles {points} {
if { [llength $points] <= 2 } {
return -code error "The set of points must have at least three points"
}
#
# The trivial case
#
if { [llength $points] == 3 } {
return [list $points]
}
set triangles {}
#
# First step:
# Find the topmost point
#
set idx -1
set idxmax 0
set ymax [lindex $points 0 1]
foreach p $points {
incr idx
set y [lindex $points 1]
if { $y > $ymax } {
set idxmax $idx
set ymax $y
}
}
#
# Second step:
# Find the two adjacent points on the convex hull
#
set idxleft [FindHullPoint $points $idxmax left]
set idxright [FindHullPoint $points $idxmax right]
#
# Third step:
# Determine which points are inside the triangle
# formed by (idxmax, idxleft, idxright) and form
# a convex arc
#
set innerpoints [FindArcPoints $points $idxmax $idxleft $idxright]
puts "$idxmax $idxleft $idxright: $innerpoints"
# Fourth step:
# Construct the triangles
#
foreach idx $innerpoints {
lappend triangles \
[list [lindex $points $idxmax] \
[lindex $points $idxleft] \
[lindex $points $idx] ]
set idxleft $idx
}
lappend triangles \
[list [lindex $points $idxmax] \
[lindex $points $idxleft] \
[lindex $points $idxright] ]
#
# Now, the last step, remove the point "idxmax" and
# repeat the procedure. Construct the list of all
# triangles found this way.
#
set triangles [concat $triangles \
[triangles [lreplace $points $idxmax $idxmax]]]
return $triangles
}
# FindHullPoint --
# Find the point on the convex hull adjacent to the given one
# Arguments:
# points List of coordinates
# idxc Index of the given point
# dir What direction (left or right)
# Result:
# Index of the point that was sought
#
proc ::math::triangulate::FindHullPoint {points idxc dir} {
foreach {xc yc} [lindex $points $idxc] {break}
set idx -1
set anglehull 360.0
if { $dir == "right" } {
set anglehull 0.0
}
set idxhull -1
foreach p $points {
incr idx
if { $idx == $idxc } {
continue
}
foreach {xp yp} [lindex $points $idx] {break}
set angle [Angle $xc $yc $xp $yp]
if { $dir == "left" } {
if { $angle < $anglehull } {
set idxhull $idx
set anglehull $angle
}
} else {
if { $angle > $anglehull } {
set idxhull $idx
set anglehull $angle
}
}
}
return $idxhull
}
# Angle --
# Determine the angle of a vector (in degrees)
# Arguments:
# x1 X coordinate of start
# y1 Y coordinate of start
# x2 X coordinate of end
# y2 Y coordinate of end
# Result:
# Angle to the x-axis
#
proc ::math::triangulate::Angle {x1 y1 x2 y2} {
variable torad
if { $x2 != $x1 || $y2 != $y1 } {
set angle [expr {$torad*atan2($y2-$y1,$x2-$x1)}]
# Corner case - important mainly for tests
if { $y1 == $y2 && $x2 > $x1 } {
set angle 360.0
}
} else {
set angle 0.0 ;# Hm, this is a problem! Better raise an error
error "Points are equal!"
}
if { $angle < 0.0 } {
set angle [expr {360.0+$angle}]
}
return $angle
}
# FindArcPoints --
# Find the points inside the given triangle that form a convex arc
# from the left to the right
# Arguments:
# points List of coordinates
# idxmax Index of the topmost point of the triangle
# idxleft Index of the leftmost point of the triangle
# idxright Index of the rightmost point of the triangle
# Result:
# List of indices
#
proc ::math::triangulate::FindArcPoints {points idxmax idxleft idxright} {
#
# The procedure is very simple:
# - Remove the topmost point from the set
# - Look for points on the convex hull of this reduced set,
# starting at the left point of the triangle
# - Continue until we reach the right point
#
# Note: here is an opportunity for optimisation - we
# need the reduced set later on.
#
set points [lreplace $points $idxmax $idxmax]
set end 0
set next $idxleft
set innerpoints {}
while { ! $end } {
set p [FindHullPoint $points $next right]
if { $p != $idxright && $p != $idxleft } {
lappend innerpoints [expr {$p>=$idxmax? $p+1 : $p}]
set next $p
} else {
set end 1
}
}
return $innerpoints
}
# test --
# Some points
#
set points {{0 1} {-1 0} {1 0} {0 -1} {-2 -1} {2 -1}}
puts [::math::triangulate::triangles $points]
