Commit 59c69498 authored by agebhard's avatar agebhard
Browse files

add bicubic interpolation for rectangular gridded data (TOMS 760)

parent 295100a3
Package: akima
Version: 0.5-9
Version: 0.5-10
Date: 2013-01-19
Title: Interpolation of irregularly spaced data
Author: Fortran code by H. Akima
......
......@@ -9,7 +9,7 @@ Copyrighted and Licensed by ACM,
see http://www.acm.org/publications/policies/softwarecrnotice
2. R interface (src/init.c R/*.R man/*.Rd data/*):
2. R interface (src/init.c src/akima.h R/*.R man/*.Rd data/*):
The R interface code has been developed as work based on the
ACM licensed code, hence it is also ACM licensed, copyright
......
useDynLib(akima)
export(aspline,interp,interpp,interp.old,interp.new,interp2xyz)
export(aspline,interp,interpp,interp.old,interp.new,interp2xyz,bicubic,bicubic.grid)
# R interface to Akimas regular grid spline interpolator
# (TOMS 464) improved version: TOMS 760
# http://www.netlib.org/toms/760
bicubic <- function(x,y,z,x0,y0){
nx <- length(x)
ny <- length(y)
if(dim(z)[1]!=nx)
stop("dim(z)[1] and length of x differs!")
if(dim(z)[2]!=ny)
stop("dim(z)[2] and length of y differs!")
n0 <- length(x0)
if(length(y0)!=n0)
stop("length of y0 and x0 differs!")
ret <- .Fortran("rgbi3p",
md=as.integer(1),
nxd=as.integer(nx),
nyd=as.integer(ny),
xd=as.double(x),
yd=as.double(y),
zd=as.double(z),
nip=as.integer(n0),
xi=as.double(x0),
yi=as.double(y0),
zi=double(n0),
ier=integer(1),
wk=double(3*nx*ny),
PACKAGE="akima")
list(x=x0,y=y0,z=ret$zi)
}
bicubic.grid <- function(x,y,z,xlim,ylim,dx,dy){
nx <- length(x)
ny <- length(y)
if(dim(z)[1]!=nx)
stop("dim(z)[1] and length of x differs!")
if(dim(z)[2]!=ny)
stop("dim(z)[2] and length of y differs!")
xi <- seq(xlim[1],xlim[2],by=dx)
nx <- length(xi)
yi <- seq(ylim[1],ylim[2],by=dy)
ny <- length(yi)
xmat <- matrix(rep(xi,ny),nrow=ny,ncol=nx,byrow=TRUE)
ymat <- matrix(rep(yi,nx),nrow=ny,ncol=nx,byrow=FALSE)
xy <- cbind(c(xmat),c(ymat))
n0 <- nx*ny
ret <- bicubic(x,y,z,xy[,1],xy[,2])
# return cell boundaries
list(x=xi,y=yi,z=t(matrix(ret$z,nrow=ny,ncol=nx,byrow=F)))
}
\name{bicubic}
\alias{bicubic}
\title{
Bivariate Interpolation for Data on a Rectangular grid
}
\description{
The description in the Fortran code says:
This subroutine performs interpolation of a bivariate function,
z(x,y), on a rectangular grid in the x-y plane. It is based on
the revised Akima method.
In this subroutine, the interpolating function is a piecewise
function composed of a set of bicubic (bivariate third-degree)
polynomials, each applicable to a rectangle of the input grid
in the x-y plane. Each polynomial is determined locally.
This subroutine has the accuracy of a bicubic polynomial, i.e.,
it interpolates accurately when all data points lie on a
surface of a bicubic polynomial.
The grid lines can be unevenly spaced.
}
\usage{
bicubic(x, y, z, x0, y0)
}
\arguments{
\item{x}{
a vector containing the \code{x} coordinates of the rectangular data grid.
}
\item{y}{
a vector containing the \code{y} coordinates of the rectangular data grid.
}
\item{z}{
a matrix containing the \code{z[i,j]} data values for the grid points (\code{x[i]},\code{y[j]}).
}
\item{x0}{
vector of \code{x} coordinates used to interpolate at.
}
\item{y0}{
vector of \code{y} coordinates used to interpolate at.
}
}
\details{
This functiuon is a R interface to Akima's Rectangular-Grid-Data
Fitting algorithm (TOMS 760). The algorithm has the accuracy of a bicubic
(bivariate third-degree) polynomial.
}
\value{
This function produces a list of interpolated points:
\item{x}{vector of \code{x} coordinates.}
\item{y}{vector of \code{y} coordinates.}
\item{z}{vector of interpolated data \code{z}.}
If you need an output grid, see \code{\link{bicubic.grid}}.
}
\references{
Akima, H. (1996) Rectangular-Grid-Data
Surface Fitting that Has the Accuracy of a
Bicubic Polynomial,
J. ACM \bold{22}(3), 357-361
}
\note{
Use \code{\link{interp}} for the general case of irregular gridded data!
}
\seealso{
\code{\link{interp}}, \code{\link{bicubic.grid}}
% maybe later:
% \code{\link[rgeostat]{bilinear}}
}
\examples{
data(akima760)
# interpolate at the diagonal of the grid [0,8]x[0,10]
akima.bic <- bicubic(akima760$x,akima760$y,akima760$z,
seq(0,8,length=50), seq(0,10,length=50))
plot(sqrt(akima.bic$x^2+akima.bic$y^2), akima.bic$z, type="l")
##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (x, y, z, x0, y0)
{
nx <- length(x)
ny <- length(y)
if (dim(z)[1] != nx)
stop("dim(z)[1] and length of x differs!")
if (dim(z)[2] != ny)
stop("dim(z)[2] and length of y differs!")
n0 <- length(x0)
if (length(y0) != n0)
stop("length of y0 and x0 differs!")
ret <- .Fortran("rgbi3p", md = as.integer(1), nxd = as.integer(nx),
nyd = as.integer(ny), xd = as.double(x), yd = as.double(y),
zd = as.double(z), nip = as.integer(n0), xi = as.double(x0),
yi = as.double(y0), zi = double(n0), ier = integer(1),
wk = double(3 * nx * ny), PACKAGE = "akima")
list(x = x0, y = y0, z = ret$zi)
}
}
\keyword{ dplot }
\name{bicubic.grid}
\alias{bicubic.grid}
\title{
Bivariate Interpolation for Data on a Rectangular grid
}
\description{
The description in the Fortran code says:
This subroutine performs interpolation of a bivariate function,
z(x,y), on a rectangular grid in the x-y plane. It is based on
the revised Akima method.
In this subroutine, the interpolating function is a piecewise
function composed of a set of bicubic (bivariate third-degree)
polynomials, each applicable to a rectangle of the input grid
in the x-y plane. Each polynomial is determined locally.
This subroutine has the accuracy of a bicubic polynomial, i.e.,
it interpolates accurately when all data points lie on a
surface of a bicubic polynomial.
The grid lines can be unevenly spaced.
}
\usage{
bicubic.grid(x,y,z,xlim,ylim,dx,dy)
}
\arguments{
\item{x}{
a vector containing the \code{x} coordinates of the rectangular data grid.
}
\item{y}{
a vector containing the \code{y} coordinates of the rectangular data grid.
}
\item{z}{
a matrix containing the \code{z[i,j]} data values for the grid points (\code{x[i]},\code{y[j]}).
}
\item{xlim}{
vector of length 2 giving lower and upper limit for range \code{x}
coordinates used for output grid.
}
\item{ylim}{
vector of length 2 giving lower and upper limit for range of \code{y}
coordinates used for output grid.
}
\item{dx}{
output grid spacing in \code{x} direction.
}
\item{dy}{
output grid spacing in \code{y} direction.
}
}
\details{
This functiuon is a R interface to Akima's Rectangular-Grid-Data
Fitting algorithm (TOMS 760). The algorithm has the accuracy of a bicubic
(bivariate third-degree) polynomial.
}
\value{
This function produces a grid of interpolated points, feasible to be
used directly with \code{\link{image}} and \code{\link{contour}}:
\item{x}{vector of \code{x} coordinates of the output grid.}
\item{y}{vector of \code{y} coordinates of the output grid.}
\item{z}{matrix of interpolated data for the output grid.}
}
\references{
Akima, H. (1996) Rectangular-Grid-Data
Surface Fitting that Has the Accuracy of a
Bicubic Polynomial,
J. ACM \bold{22}(3), 357-361
}
\note{
Use \code{\link{interp}} for the general case of irregular gridded data!
}
\seealso{
\code{\link{interp}}, \code{\link{bicubic}}
% maybe later:
% \code{\link[rgeostat]{bilinear}}
}
\examples{
data(akima760)
# interpolate at a grid [0,8]x[0,10]
akima.bic <- bicubic.grid(akima760$x,akima760$y,akima760$z,
c(0,8),c(0,10),0.1,0.1)
image(akima.bic)
contour(akima.bic, add=TRUE)
}
\keyword{ dplot }
......@@ -28,3 +28,5 @@ int F77_NAME(uvip3p) (int *np, int *nd, double *xd, double *yd,
int F77_NAME(intrpl) (int *l,double *x, double *y, int *n,
double *u, double *v, int *err);
int F77_NAME(rgbi3p) (int *md, int *nxd, int *nyd, double *xd, double *yd, double *zd,
int *nip, double *xi, double *yi, double *zi, int *err);
......@@ -24,6 +24,20 @@
/* Fortran interface descriptions: */
static R_NativePrimitiveArgType rgbi3p_t[12] = {
INTSXP, /* MD, */
INTSXP, /* NXP, */
INTSXP, /* NYP, */
REALSXP, /* XD, */
REALSXP, /* YD, */
REALSXP, /* ZD, */
INTSXP, /* NIP, */
REALSXP, /* XI, */
REALSXP, /* YI, */
REALSXP, /* ZI, */
INTSXP, /* IER, */
REALSXP /* WK, */
};
static R_NativePrimitiveArgType idbvip_t[13] = {
INTSXP, /* MD, */
INTSXP, /* NCP, */
......@@ -124,6 +138,7 @@ static R_FortranMethodDef fortranMethods[] = {
{"sdsf3p", (DL_FUNC) &F77_SUB(sdsf3p), 17, sdsf3p_t}, /* interp.new */
{"uvip3p", (DL_FUNC) &F77_SUB(uvip3p), 8, uvip3p_t}, /* aspline */
{"intrpl", (DL_FUNC) &F77_SUB(intrpl), 7, intrpl_t}, /* aspline */
{"rgbi3p", (DL_FUNC) &F77_SUB(rgbi3p), 12, sdbi3p_t}, /* bicubic */
{NULL, NULL, 0}
};
......
c http://gams.nist.gov/serve.cgi/ModuleComponent/13334/Fullsource/NETLIB/760
SUBROUTINE RGBI3P(MD,NXD,NYD,XD,YD,ZD,NIP,XI,YI, ZI,IER, WK)
*
* Rectangular-grid bivariate interpolation
* (a master subroutine of the RGBI3P/RGSF3P subroutine package)
*
* Hiroshi Akima
* U.S. Department of Commerce, NTIA/ITS
* Version of 1995/08
*
* This subroutine performs interpolation of a bivariate function,
* z(x,y), on a rectangular grid in the x-y plane. It is based on
* the revised Akima method.
*
* In this subroutine, the interpolating function is a piecewise
* function composed of a set of bicubic (bivariate third-degree)
* polynomials, each applicable to a rectangle of the input grid
* in the x-y plane. Each polynomial is determined locally.
*
* This subroutine has the accuracy of a bicubic polynomial, i.e.,
* it interpolates accurately when all data points lie on a
* surface of a bicubic polynomial.
*
* The grid lines can be unevenly spaced.
*
* The input arguments are
* MD = mode of computation
* = 1 for new XD, YD, or ZD data (default)
* = 2 for old XD, YD, and ZD data,
* NXD = number of the input-grid data points in the x
* coordinate (must be 2 or greater),
* NYD = number of the input-grid data points in the y
* coordinate (must be 2 or greater),
* XD = array of dimension NXD containing the x coordinates
* of the input-grid data points (must be in a
* monotonic increasing order),
* YD = array of dimension NYD containing the y coordinates
* of the input-grid data points (must be in a
* monotonic increasing order),
* ZD = two-dimensional array of dimension NXD*NYD
* containing the z(x,y) values at the input-grid data
* points,
* NIP = number of the output points at which interpolation
* of the z value is desired (must be 1 or greater),
* XI = array of dimension NIP containing the x coordinates
* of the output points,
* YI = array of dimension NIP containing the y coordinates
* of the output points.
*
* The output arguments are
* ZI = array of dimension NIP where the interpolated z
* values at the output points are to be stored,
* IER = error flag
* = 0 for no errors
* = 1 for NXD = 1 or less
* = 2 for NYD = 1 or less
* = 3 for identical XD values or
* XD values out of sequence
* = 4 for identical YD values or
* YD values out of sequence
* = 5 for NIP = 0 or less.
*
* The other argument is
* WK = three dimensional array of dimension 3*NXD*NYD used
* internally as a work area.
*
* The very fisrt call to this subroutine and the call with a new
* XD, YD, and ZD array must be made with MD=1. The call with MD=2
* must be preceded by another call with the same XD, YD, and ZD
* arrays. Between the call with MD=2 and its preceding call, the
* WK array must not be disturbed.
*
* The constant in the PARAMETER statement below is
* NIPIMX = maximum number of output points to be processed
* at a time.
* The constant value has been selected empirically.
*
* This subroutine calls the RGPD3P, RGLCTN, and RGPLNL subroutines.
*
*
* Specification statements
* .. Parameters ..
INTEGER NIPIMX
PARAMETER (NIPIMX=51)
* ..
* .. Scalar Arguments ..
INTEGER IER,MD,NIP,NXD,NYD
* ..
* .. Array Arguments ..
REAL*8 WK(3,NXD,NYD),XD(NXD),XI(NIP),YD(NYD),YI(NIP),
+ ZD(NXD,NYD),ZI(NIP)
* ..
* .. Local Scalars ..
INTEGER IIP,IX,IY,NIPI
* ..
* .. Local Arrays ..
INTEGER INXI(NIPIMX),INYI(NIPIMX)
* ..
* .. External Subroutines ..
EXTERNAL RGLCTN,RGPD3P,RGPLNL
* ..
* .. Intrinsic Functions ..
INTRINSIC MIN
* ..
* Preliminary processing
* Error check
IF (NXD.LE.1) GO TO 40
IF (NYD.LE.1) GO TO 50
DO 10 IX = 2,NXD
IF (XD(IX).LE.XD(IX-1)) GO TO 60
10 CONTINUE
DO 20 IY = 2,NYD
IF (YD(IY).LE.YD(IY-1)) GO TO 70
20 CONTINUE
IF (NIP.LE.0) GO TO 80
IER = 0
* Calculation
* Estimates partial derivatives at all input-grid data points
* (for MD=1).
IF (MD.NE.2) THEN
CALL RGPD3P(NXD,NYD,XD,YD,ZD, WK)
* CALL RGPD3P(NXD,NYD,XD,YD,ZD, PDD)
END IF
* DO-loop with respect to the output point
* Processes NIPIMX output points, at most, at a time.
DO 30 IIP = 1,NIP,NIPIMX
NIPI = MIN(NIP- (IIP-1),NIPIMX)
* Locates the output points.
CALL RGLCTN(NXD,NYD,XD,YD,NIPI,XI(IIP),YI(IIP), INXI,INYI)
* CALL RGLCTN(NXD,NYD,XD,YD,NIP,XI,YI, INXI,INYI)
* Calculates the z values at the output points.
CALL RGPLNL(NXD,NYD,XD,YD,ZD,WK,NIPI,XI(IIP),YI(IIP),INXI,
+ INYI, ZI(IIP))
* CALL RGPLNL(NXD,NYD,XD,YD,ZD,PDD,NIP,XI,YI,INXI,INYI, ZI)
30 CONTINUE
RETURN
* Error exit
40 CONTINUE
C WRITE (*,FMT=9000)
IER = 1
GO TO 90
50 CONTINUE
C WRITE (*,FMT=9010)
IER = 2
GO TO 90
60 CONTINUE
C WRITE (*,FMT=9020) IX,XD(IX)
IER = 3
GO TO 90
70 CONTINUE
C WRITE (*,FMT=9030) IY,YD(IY)
IER = 4
GO TO 90
80 CONTINUE
C WRITE (*,FMT=9040)
IER = 5
90 CONTINUE
C WRITE (*,FMT=9050) NXD,NYD,NIP
RETURN
* Format statements for error messages
9000 FORMAT (1X,/,'*** RGBI3P Error 1: NXD = 1 or less')
9010 FORMAT (1X,/,'*** RGBI3P Error 2: NYD = 1 or less')
9020 FORMAT (1X,/,'*** RGBI3P Error 3: Identical XD values or',
+ ' XD values out of sequence',/,' IX =',I6,', XD(IX) =',
+ E11.3)
9030 FORMAT (1X,/,'*** RGBI3P Error 4: Identical YD values or',
+ ' YD values out of sequence',/,' IY =',I6,', YD(IY) =',
+ E11.3)
9040 FORMAT (1X,/,'*** RGBI3P Error 5: NIP = 0 or less')
9050 FORMAT (' NXD =',I5,', NYD =',I5,', NIP =',I5,/)
END
SUBROUTINE RGSF3P(MD,NXD,NYD,XD,YD,ZD,NXI,XI,NYI,YI, ZI,IER, WK)
*
* Rectangular-grid surface fitting
* (a master subroutine of the RGBI3P/RGSF3P subroutine package)
*
* Hiroshi Akima
* U.S. Department of Commerce, NTIA/ITS
* Version of 1995/08
*
* This subroutine performs surface fitting by interpolating
* values of a bivariate function, z(x,y), on a rectangular grid
* in the x-y plane. It is based on the revised Akima method.
*
* In this subroutine, the interpolating function is a piecewise
* function composed of a set of bicubic (bivariate third-degree)
* polynomials, each applicable to a rectangle of the input grid
* in the x-y plane. Each polynomial is determined locally.
*
* This subroutine has the accuracy of a bicubic polynomial, i.e.,
* it fits the surface accurately when all data points lie on a
* surface of a bicubic polynomial.
*
* The grid lines of the input and output data can be unevenly
* spaced.
*
* The input arguments are
* MD = mode of computation
* = 1 for new XD, YD, or ZD data (default)
* = 2 for old XD, YD, and ZD data,
* NXD = number of the input-grid data points in the x
* coordinate (must be 2 or greater),
* NYD = number of the input-grid data points in the y
* coordinate (must be 2 or greater),
* XD = array of dimension NXD containing the x coordinates
* of the input-grid data points (must be in a
* monotonic increasing order),
* YD = array of dimension NYD containing the y coordinates
* of the input-grid data points (must be in a
* monotonic increasing order),
* ZD = two-dimensional array of dimension NXD*NYD
* containing the z(x,y) values at the input-grid data
* points,
* NXI = number of output grid points in the x coordinate
* (must be 1 or greater),
* XI = array of dimension NXI containing the x coordinates
* of the output grid points,
* NYI = number of output grid points in the y coordinate
* (must be 1 or greater),
* YI = array of dimension NYI containing the y coordinates
* of the output grid points.
*
* The output arguments are
* ZI = two-dimensional array of dimension NXI*NYI, where
* the interpolated z values at the output grid
* points are to be stored,
* IER = error flag
* = 0 for no error
* = 1 for NXD = 1 or less
* = 2 for NYD = 1 or less
* = 3 for identical XD values or
* XD values out of sequence
* = 4 for identical YD values or
* YD values out of sequence
* = 5 for NXI = 0 or less
* = 6 for NYI = 0 or less.
*
* The other argument is
* WK = three-dimensional array of dimension 3*NXD*NYD used
* internally as a work area.
*
* The very first call to this subroutine and the call with a new
* XD, YD, or ZD array must be made with MD=1. The call with MD=2
* must be preceded by another call with the same XD, YD, and ZD
* arrays. Between the call with MD=2 and its preceding call, the
* WK array must not be disturbed.
*
* The constant in the PARAMETER statement below is
* NIPIMX = maximum number of output points to be processed
* at a time.
* The constant value has been selected empirically.
*
* This subroutine calls the RGPD3P, RGLCTN, and RGPLNL subroutines.
*
*
* Specification statements
* .. Parameters ..
INTEGER NIPIMX
PARAMETER (NIPIMX=51)
* ..
* .. Scalar Arguments ..
INTEGER IER,MD,NXD,NXI,NYD,NYI
* ..
* .. Array Arguments ..
REAL*8 WK(3,NXD,NYD),XD(NXD),XI(NXI),YD(NYD),YI(NYI),
+ ZD(NXD,NYD),ZI(NXI,NYI)
* ..
* .. Local Scalars ..
INTEGER IX,IXI,IY,IYI,NIPI
* ..
* .. Local Arrays ..
REAL*8 YII(NIPIMX)
INTEGER INXI(NIPIMX),INYI(NIPIMX)
* ..
* .. External Subroutines ..
EXTERNAL RGLCTN,RGPD3P,RGPLNL