3-D interpolation
- distance-weighted averaging
Distance-weighted averaging involves the calculation of weighted
averages of Z at the target points, where the weights are some
function of distance between a target point and nearby control
points. The number of control points to use is an issue; more
points will yield a smoother surface, but one that may not
"respect" the values of Z at the control points,
unless specific "fixups" are done.
# distance-weighted averaging
interpolation
x <- c( 8,38,46)
y <- c(12,42, 2)
z <- c( 5, 9, 1)
xtarg <- seq(1,50,1)
ytarg <- seq(1,45,1)
zhat <- matrix(0, 50, 45)
for (j in 1:45) {
for (i in 1:50) {
d1 <- sqrt((x[1]-xtarg[i])^2 + (y[1]-ytarg[j])^2)
d2 <- sqrt((x[2]-xtarg[i])^2 + (y[2]-ytarg[j])^2)
d3 <- sqrt((x[3]-xtarg[i])^2 + (y[3]-ytarg[j])^2)
zhat[i,j]=(z[1]/d1+z[2]/d2+z[3]/d3)/(1.0/d1+1.0/d2+1.0/d3)
if (d1 <= 0.00001) zhat[i,j] <- z[1]
if (d2 <= 0.00001) zhat[i,j] <- z[2]
if (d3 <= 0.00001) zhat[i,j] <- z[3]
}
}
image(xtarg,ytarg,zhat)
contour(xtarg,ytarg,zhat, nlevels=9, xlab="", ylab="", add=TRUE)
points(x, y); text(x+2, y, z)
title ("Distance-weighted averaging")
Another approach involves fitting a plane to three points, and
evaluating the fitted function at the control points contained within
the triangle formed by the points. A plane can be fit in R
using the lm (linear
model) function, the one most often used in regression analysis.
One problem that arises in fitting planes (not illustrated by this
example, however) is the sometimes abrupt changes in interpolated
values occur as one passes from one triangular region to another.
# fit a plane to three points
x <- c( 8,38,46)
y <- c(12,42, 2)
z <- c( 5, 9, 1)
xmar <- seq(1,50,1)
ymar <- seq(1,45,1)
target.mar <- list(x=xmar,y=ymar)
zpoly <- lm(z ~ x + y )
zpoly
zhat <- predict.lm(zpoly, newdata=expand.grid(target.mar))
zhat <- matrix(zhat, 50, 45)
image(xmar,ymar,zhat)
contour(xmar,ymar,zhat, nlevels=9, xlab="", ylab="", add=T)
points(x, y); text(x+2, y, z)
title ("First-order (planar) polynomial")
Higher-order (i.e. polynomial) functions can also be fit to a small
number of control points and then evaluated at the target
points, Higher-order models require more control points, and as
is the case with lower-order models, at times unusual interpolated
values may occur as the control points change across the region under
consideration.
# fit a quadratic polynomial to six
points
x <- c( 8,38,46,10,20,48)
y <- c(12,42, 2,27,32,15)
z <- c( 5, 9, 1, 7, 9, 3)
xmar <- seq(1,50,1)
ymar <- seq(1,45,1)
target.mar <- list(x=xmar,y=ymar)
zpoly <- lm(z ~ x * y + I(x^2) + I(y^2))
zpoly
zhat <- predict.lm(zpoly, newdata=expand.grid(target.mar))
zhat <- matrix(zhat, 50, 45)
image(xmar,ymar,zhat)
contour(xmar,ymar,zhat, nlevels=9, xlab="", ylab="", add=T)
points(x, y); text(x+2, y, z)
title ("Second-order (quadtratic) polynomial")
|