* Fix ACM license URL

* remove old Akima code and interface, everything is now handled by
  Akima new code
* Added Frankes interpolation example data (1979, see acording man page)
* Added jitter option to interp to allow for collinear data, see examples

.Rbuildignore

@@ -2,3 +2,7 @@ orig
 config.log 
 config.status 
 autom4te.cache
+src-c
+src-i386
+build-debug.sh 
+work

DESCRIPTION

 Package: akima
 Version: 0.6-1
-Date: 2016-12-09
+Date: 2016-12-16
 Title: Interpolation of Irregularly and Regularly Spaced Data
 Authors@R: c(person("Hiroshi", "Akima", role = c("aut", "cph"),
 		     comment = "Fortran code (TOMS 526, 761, 697 and 433)"),
@@ -14,6 +14,7 @@ Authors@R: c(person("Hiroshi", "Akima", role = c("aut", "cph"),
               person("YYYY Association for Computing Machinery, Inc.",
                       role = c("cph"), comment = "covers code from TOMS 526, 761, 697 and 433")
             )
+Maintainer: Albrecht Gebhardt <albrecht.gebhardt@aau.at>
 Description: Several cubic spline interpolation methods of H. Akima for irregular and
   regular gridded data are available through this package, both for the bivariate case
   (irregular data: ACM 526 and ACM 761, regular data: ACM 433) and univariate case (ACM 697).

NAMESPACE

 useDynLib(akima)
 
-export(aspline,interp,interpp,interp2xyz,bicubic,bicubic.grid,bilinear,bilinear.grid)
+export(aspline,interp,interpp,interp2xyz,bicubic,bicubic.grid,bilinear,bilinear.grid,franke.fn,franke.data)
 
 importFrom("sp","coordinates")
 importFrom("sp","coordinates<-")

R/franke.R

+franke.fn <- function(x,y,fn=1){
+    switch(fn,
+           "1"=0.75*exp(-((9*x-2)^2+(9*y-2)^2)/4)+
+               0.75*exp(-((9*x+1)^2)/49-(9*y+1)/10)+
+               0.5*exp(-((9*x-7)^2+(9*y-3)^2)/4)-
+               0.2*exp(-(9*x-4)^2-(9*y-7)^2),
+           "2"=(tanh(9*y-9*x)+1)/9,
+           "3"=(1.25+cos(5.4*y))/(6*(1+(3*x-1)^2)),
+           "4"=exp(-81*((x-0.5)^2+(y-0.5)^2)/16)/3,
+           "5"=exp(-81*((x-0.5)^2+(y-0.5)^2)/4)/3,
+           "6"=sqrt(64-81*((x-0.5)^2+(y-0.5)^2))/9-0.5)
+}
+
+franke.data <- function(fn=1,ds=1,data){
+    ret <- cbind(x=data[[ds]]$x,y=data[[ds]]$y,
+                 z=franke.fn(data[[ds]]$x,y=data[[ds]]$y,fn))
+    list(x=ret[,"x"],y=ret[,"y"],z=ret[,"z"])
+}

R/interp.R

@@ -106,10 +106,29 @@ interp <-
         jitter.trials <- 1
         success <- FALSE
         while(jitter.trials<jitter.iter & !success){
-            ## dont use random jitter for reproducabilty by default:
+            ## dont use random jitter for reproducabilty by default,
+            ## same as in tri.mesh in package tripack:
             if(jitter.random){
-                xj <- x+rnorm(n,0,diff(range(x))*jitter*jitter.trials^1.5)
-                yj <- y+rnorm(n,0,diff(range(y))*jitter*jitter.trials^1.5)
+                ## the randomness is only contained in a random shift
+                ## of a regular -1,+/-0,+1,-1,+/-0,+1,... scheme
+                ## determining the fixed amounts of jitter to be added
+                ## to the x and y values separately.
+                ## Really random jitter like this
+                ## xj <- x+rnorm(n,0,diff(range(x))*jitter*jitter.trials^1.5)
+                ## yj <- y+rnorm(n,0,diff(range(y))*jitter*jitter.trials^1.5)
+                ## still triggers spurious errors in the triangulation.
+                ## j <- sample(1:length(x), length(x))
+                ## xj <- x[j]
+                ## yj <- y[j]
+                ## but it needs afterwards back ordering!
+                j <- list()
+                j[[1]] <- rep(c(-1,0,1),length.out=length(x))
+                j[[2]] <- rep(c(0,1,-1),length.out=length(x))
+                j[[3]] <- rep(c(1,-1,0),length.out=length(x))
+                jx <- sample(1:3,1)
+                jy <- sample(1:3,1)
+                xj <- x+j[[jx]]*diff(range(x))*jitter*jitter.trials^1.5
+                yj <- y+j[[jy]]*diff(range(y))*jitter*jitter.trials^1.5
             } else {
                 xj <- x+rep(c(-1,0,1),length.out=length(x))*diff(range(x))*jitter*jitter.trials^1.5
                 yj <- y+rep(c(0,1,-1),length.out=length(y))*diff(range(y))*jitter*jitter.trials^1.5
@@ -117,8 +136,8 @@ interp <-
             ans <- .Fortran("sdsf3p",
                             as.integer(1),
                             as.integer(n),
-                            as.double(xj),
-                            as.double(yj),
+                            xd=as.double(xj),
+                            yd=as.double(yj),
                             as.double(z),
                             as.integer(nx),
                             x = as.double(xo),
@@ -136,6 +155,8 @@ interp <-
                             PACKAGE = "akima")
             if(miss)
                 ans$z[ans$extrap] <- NA
+            if(linear)
+                ans$z[ans$extrap] <- NA
             success <- (ans$ier==0)
             if(success)
                 warning("success: collinearities reduced through jitter")

R/interpp.R

@@ -84,7 +84,7 @@
             stop("duplicate data points")
     zo <- rep(0, np)
     miss <- !extrap	#if not extrapolating use missing values
-    extrap <- seq(extrap, np)
+    extrap <- rep(extrap, np)
 
     ans <- .Fortran("sdbi3p",
                     md = as.integer(1),
@@ -129,8 +129,8 @@
         ans <- .Fortran("sdbi3p",
                         md = as.integer(1),
                         ndp = as.integer(n),
-                        xd = as.double(x),
-                        yd = as.double(y),
+                        xd = as.double(xj),
+                        yd = as.double(yj),
                         zd = as.double(z),
                         nx = as.integer(np),
                         x = as.double(xo),
@@ -147,6 +147,9 @@
                         PACKAGE = "akima")
         if(miss)
             ans$z[ans$extrap] <- NA
+        if(linear)
+            ans$z[ans$extrap] <- NA
+
         success <- (ans$ier==0)
         if(success)
             warning("success: collinearities reduced through jitter")

man/franke.data.Rd

+\name{franke.data}
+\alias{franke.data}
+\alias{franke.fn}
+\alias{franke}
+\title{
+Test datasets from Franke for interpolation of scattered data
+}
+\description{
+\code{franke.data} generates the test datasets from Franke, 1979, see references.
+}
+\usage{
+franke.data(fn = 1, ds = 1, data)
+franke.fn(x, y, fn = 1)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{fn}{
+    function number, from 1 to 5.
+  }
+  \item{x}{'x' value}
+  \item{y}{'y' value}
+  \item{ds}{
+    data set number, from 1 to 3. Dataset 1 consists of 100 points,
+    dataset 2 of 33 points and dataset 3 of 25 points scattered in the
+    square \eqn{[0,1]\times[0,1]}{[0,1]x[0,1]}. (and partially slightly
+    outside).
+  }
+  \item{data}{
+    A list of dataframes with 'x' and 'y' to choose from, dataset
+    \code{franke} should be used here.
+  }
+}
+\details{
+  These datasets are mentioned in [Akima 1996] as testbed for the
+  irregular scattered data interpolator.
+
+  Franke used the five functions:
+
+  \deqn{0.75e^{-\frac{(9x-2)^2+(9y-2)^2}{4}}+
+    0.75e^{-\frac{(9x+1)^2}{49}-\frac{9y+1}{10}}+
+    0.5e^{-\frac{(9x-7)^2+(9y-3)^2}{4}}-
+      0.2e^{-((9x-4)^2-(9y-7)^2)}
+    }{0.75*exp(-((9*x-2)^2+(9*y-2)^2)/4)+
+      0.75*exp(-((9*x+1)^2)/49-(9*y+1)/10)+
+      0.5*exp(-((9*x-7)^2+(9*y-3)^2)/4)-
+      0.2*exp(-(9*x-4)^2-(9*y-7)^2)}
+
+
+  \deqn{\frac{\mbox{tanh}(9y-9x)+1}{9}}{(tanh(9*y-9*x)+1)/9}
+
+  \deqn{\frac{1.25+\cos(5.4y)}{6(1+(3x-1)^2)}}{(1.25+cos(5.4*y))/(6*(1+(3*x-1)^2))}
+
+  \deqn{e^{-\frac{81((x-0.5)^2+\frac{(y-0.5)^2}{16})}{3}}}{exp(-81*((x-0.5)^2+(y-0.5)^2)/16)/3}
+
+  \deqn{e^{-\frac{81((x-0.5)^2+\frac{(y-0.5)^2}{4})}{3}}}{exp(-81*((x-0.5)^2+(y-0.5)^2)/4)/3}
+
+  \deqn{\frac{\sqrt{64-81((x-0.5)^2+(y-0.5)^2)}}{9}-0.5}{sqrt(64-81*((x-0.5)^2+(y-0.5)^2))/9-0.5}
+
+  and evaluated them on different more or less dense grids over \eqn{[0,1]\times[0,1]}{[0,1]x[0,1]}.
+}
+\value{
+A data frame with components
+  \item{x }{'x' coordinate}
+  \item{y }{'y' coordinate}
+  \item{z }{'z' value}
+}
+\note{
+  The datasets have to be generated via \code{franke.data} before
+  use, the dataset \code{franke} only contains a list of 3 dataframes of
+  'x' and 'y' coordinates for the above mentioned irregular grids.
+  Do not forget to load the \code{franke} dataset first.
+
+  The 'x' and 'y' values have been taken out of [Akima 1996].
+}
+\references{
+FRANKE, R., (1979). A critical comparison of some methods for interpolation
+of scattered data. Tech. Rep. NPS-53-79-003, Dept. of Mathematics, Naval
+Postgraduate School, Monterey, Calif.
+
+Akima, H. (1996). Algorithm 761: scattered-data surface fitting that has
+the accuracy of a cubic polynomial.
+ACM Transactions on Mathematical Software \bold{22}, 362--371.
+}
+\author{
+A. Gebhardt,
+}
+
+\seealso{
+\code{\link{interp}}
+}
+\examples{
+## generate Frankes data set for function 2 and dataset 3:
+data(franke)
+F23 <- franke.data(2,3,franke)
+str(F23)
+}
+\keyword{ datagen }
+

man/interp.Rd

@@ -92,7 +92,7 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
     same triangulation as used for interpolation, see examples below.
   }
   \item{jitter.iter}{number of iterations to retry with jitter, amount
-    will be increased in each iteration by \code{iter^1.5}}
+    will be multiplied in each iteration by \code{iter^1.5}}
   \item{jitter.random}{logical, see \code{jitter}, defaults to
     \code{FALSE}
   }
@@ -113,9 +113,12 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
   \code{SpatialPixelssDataFrame} is returned.
 }
 \note{
-  \code{interp} uses the same Fortran code from Akima 1978 as the S-Plus
-  version for linear interpolation and his newer Fortran code from 1996
-  for spline interpolation.
+  \code{interp} uses Akimas new Fortran code from 1996
+  for spline interpolation, the triangulation (based on Renkas tripack)
+  is reused for linear interpolation. In this newer version Akima
+  switched from his own triangulation to Renkas tripack (=TOMS 751).
+
+  Note that S-Plus uses (used?) the old Fortran code from Akima 1978.
 
   The resulting structure is suitable for input to the
   functions \code{\link{contour}} and \code{\link{image}}.  Check
@@ -149,12 +152,18 @@ interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
   Akima, H. (1996). Algorithm 761: scattered-data surface fitting that has
   the accuracy of a cubic polynomial.
   ACM Transactions on Mathematical Software \bold{22}, 362--371.
+
+  R. J. Renka (1996). Algorithm 751: TRIPACK: a constrained
+  two-dimensional Delaunay triangulation package.
+  ACM Transactions on Mathematical Software.
+  \bold{22}, 1-8.
 }
 \seealso{
   \code{\link{contour}}, \code{\link{image}},
   \code{\link{approx}}, \code{\link{spline}},
   \code{\link{aspline}},
-  \code{\link{outer}}, \code{\link{expand.grid}}.
+  \code{\link{outer}}, \code{\link{expand.grid}},
+  \code{link{franke.data}}.
 }
 \examples{
 data(akima)
@@ -217,11 +226,6 @@ filled.contour(akima.spl, color.palette = full.pal,
                       points(akima, pch = 3, col= hcl(c=100, l = 20))})
 ## no extrapolation!
 
-## example with duplicate points :
-
-### ... TODO
-
-
 \dontrun{
     ## interp can handle spatial point dataframes created by the sp package:
     library(sp)
@@ -238,6 +242,7 @@ filled.contour(akima.spl, color.palette = full.pal,
 ### An example demonstrating the problems that occur for rectangular
 ### gridded data.
 ###
+require(tripack)
 ### Create irregularly spaced sample data on even values of x and y
 ### (the "14" makes it irregular spacing).
 x   <- c(seq(2,10,2),14)
@@ -271,7 +276,8 @@ title(main="bicubic interpolation",
       xlab="bcubic.grid(...)",
       sub="Akimas regular grid version, ACM 760")
 
-### Now Akima bicubic splines for irregular gridded data:
+### Now Akima splines with accurracy of bicubic polynomial
+### for irregular gridded data:
 iRspl <- with(df,interp(x,y,z,linear=FALSE,nx=250,ny=250))
 ### Note that the triangulation is created by adding small amounts
 ### of jitter to the coordinates, resulting in an unique triangulation.
@@ -281,8 +287,9 @@ iRspl <- with(df,interp(x,y,z,linear=FALSE,nx=250,ny=250))
 image(iRspl,breaks=breaks,col = colors)
 contour(iRspl,col="black",levels=breaks,add=TRUE)
 plot(tri.mesh(xy$x,xy$y),col="white",add=TRUE)
-title(main="bicubic interpolation",
+title(main="bicubic* interpolation",
       xlab="interp(...,linear=FALSE)",
+      ylab="*: accuracy of bicubic polynomial"
       sub="Akimas irregular grid version, ACM 761")
 
 ### Just for comparison an implementation of bilinear interpolation,
@@ -311,12 +318,10 @@ plot(tri.mesh(xy$x,xy$y),col="white",add=TRUE)
 title(main="linear interpolation",
       xlab="interp(...,linear=TRUE)",
       sub="same triangulation as Akima irregular grid")
+
 ### And now four times Akima 761 with random jitter for
 ### triangulation correction, note that now interp() and tri.mesh()
-### need the same random seed to produce identical triangualtions!
-### Note that even identical (on the first look) triangulations deliver
-### different interpolations (as the jittered coordinates now differ by
-### some small amount).
+### need the same random seed to produce identical triangulations!
 for(i in 1:4){
     set.seed(42+i)
     iRspl <- with(df,interp(x,y,z,linear=FALSE,nx=250,ny=250,jitter.random=TRUE))
@@ -324,7 +329,7 @@ for(i in 1:4){
     contour(iRspl,col="black",levels=breaks,add=TRUE)
     set.seed(42+i)
     plot(tri.mesh(xy$x,xy$y,jitter.random=TRUE),col="white",add=TRUE)
-    title(main="bicubic interpolation",
+    title(main="bicubic* interpolation",
           xlab="interp(...,linear=FALSE)",
           ylab="random jitter added",
           sub="Akimas irregular grid version, ACM 761")
@@ -332,6 +337,16 @@ for(i in 1:4){
 par(old.par)
 
 }
+### Use all datasets from Franke, 1979:
+data(franke)
+for(i in 1:5)
+    for(j in 1:3){
+        FR <- franke.data(i,j,franke)
+        IL <- with(FR, interp(x,y,z,linear=FALSE))
+        image(IL)
+        contour(IL,add=TRUE)
+        with(FR,points(x,y))
+    }
 }
 \keyword{dplot}
 

man/interpp.Rd

@@ -136,6 +136,11 @@ interpp(x, y=NULL, z, xo, yo=NULL, linear=TRUE, extrap=FALSE,
   the accuracy of a cubic polynomial.
   ACM Transactions on Mathematical Software,
   \bold{22}, 362-371.
+
+  R. J. Renka (1996). Algorithm 751: TRIPACK: a constrained
+  two-dimensional Delaunay triangulation package.
+  ACM Transactions on Mathematical Software.
+  \bold{22}, 1-8.
 }
 \seealso{
   \code{\link[graphics]{contour}}, \code{\link[graphics]{image}},

src/akima433.f

@@ -4,6 +4,7 @@ C     ALGORITHM 433 COLLECTED ALGORITHMS FROM ACM.
 C     ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 10,
 C     P. 914.
       SUBROUTINE  INTRPL(L,X,Y,N,U,V,ERR)
+      IMPLICIT INTEGER (I-N), DOUBLE PRECISION (A-H,O-Z)
 C INTERPOLATION OF A SINGLE-VALUED FUNCTION
 C THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION
 C GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE
@@ -46,7 +47,9 @@ C PRELIMINARY PROCESSING
       IF(LM2.LT.0)        GO TO 90
       IF(N0.LE.0)         GO TO 91
       DO 11  I=2,L0
-        IF(X(I-1)-X(I))   11,95,96
+        IF(X(I-1)-X(I).LT.0.0D0) GO TO 11
+        IF(X(I-1)-X(I).EQ.0.0D0) GO TO 95
+        IF(X(I-1)-X(I).GT.0.0D0) GO TO 96
    11   CONTINUE
       IPV=0
 C MAIN DO-LOOP
@@ -152,7 +155,8 @@ C DETERMINATION OF THE COEFFICIENTS
         Q3=(-M3-M3+T3+T4)/(A3*A3)
 C COMPUTATION OF THE POLYNOMIAL
    70   DX=UK-P0
-   80   V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3))
+        V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3))
+   80 CONTINUE
       RETURN
 C ERROR EXIT
    90 ERR = 1

src/akima697.f

@@ -5,6 +5,7 @@ C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
 C      VOL. 17, NO. 3, SEPTEMBER, 1991, PP. 367.
 C
       SUBROUTINE  UVIP3P(NP,ND,XD,YD,NI,XI, YI, ERR)
+      IMPLICIT INTEGER (I-N), DOUBLE PRECISION (A-H,O-Z)
 C
 C Univariate Interpolation (Improved Akima Method)
 C

src/akima761.f

@@ -41,7 +41,7 @@
 *       = 2 for NDP not equal to NDPPV
 *       = 3 for NIP = 0 or less
 *       = 9 for errors in SDTRAN called by this subroutine, except the next one
-*       =10 for error 2 in SDTRAN (first three points collinear), as this can 
+*       =10 for error 2 in SDTRAN (first three points collinear), as this can
 *           be fixed by adding jitter to the locations in the calling routine.
 *
 * The other arguments are
@@ -226,7 +226,7 @@ C first three points collinear:
 *       = 3 for NXI = 0 or less
 *       = 4 for NYI = 0 or less
 *       = 9 for errors in SDTRAN called by this subroutine, except the next one
-*       =10 for error 2 in SDTRAN (first three points collinear), as this can 
+*       =10 for error 2 in SDTRAN (first three points collinear), as this can
 *           be fixed by adding jitter to the locations in the calling routine.
 *
 * The other arguments are
@@ -2339,15 +2339,15 @@ C     WRITE (*,FMT=9040)
 *     .. Array Arguments ..
       DOUBLE PRECISION             XD(NDP),XI(NIP),YD(NDP),
      +     YI(NIP),ZD(NDP),ZI(NIP)
-      INTEGER          IPL(2,NL),IPT(3,NT),ITLI(NIP),KTLI(NIP)
+      INTEGER          IPT(3,NT),ITLI(NIP),KTLI(NIP)
       LOGICAL          EXTRPI(NIP)
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION             A,B,C,EA,EB,EC,DV
-      INTEGER          I,IDP,IIP,ILI,IR,ITLII,ITLIPV,K,KTLII,KTLIPV
+      INTEGER          I,IDP,IIP,ITLII,ITLIPV,KTLII,KTLIPV
 *     ..
 *     .. Local Arrays ..
-      DOUBLE PRECISION             PD(5,3),X(3),Y(3),Z(3),ZU(3),ZUU(3),
+      DOUBLE PRECISION             X(3),Y(3),Z(3),
      +     ZUV(3),ZV(3),ZVV(3)
 *     ..
 *     .. Intrinsic Functions ..
@@ -2375,14 +2375,17 @@ C     WRITE (*,FMT=9040)
                   Y(I) = YD(IDP)
                   Z(I) = ZD(IDP)
  21            CONTINUE
+*     Solve the equations for the plane throuh (x1,y1,z1),
+*     (x2,y2,z2), (x3,y3,z3):
 *     Maxima:
-*     (%i1) g1:a*x1+b*y1+c=z1;
+*     (%i1) eq1:a*x1+b*y1+c=z1;
 *     (%o1)                        b y1 + a x1 + c = z1
-*     (%i2) g2:a*x2+b*y2+c=z2;
+*     (%i2) eq2:a*x2+b*y2+c=z2;
 *     (%o2)                        b y2 + a x2 + c = z2
-*     (%i3) g3:a*x3+b*y3+c=z3;
+*     (%i3) eq3:a*x3+b*y3+c=z3;
 *     (%o3)                        b y3 + a x3 + c = z3
-*     solve([g1,g2,g3],[a,b,c]);
+*     solve([eq1,eq2,eq3],[a,b,c]);
+*     gives:
 
              DV =  X(1)*(Y(3)-Y(2))-X(2)*Y(3)+X(3)*Y(2)+(X(2)-X(3))*Y(1)
              EA =  Y(1)*(Z(3)-Z(2))-Y(2)*Z(3)+Y(3)*Z(2)+(Y(2)-Y(3))*Z(1)
@@ -2394,17 +2397,29 @@ C     WRITE (*,FMT=9040)
                   A = -EA/DV
                   B =  EB/DV
                   C =  EC/DV
+                  ZI(IIP) = A*XI(IIP)+B*YI(IIP)+C
+                  EXTRPI(IIP) = .FALSE.
+               ELSE
+*     Use the extrapolation marker as error message, as no extrapolation
+*     is possible:
+                  ZI(IIP) = 0.0D0
+                  EXTRPI(IIP) = .TRUE.
+               END IF
+            ELSE
+*     reuse A,B,C from last (same) triangle:
+               IF (DABS(DV) .GT. 1.0D-10) THEN
+                  ZI(IIP) = A*XI(IIP)+B*YI(IIP)+C
+                  EXTRPI(IIP) = .FALSE.
                ELSE
-*     TODO:
-                  ZI(IIP)=0
+*     Use the extrapolation marker as error message, as no extrapolation
+*     is possible:
+                  ZI(IIP) = 0.0D0
                   EXTRPI(IIP) = .TRUE.
                END IF
             END IF
-            ZI(IIP) = A*XI(IIP)+B*YI(IIP)+C
-c     ZI(IIP) = (Z(1)+Z(2)+Z(3))/3.0D0
-            EXTRPI(IIP) = .FALSE.
          ELSE
-            ZI(IIP)=0.0D0
+*     outside of triangles:
+            ZI(IIP) = 0.0D0
             EXTRPI(IIP) = .TRUE.
          END IF
  121  CONTINUE