Arjen Markus (5 may 2004) I wrote this script to provide yet another interpolation method, this time focussing on two-dimensional problems. You can in principle extend it to an arbitrary number of dimensions but it will get tedious and complicated very fast.
Note that the name "Digital Elevation Model" or "Digital Terrain Model" refers to a technique used by geographers to describe the height (elevation above a reference level) of the terrain they are studying. In this script, the parameter can be any (scalar) quantity.
A prerequisite for using it is that the data points have been organised in non-overlapping triangles. An external program that does this quickly and reliably is qhull (http://www.qhull.org ).
Enjoy
# digelev.tcl --
#
# Interpolation in two dimensions using a "digital elevation model"
#
# The interpolation technique used here relies on the presence of
# arbitrarily positioned data points and a network of triangles.
# Each data point is the vertex of at least one such triangle. Within
# each triangle the quantity is supposed to vary linearly with the
# two coordinates.
# Determining an appropriate triangulation is not trivial, but with
# software like qhull (http://www.qhull.org) the burden is put on
# an external program. Interfacing to that program (or if need be
# to the library) is easy enough. Here we assume this has been
# done.
# The data are organised as follows:
# - the data points are stored as x, y, z-value
# - the triangles are stored as the indices of three data points
# (the index being the index into the list of data points)
# This script does not optimise the search for the point to be
# interpolated: you can think of a number of ways to do that,
# but the basic operations remain the same:
# - Determine which triangle (if any) the point is contained in
# - Interpolate the z values on the three vertices
# (Optimisation makes sense if you have a lot of points to
# interpolate, for instance, if you are going to interpolate on a
# rectangular grid)
#
namespace eval ::math::interpolation {
}
# InterpolatePlane --
# Interpolate the z values given the three vertices
# Arguments:
# datapoints The list of data points
# vertices The indices of the three vertices
# x The X coordinate of the point in question
# y The Y coordinate of the point in question
# Returns:
# zvalue - the estimated value at the given point
#
proc ::math::interpolation::InterpolatePlane {datapoints vertices x y} {
foreach {v1 v2 v3} $vertices {break}
foreach {x1 y1 z1} [lindex $datapoints $v1] {break}
foreach {x2 y2 z2} [lindex $datapoints $v2] {break}
foreach {x3 y3 z3} [lindex $datapoints $v3] {break}
set x12 [expr {$x2-$x1}]
set y12 [expr {$y2-$y1}]
set z12 [expr {$z2-$z1}]
set x13 [expr {$x3-$x1}]
set y13 [expr {$y3-$y1}]
set z13 [expr {$z3-$z1}]
set dx [expr {$x-$x1}]
set dy [expr {$y-$y1}]
#
# Parametrisation of the plane:
#
# /x\ /x1\ /x12\ /x13\
# |y| = |y1| + L |y12| + M |y13|
# \z/ \z1/ \z12/ \z13/
#
# Now: find the values of L and M for (x,y) and fill in
#
# L (x12*y13-x13*y12) = (x-x1)*y13 - (y-y1)*x13
# M (x13*y12-x12*y13) = (x-x1)*y12 - (y-y1)*x12
#
set det [expr {$x12*$y13-$x13*$y12}]
set L [expr {($dx*$y13-$dy*$x13)/$det}]
set M [expr {-($dx*$y12-$dy*$x12)/$det}]
set z [expr {$z1+$z12*$L+$z13*$M}]
return $z
}
# InsideTriangle --
# Find out if the point lies in the given triangle
# Arguments:
# datapoints The list of data points
# vertices The indices of the triangle
# x The X coordinate of the point in question
# y The Y coordinate of the point in question
# Returns:
# 1 if inside, 0 otherwise
# Note:
# Optimised for triangles!
#
proc ::math::interpolation::InsideTriangle {datapoints vertices x y} {
foreach {v1 v2 v3} $vertices {break}
foreach {x1 y1 z1} [lindex $datapoints $v1] {break}
foreach {x2 y2 z2} [lindex $datapoints $v2] {break}
foreach {x3 y3 z3} [lindex $datapoints $v3] {break}
set xmin [expr {$x1 < $x2 ? $x1 : $x2}]
set xmin [expr {$xmin < $x3 ? $xmin : $x3}]
set xmax [expr {$x1 > $x2 ? $x1 : $x2}]
set xmax [expr {$xmax > $x3 ? $xmax : $x3}]
set ymin [expr {$y1 < $y2 ? $y1 : $y2}]
set ymin [expr {$ymin < $y3 ? $ymin : $y3}]
set ymax [expr {$y1 > $y2 ? $y1 : $y2}]
set ymax [expr {$ymax > $y3 ? $ymax : $y3}]
if { $x < $xmin } {return 0}
if { $x > $xmax } {return 0}
if { $y < $ymin } {return 0}
if { $y > $ymax } {return 0}
#
# Within the bounding box, so now for some more
# serious work:
# The point must lie on the same side of each edge
# spanned by two vertices as the third vertex
#
foreach {xn1 yn1 xn2 yn2 xn3 yn3} \
[list $x1 $y1 $x2 $y2 $x3 $y3 \
$x2 $y2 $x3 $y3 $x1 $y1 \
$x3 $y3 $x1 $y1 $x2 $y2 ] {
set inpp [expr {($x -$xn1)*($yn2-$yn1)-($y -$yn1)*($xn2-$xn1)}]
set inp3 [expr {($xn3-$xn1)*($yn2-$yn1)-($yn3-$yn1)*($xn2-$xn1)}]
if { $inpp < 0.0 && $inp3 > 0.0 } {return 0}
if { $inpp > 0.0 && $inp3 < 0.0 } {return 0}
}
#
# The point lies inside the triangle ...
#
return 1
}
# WhichTriangle --
# Find out which triangle a given point lies in
# Arguments:
# datapoints The list of data points
# network The network of triangles
# x The X coordinate of the point in question
# y The Y coordinate of the point in question
# Returns:
# idx - the index of the triangle or {} if outside all
#
proc ::math::interpolation::WhichTriangle {datapoints network x y} {
set idx 0
foreach triangle $network {
if { [InsideTriangle $datapoints $triangle $x $y] } {
return $idx
}
incr idx
}
# Not found, then return an empty value
#
return {}
}
# interpolateTriangles --
# Interpolate using a triangulation network
# Arguments:
# datapoints The list of data points
# network The network of triangles
# x The X coordinate of the point in question
# y The Y coordinate of the point in question
# Returns:
# Interpolated value for the given point or {} if outside
# any triangle
#
proc ::math::interpolation::interpolateTriangles {datapoints network x y} {
set idx [WhichTriangle $datapoints $network $x $y]
set value {}
if { $idx != {} } {
set vertices [lindex $network $idx]
set value [InterpolatePlane $datapoints $vertices $x $y]
}
return $value
}
# main --
# Test the procedures. The data points are arranged like this:
#
# 2
#
# 4 0 1
#
# 3
#
# The first set consists of five points on a sloping plane
# (equation: x+y=1), the second looks like a four-sided pyramid
#
set datapoints {
{ 0 0 0 }
{ 1 0 1 }
{ 0 1 1 }
{ 0 -1 -1 }
{ -1 0 -1 }
}
set hatdata {
{ 0 0 1 }
{ 1 0 0 }
{ 0 1 0 }
{ 0 -1 0 }
{ -1 0 0 }
}
set vertices { 0 1 2 }
set network {
{ 0 1 2 }
{ 0 1 3 }
{ 0 4 3 }
{ 0 4 2 } }
puts "[::math::interpolation::InterpolatePlane $datapoints $vertices 0 0] - expected: 0"
puts "[::math::interpolation::InterpolatePlane $datapoints $vertices 1 0] - expected: 1"
puts "[::math::interpolation::InterpolatePlane $datapoints $vertices 0 1] - expected: 1"
puts "[::math::interpolation::InterpolatePlane $datapoints $vertices 0.5 0.5] - expected: 1"
puts "[::math::interpolation::interpolateTriangles \
$datapoints $network 0 0] - expected: 0"
puts "[::math::interpolation::interpolateTriangles \
$datapoints $network 1 1] - expected: --"
puts "[::math::interpolation::interpolateTriangles \
$datapoints $network 0.5 0.5] - expected: 1"
puts "[::math::interpolation::interpolateTriangles \
$datapoints $network 0.5 0.0] - expected: 0.5"
puts "[::math::interpolation::interpolateTriangles \
$hatdata $network 0.0 0.0] - expected: 1.0"
puts "[::math::interpolation::interpolateTriangles \
$hatdata $network 0.0 0.5] - expected: 0.5"
puts "[::math::interpolation::interpolateTriangles \
$hatdata $network -0.5 0.5] - expected: 0.0"
puts "[::math::interpolation::interpolateTriangles \
$hatdata $network -0.5 0.5001] - expected: --"
puts "[::math::interpolation::interpolateTriangles \
$hatdata $network -0.4999 0.5001] - expected: 0.0" ;# Tricky this one