Commit d9dd2590 authored by agebhard's avatar agebhard
Browse files

* Finally remove AMS 526 code, completely switch to AMS 761 for irregular

  grids.
* Now both inter{p|pp}.{old|new} print warnings and divert to inter{p|pp}.
* ACM 761 code updated to the revised version by Renka and Brown from 1998.
* Fixed two array overruns in ACM 761 code (work array size increased,
  thin triangle removing step cancelled if necessary).
* Thin triangle removing can now be omitted by paramater remove=FALSE
  (this is default for linear interpolation).
parent b9de5f5f
...@@ -2,6 +2,7 @@ orig ...@@ -2,6 +2,7 @@ orig
config.log config.log
config.status config.status
autom4te.cache autom4te.cache
src-old
src-c src-c
src-i386 src-i386
build-debug.sh build-debug.sh
......
...@@ -3,13 +3,14 @@ http://www.acm.org/publications/policies/software-copyright-notice ...@@ -3,13 +3,14 @@ http://www.acm.org/publications/policies/software-copyright-notice
In detail: In detail:
1. Fortran code (src/*.f): 1. Fortran code (src/*.f except src/bilinear.f):
Copyrighted and Licensed by ACM, Copyrighted and Licensed by ACM,
see http://www.acm.org/publications/policies/software-copyright-notice see http://www.acm.org/publications/policies/software-copyright-notice
2. R interface (src/init.c src/akima.h R/*.R man/*.Rd data/*): 2. R interface (src/init.c src/akima.h R/*.R man/*.Rd data/*
except */bilinear.*)
The R interface code has been developed as work based on the The R interface code has been developed as work based on the
ACM licensed code, hence it is also ACM licensed, copyright ACM licensed code, hence it is also ACM licensed, copyright
...@@ -21,6 +22,7 @@ and to fulfill this, the modified work including the R interface ...@@ -21,6 +22,7 @@ and to fulfill this, the modified work including the R interface
is available free to secondary users, and no charge is associated is available free to secondary users, and no charge is associated
with such copies. with such copies.
3. Bilinear interpolation (*/bilinear*): Copyright A. Gebhardt,
free to reuse whithout restrictions.
...@@ -19,9 +19,12 @@ bilinear <- function(x,y,z,x0,y0){ ...@@ -19,9 +19,12 @@ bilinear <- function(x,y,z,x0,y0){
as.double(z), as.double(z),
as.integer(nx), as.integer(nx),
as.integer(ny), as.integer(ny),
ier=integer(1),
PACKAGE="akima") PACKAGE="akima")
if(ret$ier==1)
list(x=x0,y=y0,z=ret$z0) stop("duplicate coordinates in input grid!")
else
list(x=x0,y=y0,z=ret$z0)
} }
......
...@@ -4,7 +4,8 @@ interp <- ...@@ -4,7 +4,8 @@ interp <-
yo = seq(min(y), max(y), length = ny), linear = TRUE, yo = seq(min(y), max(y), length = ny), linear = TRUE,
extrap = FALSE, duplicate = "error", dupfun = NULL, extrap = FALSE, duplicate = "error", dupfun = NULL,
nx=40, ny=40, nx=40, ny=40,
jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE) jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE,
remove = !linear)
{ {
## handle sp data, save coordinate and value names ## handle sp data, save coordinate and value names
is.sp <- FALSE is.sp <- FALSE
...@@ -25,12 +26,6 @@ interp <- ...@@ -25,12 +26,6 @@ interp <-
} else } else
stop("either x,y,z are numerical or x is SpatialPointsDataFrame and z a name of a data column in x") stop("either x,y,z are numerical or x is SpatialPointsDataFrame and z a name of a data column in x")
} }
## FIXME: drop old akima code
if(linear)
ret <- interp.old(x,y,z,xo,yo,ncp=0,extrap=FALSE,
duplicate=duplicate,dupfun=dupfun)
else {
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z)))) if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed") stop("missing values and Infs not allowed")
...@@ -77,6 +72,12 @@ interp <- ...@@ -77,6 +72,12 @@ interp <-
## storage.mode(zo) <- "double" ## storage.mode(zo) <- "double"
miss <- !extrap # if not extrapolating, set missing values miss <- !extrap # if not extrapolating, set missing values
## Defaults from Fortran code for triangle removal moved to R code,
## see SDTRTT:
hbrmn <- 0.1 # height to base ratio for thin triangles
nrrtt <- 5 # recursion depth for triangle removal, 0 disables
if(!remove) nrrtt <- 0
ans <- .Fortran("sdsf3p", ans <- .Fortran("sdsf3p",
md = as.integer(1), md = as.integer(1),
ndp = as.integer(n), ndp = as.integer(n),
...@@ -89,17 +90,21 @@ interp <- ...@@ -89,17 +90,21 @@ interp <-
y = as.double(yo), y = as.double(yo),
z = as.double(matrix(0,nx,ny)), z = as.double(matrix(0,nx,ny)),
ier = integer(1), ier = integer(1),
wk = double(36 * n), wk = double(36 * 5* n),
iwk = integer(25 * n), iwk = integer(39 * n),
extrap = as.logical(matrix(extrap,nx,ny)), extrap = as.logical(matrix(extrap,nx,ny)),
near = integer(n),
nxt = integer(n),
dist = double(n),
linear = as.logical(linear), linear = as.logical(linear),
hbrmn = as.double(hbrmn),
nrrtt = as.integer(nrrtt),
PACKAGE = "akima") PACKAGE = "akima")
if(miss) if(miss)
ans$z[ans$extrap] <- NA ans$z[ans$extrap] <- NA
## Error code 11 from sdsf3p indicates failure of thin triangle removal.
if(ans$ier==11){
stop("removal of thin triangles from border failed. try to re-run with remove=FALSE")
}
## Error code 10 from sdsf3p indicates error code -2 from trmesh: ## Error code 10 from sdsf3p indicates error code -2 from trmesh:
## first three points collinear. ## first three points collinear.
## Try to add jitter to data locations to avoid collinearities, ## Try to add jitter to data locations to avoid collinearities,
...@@ -150,13 +155,12 @@ interp <- ...@@ -150,13 +155,12 @@ interp <-
y = as.double(yo), y = as.double(yo),
z = as.double(matrix(0,nx,ny)), z = as.double(matrix(0,nx,ny)),
ier = integer(1), ier = integer(1),
double(36 * n), double(36 * 5* n),
integer(25 * n), integer(39 * n),
extrap = as.logical(matrix(extrap,nx,ny)), extrap = as.logical(matrix(extrap,nx,ny)),
near = integer(n),
nxt = integer(n),
dist = double(n),
linear = as.logical(linear), linear = as.logical(linear),
hbrmn = as.double(hbrmn),
nrrtt = as.integer(nrrtt),
PACKAGE = "akima") PACKAGE = "akima")
if(miss) if(miss)
ans$z[ans$extrap] <- NA ans$z[ans$extrap] <- NA
...@@ -170,8 +174,8 @@ interp <- ...@@ -170,8 +174,8 @@ interp <-
} }
## prepare return value ## prepare return value
if(is.sp){ if(is.sp){
zm <- dim(ans$z)[1] zm <- nx
zn <- dim(ans$z)[2] zn <- ny
zvec <- c(ans$z) zvec <- c(ans$z)
xvec <- c(matrix(rep(ans$x,zn),nrow=zm,ncol=zn,byrow=FALSE)) xvec <- c(matrix(rep(ans$x,zn),nrow=zm,ncol=zn,byrow=FALSE))
yvec <- c(matrix(rep(ans$y,zm),nrow=zm,ncol=zn,byrow=TRUE)) yvec <- c(matrix(rep(ans$y,zm),nrow=zm,ncol=zn,byrow=TRUE))
...@@ -184,6 +188,5 @@ interp <- ...@@ -184,6 +188,5 @@ interp <-
} else { } else {
ret <- list(x=ans$x,y=ans$y,z=matrix(ans$z,nx,ny)) ret <- list(x=ans$x,y=ans$y,z=matrix(ans$z,nx,ny))
} }
} ## END FIXME
ret ret
} }
...@@ -2,76 +2,9 @@ ...@@ -2,76 +2,9 @@
yo = seq(min(y), max(y), length = 40), yo = seq(min(y), max(y), length = 40),
ncp = 0, extrap = FALSE, duplicate = "error", dupfun = NULL) ncp = 0, extrap = FALSE, duplicate = "error", dupfun = NULL)
{ {
warning("interp.old() is deprecated, use interp()")
## warning("interp.old() is deprecated, future versions will only provide interp()") interp(x, y, z,
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z)))) xo,
stop("missing values and Infs not allowed") yo, linear=(ncp==0),
if(ncp > 25) { extrap, duplicate, dupfun)
ncp <- 25
cat("ncp too large, using ncp=25\n")
}
drx <- diff(range(x))
dry <- diff(range(y))
if(drx == 0 || dry == 0)
stop("all data collinear") # other cases caught in Fortran code
if(drx/dry > 10000 || drx/dry < 0.0001)
stop("scales of x and y are too dissimilar")
n <- length(x)
nx <- length(xo)
ny <- length(yo)
if(length(y) != n || length(z) != n)
stop("Lengths of x, y, and z do not match")
xy <- paste(x, y, sep = ",")# trick for 'duplicated' (x,y)-pairs
if(duplicate == "error") {
if(any(duplicated(xy)))
stop("duplicate data points: need to set 'duplicate = ..' ")
}
else { ## duplicate != "error"
i <- match(xy, xy)
if(duplicate == "user")
dupfun <- match.fun(dupfun)#> error if it fails
ord <- !duplicated(xy)
if(duplicate != "strip") {
centre <- function(x)
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x))
z <- unlist(lapply(split(z,i), centre))
} else {
z <- z[ord]
}
x <- x[ord]
y <- y[ord]
n <- length(x)
}
zo <- matrix(0, nx, ny)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
misso <- matrix(miss, nx, ny)
if(extrap & ncp == 0)
warning("Cannot extrapolate with linear option")
ans <- .Fortran("idsfft",
as.integer(1),
as.integer(ncp),
as.integer(n),
as.double(x),
as.double(y),
as.double(z),
as.integer(nx),
as.integer(ny),
x = as.double(xo),
y = as.double(yo),
z = zo,
integer((31 + ncp) * n + nx * ny),
double(5 * n),
misso = as.logical(misso),
PACKAGE = "akima")[c("x", "y", "z", "misso")]
ans$z[ans$misso] <- NA
ans[c("x", "y", "z")]
} }
"interpp"<-function(x, y=NULL, z, xo, yo=NULL, linear = TRUE, extrap = FALSE, "interpp"<-function(x, y=NULL, z, xo, yo=NULL, linear = TRUE, extrap = FALSE,
duplicate = "error", dupfun = NULL, duplicate = "error", dupfun = NULL,
jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE) jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE,
remove = !linear)
{ {
## handle sp data, save coordinate and value names ## handle sp data, save coordinate and value names
...@@ -29,12 +30,6 @@ ...@@ -29,12 +30,6 @@
"or both x and xo have to be SpatialPointsDataFrames", "or both x and xo have to be SpatialPointsDataFrames",
"and z a name of a data column in x")) "and z a name of a data column in x"))
} }
## FIXME: drop old akima code
if(linear)
ret <- interpp.old(x,y,z,xo,yo,ncp=0,extrap=FALSE,
duplicate=duplicate,dupfun=dupfun)
else {
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z)))) if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed") stop("missing values and Infs not allowed")
if(is.null(xo)) if(is.null(xo))
...@@ -91,6 +86,11 @@ ...@@ -91,6 +86,11 @@
zo <- rep(0, np) zo <- rep(0, np)
miss <- !extrap #if not extrapolating use missing values miss <- !extrap #if not extrapolating use missing values
extrap <- rep(extrap, np) extrap <- rep(extrap, np)
## Defaults from Fortran code for triangle removal moved to R code,
## see SDTRTT:
hbrmn <- 0.1 # height to base ratio for thin triangles
nrrtt <- 5 # recursion depth for triangle removal, 0 disables
if(!remove) nrrtt <- 0
ans <- .Fortran("sdbi3p", ans <- .Fortran("sdbi3p",
md = as.integer(1), md = as.integer(1),
...@@ -103,23 +103,29 @@ ...@@ -103,23 +103,29 @@
y = as.double(yo), y = as.double(yo),
z = as.double(zo), z = as.double(zo),
ier = integer(1), ier = integer(1),
wk = double(17 * n), wk = double(17 * 5 * n),
iwk = integer(25 * n), iwk = integer(39 * n),
extrap = as.logical(extrap), extrap = as.logical(extrap),
near = integer(n),
nxt = integer(n),
dist = double(n),
linear = as.logical(linear), linear = as.logical(linear),
hbrmn = as.double(hbrmn),
nrrtt = as.integer(nrrtt),
PACKAGE = "akima") PACKAGE = "akima")
if(miss) if(miss)
ans$z[ans$extrap]<-NA ans$z[ans$extrap]<-NA
## Error code 10 from sdsf3p indicates error code -2 from trmesh:
## Error code 11 from sdbi3p indicates failure of thin triangle removal.
if(ans$ier==11){
stop("removal of thin triangles from border failed. try to re-run with remove=FALSE")
}
## Error code 10 from sdbi3p indicates error code -2 from trmesh:
## first three points collinear. ## first three points collinear.
## Try to add jitter to data locations to avoid collinearities, ## Try to add jitter to data locations to avoid collinearities,
## start with diff(range(x))*jitter*jitter.trials^1.5 and repeat for ## start with diff(range(x))*jitter*jitter.trials^1.5 and repeat for
## jitter.trials steps until success (ier=0) ## jitter.trials steps until success (ier=0)
if(ans$ier==10) if(ans$ier==10){
warning("collinear points, trying to add some jitter to avoid collinearities!") warning("collinear points, trying to add some jitter to avoid collinearities!")
jitter.trials <- 1 jitter.trials <- 1
success <- FALSE success <- FALSE
...@@ -143,13 +149,12 @@ ...@@ -143,13 +149,12 @@
y = as.double(yo), y = as.double(yo),
z = as.double(zo), z = as.double(zo),
ier = integer(1), ier = integer(1),
wk = double(17 * n), wk = double(17 * 5 * n),
iwk = integer(25 * n), iwk = integer(39 * n),
extrap = as.logical(extrap), extrap = as.logical(extrap),
near = integer(n),
nxt = integer(n),
dist = double(n),
linear = as.logical(linear), linear = as.logical(linear),
hbrmn = as.double(hbrmn),
nrrtt = as.integer(nrrtt),
PACKAGE = "akima") PACKAGE = "akima")
if(miss) if(miss)
ans$z[ans$extrap] <- NA ans$z[ans$extrap] <- NA
...@@ -160,6 +165,7 @@ ...@@ -160,6 +165,7 @@
if(success) if(success)
warning("success: collinearities reduced through jitter") warning("success: collinearities reduced through jitter")
jitter.trials <- jitter.trials+1 jitter.trials <- jitter.trials+1
}
} }
if(is.sp){ if(is.sp){
nona <- !is.na(ans$z) nona <- !is.na(ans$z)
...@@ -170,6 +176,5 @@ ...@@ -170,6 +176,5 @@
} else { } else {
ret <- list(x=ans$x,y=ans$y,z=ans$z) ret <- list(x=ans$x,y=ans$y,z=ans$z)
} }
} ## END FIXME
ret ret
} }
"interpp.new"<-function(x, y, z, xo, yo, extrap = FALSE, "interpp.new"<-function(x, y, z, xo, yo, extrap = FALSE,
duplicate = "error", dupfun = NULL) duplicate = "error", dupfun = NULL)
{ {
warning("interp.new() is deprecated, use interp()") warning("interpp.new() is deprecated, use interpp()")
interpp(x, y, z, interpp(x, y, z,
xo, xo,
yo, linear=FALSE, yo, linear=FALSE,
......
"interpp.old"<-function(x, y, z, xo, yo, ncp = 0, extrap = FALSE, "interpp.old"<-function(x, y, z, xo, yo, ncp = 0, extrap = FALSE,
duplicate = "error", dupfun = NULL) duplicate = "error", dupfun = NULL)
{ {
warning("interpp.old() is deprecated, future versions will only provide interpp()") warning("interpp.old() is deprecated, use interpp()")
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z)))) interpp(x, y, z,
stop("missing values and Infs not allowed") xo,
if(is.null(xo)) yo, linear=(ncp==0),
stop("xo missing") extrap=extrap, duplicate=duplicate, dupfun=dupfun)
if(is.null(yo))
stop("yo missing")
if(ncp>25){
ncp <- 25
cat("ncp too large, using ncp=25\n")
}
drx <- diff(range(x))
dry <- diff(range(y))
if(drx == 0 || dry == 0)
stop("all data collinear") # other cases caught in Fortran code
if(drx/dry > 10000 || drx/dry < 0.0001)
stop("scales of x and y are too dissimilar")
n <- length(x)
np <- length(xo)
if(length(yo)!=np)
stop("length of xo and yo differ")
if(length(y) != n || length(z) != n)
stop("Lengths of x, y, and z do not match")
xy <- paste(x, y, sep =",")
i <- match(xy, xy)
if(duplicate=="user" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x)
)
}
if(duplicate!="strip"){
z <- unlist(lapply(split(z,i), centre))
ord <- !duplicated(xy)
x <- x[ord]
y <- y[ord]
n <- length(x)
}
else{
ord <- (hist(i,plot=FALSE,freq=TRUE,breaks=seq(0.5,max(i)+0.5,1))$counts==1)
x <- x[ord]
y <- y[ord]
z <- z[ord]
n <- length(x)
}
}
else
if(any(duplicated(xy)))
stop("duplicate data points")
zo <- rep(0, np)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
misso <- seq(miss, np)
if(extrap & ncp == 0)
warning("Cannot extrapolate with linear option")
ans <- .Fortran("idbvip",
as.integer(1),
as.integer(ncp),
as.integer(n),
as.double(x),
as.double(y),
as.double(z),
as.integer(np),
x = as.double(xo),
y = as.double(yo),
z = zo,
integer((31 + ncp) * n + np),
double(8 * n),
misso = as.logical(misso),
PACKAGE = "akima")
temp <- ans[c("x", "y", "z", "misso")]
temp$z[temp$misso]<-NA
temp[c("x", "y", "z")]
} }
...@@ -12,7 +12,8 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx), ...@@ -12,7 +12,8 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
yo=seq(min(y), max(y), length = ny), yo=seq(min(y), max(y), length = ny),
linear = TRUE, extrap=FALSE, duplicate = "error", dupfun = NULL, linear = TRUE, extrap=FALSE, duplicate = "error", dupfun = NULL,
nx = 40, ny = 40, nx = 40, ny = 40,
jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE) jitter = 10^-12, jitter.iter = 6, jitter.random = FALSE,
remove = !linear)
} }
\arguments{ \arguments{
\item{x}{ \item{x}{
...@@ -96,6 +97,12 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx), ...@@ -96,6 +97,12 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
\item{jitter.random}{logical, see \code{jitter}, defaults to \item{jitter.random}{logical, see \code{jitter}, defaults to
\code{FALSE} \code{FALSE}
} }
\item{remove}{logical, indicates whether Akimas removal of thin triangles along
the border of the convex hull should be performed, experimental setting!
defaults to \code{!linear}, so it will be left out for linear interpolation
by default. For some point configurations it is the only
available option to skip this removal step.
}
} }
\value{ \value{
list with 3 components: list with 3 components:
...@@ -113,12 +120,14 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx), ...@@ -113,12 +120,14 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
\code{SpatialPixelssDataFrame} is returned. \code{SpatialPixelssDataFrame} is returned.
} }
\note{ \note{
\code{interp} uses Akimas new Fortran code from 1996 \code{interp} uses Akimas new Fortran code (AMS 761) from 1996 in the revised
for spline interpolation, the triangulation (based on Renkas tripack) version by Renka from 1998 for spline interpolation, the triangulation
is reused for linear interpolation. In this newer version Akima (based on Renkas tripack) is reused for linear interpolation. In this
switched from his own triangulation to Renkas tripack (=TOMS 751). newer version Akima switched from his own triangulation to Renkas
tripack (AMS 751).
Note that S-Plus uses (used?) the old Fortran code from Akima 1978. Note that earlier versions (up to version 0.5-12) as well as S-Plus
used the old Fortran code from Akima 1978 (AMS 526).
The resulting structure is suitable for input to the The resulting structure is suitable for input to the
functions \code{\link{contour}} and \code{\link{image}}. Check functions \code{\link{contour}} and \code{\link{image}}. Check
...@@ -157,6 +166,10 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx), ...@@ -157,6 +166,10 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
two-dimensional Delaunay triangulation package. two-dimensional Delaunay triangulation package.
ACM Transactions on Mathematical Software. ACM Transactions on Mathematical Software.
\bold{22}, 1-8. \bold{22}, 1-8.