%% -*- mode: Rnw; coding: utf-8; -*- %\VignetteIndexEntry{Interpolation} %\VignetteDepends{scatterplot3d,MASS} %\VignetteKeywords{nonparametric} %\VignettePackage{interp} \documentclass[nojss]{jss} \usepackage[utf8]{inputenc} %\usepackage{Sweave} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsmath} \usepackage{amsthm} \usepackage{flexisym} \usepackage{breqn} \usepackage{bm} \usepackage{graphicx} % put floats before next section: \usepackage[section]{placeins} % collect appendices as subsections \usepackage[toc,page]{appendix} % customize verbatim parts \usepackage{listings} \lstdefinestyle{Sstyle}{ basicstyle=\ttfamily\rsize, columns=fixed, breaklines=true, % sets automatic line breaking breakatwhitespace=false, postbreak=\raisebox{0ex}[0ex][0ex]{\ensuremath{\color{red}\hookrightarrow\space}}, fontadjust=true, basewidth=0.5em, inputencoding=utf8, extendedchars=true, literate={‘}{{'}}1 {’}{{'}}1 % Zeichencodes für Ausgabe von lm() ! {á}{{\'a}}1 {é}{{\'e}}1 {í}{{\'i}}1 {ó}{{\'o}}1 {ú}{{\'u}}1 {Á}{{\'A}}1 {É}{{\'E}}1 {Í}{{\'I}}1 {Ó}{{\'O}}1 {Ú}{{\'U}}1 {à}{{\`a}}1 {è}{{\`e}}1 {ì}{{\`i}}1 {ò}{{\`o}}1 {ù}{{\`u}}1 {À}{{\`A}}1 {È}{{\'E}}1 {Ì}{{\`I}}1 {Ò}{{\`O}}1 {Ù}{{\`U}}1 {ä}{{\"a}}1 {ë}{{\"e}}1 {ï}{{\"i}}1 {ö}{{\"o}}1 {ü}{{\"u}}1 {Ä}{{\"A}}1 {Ë}{{\"E}}1 {Ï}{{\"I}}1 {Ö}{{\"O}}1 {Ü}{{\"U}}1 {â}{{\^a}}1 {ê}{{\^e}}1 {î}{{\^i}}1 {ô}{{\^o}}1 {û}{{\^u}}1 {Â}{{\^A}}1 {Ê}{{\^E}}1 {Î}{{\^I}}1 {Ô}{{\^O}}1 {Û}{{\^U}}1 {œ}{{\oe}}1 {Œ}{{\OE}}1 {æ}{{\ae}}1 {Æ}{{\AE}}1 {ß}{{\ss}}1 {ű}{{\H{u}}}1 {Ű}{{\H{U}}}1 {ő}{{\H{o}}}1 {Ő}{{\H{O}}}1 {ç}{{\c c}}1 {Ç}{{\c C}}1 {ø}{{\o}}1 {å}{{\r a}}1 {Å}{{\r A}}1 {€}{{\euro}}1 {£}{{\pounds}}1 {«}{{\guillemotleft}}1 {»}{{\guillemotright}}1 {ñ}{{\~n}}1 {Ñ}{{\~N}}1 {¿}{{?`}}1 } % switch to above defined style \lstset{style=Sstyle} % nice borders for code blocks \usepackage{tcolorbox} % enable boxes over several pages: \tcbuselibrary{breakable,skins} \tcbset{breakable,enhanced} \definecolor{grey2}{rgb}{0.6,0.6,0.6} \definecolor{grey1}{rgb}{0.8,0.8,0.8} % some abbreviations: \newcommand{\R}{\mathbb{R}} \newcommand{\EV}{\mathbb{E}} \newcommand{\Vect}[1]{\underline{#1}} \newcommand{\Mat}[1]{\boldsymbol{#1}} \newcommand{\Var}{\mbox{Var}} \newcommand{\Cov}{\mbox{Cov}} % lstinline can break code across lines \def\cmd{\lstinline[basicstyle=\ttfamily,keywordstyle={},breaklines=true,breakatwhitespace=false]} % but lstinline generates ugly sectionnames in PDF TOC, so use \texttt there \newcommand{\cmdtxt}[1]{\texttt{#1}} \newtheorem{definition}{Definition}[section] \newtheorem{remark}{Remark}[section] \newtheorem{lemma}{Lemma}[section] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% almost as usual \author{ Albrecht Gebhardt\\ %Department of Statistics, University Klagenfurt \And Roger Bivand\\ %Department of Economics, Norwegian School of Economics} \title{A Re-Implementation of Akima's Spline Interpolation for Scattered Data} %% for pretty printing and a nice hypersummary also set: \Plainauthor{Albrecht Gebhardt, Roger Bivand} %% comma-separated \Plaintitle{A Reimplementation of Akima's Spline Interpolation for Scattered Data} %% a short title (if necessary) %% an abstract and keywords \Abstract{ This vignette presents the \proglang{R} package \pkg{interp} and focuses on interpolation of irregular spaced data. This is the second of planned three vignettes for this package (not yet finished). } \Keywords{interpolation, spline, \proglang{R} software} \Plainkeywords{interpolation, spline, R software} %% without formatting %% at least one keyword must be supplied %% publication information %% NOTE: Typically, this can be left commented and will be filled out by the technical editor % \Volume{XX} %% \Issue{X} %% \Month{XXXXXXX} %% \Year{XXXX} %% \Submitdate{XXXX-XX-XX} %% \Acceptdate{XXXX-XX-XX} %% The address of (at least) one author should be given %% in the following format: \Address{ Albrecht Gebhardt\ Institut für Statistik\\ Universität Klagenfurt\ 9020 Klagenfurt, Austria\\ E-mail: \email{albrecht.gebhardt@aau.at}\ %URL: \url{http://statmath.wu-wien.ac.at/~zeileis/} } %% It is also possible to add a telephone and fax number %% before the e-mail in the following format: %% Telephone: +43/1/31336-5053 %% Fax: +43/1/31336-734 %% for those who use Sweave please include the following line (with % symbols): %% need no \usepackage{Sweave.sty} %% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % for Sinput to set font size of R input code: \newcommand\rsize{% \fontsize{8.5pt}{9.1pt}\selectfont% } \begin{document} % undefine Sinput, Soutput, Scode to be able to redefine them as % \lstnewenvironment{Sinput}... \makeatletter \let\Sinput\@undefined \let\endSinput\@undefined \let\Soutput\@undefined \let\endSoutput\@undefined \let\Scode\@undefined \let\endScode\@undefined \makeatother \hypersetup{pdftitle={Interpolation},pdfauthor={Albrecht Gebhardt and Roger Bivand}, pdfborder=1 1 1 1 1} % Sweave stuff: % graphics dimension: \setkeys{Gin}{width=0.8\textwidth} %\setkeys{Gin}{width=1in} % all in- and output black: \definecolor{Sinput}{rgb}{0,0,0} \definecolor{Soutput}{rgb}{0,0,0} \definecolor{Scode}{rgb}{0,0,0} % redefine Sinput, Soutput, Scode, variant 1 use fancy verbatim % %\DefineVerbatimEnvironment{Sinput}{Verbatim} % gobble=0 !!! otherwise 2 characters of S lines are hidden !!! %{formatcom = {\color{Sinput}},fontsize=\rsize,xleftmargin=2em,gobble=0} %\DefineVerbatimEnvironment{Soutput}{Verbatim} %{formatcom = {\color{Soutput}},fontsize=\rsize,xleftmargin=2em,gobble=0} %\DefineVerbatimEnvironment{Scode}{Verbatim} %{formatcom = {\color{Scode}},fontsize=\rsize,xleftmargin=2em,gobble=0} %\fvset{listparameters={\setlength{\topsep}{0pt}}} %\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}} % % redefine Sinput, Soutput, Scode, variant 2, use color boxes (tcb) \lstnewenvironment{Sinput}{\lstset{style=Sstyle}}{}% \lstnewenvironment{Soutput}{\lstset{style=Sstyle}}{}% \lstnewenvironment{Scode}{\lstset{style=Sstyle}}{}% \renewenvironment{Schunk}{\vspace{\topsep}\begin{tcolorbox}[breakable,colback=grey1]}{\end{tcolorbox}\vspace{\topsep}} % see http://www.stat.auckland.ac.nz/~ihaka/downloads/Sweave-customisation.pdf % % all in one line!!! setting for direct PDF output ! \SweaveOpts{keep.source=TRUE,engine=R,eps=FALSE,pdf=TRUE,strip.white=all,prefix=TRUE,prefix.string=fig-,include=TRUE,concordance=FALSE,width=6,height=6.5} % Sweave initialization: % restrict line length of R output, no "+" for continued lines, % set plot margins: % initialize libraries and RNG if necessary <>= set.seed(42) options(width=80) options(continue=" ") options(SweaveHooks=list(fig=function() par(mar=c(5.1, 4.1, 1.1, 2.1)))) library(interp) @ \section[Note]{Note} \label{sec:note} Notice: This is a preliminary and not yet complete version of this vignette. Finally three vignettes will be available for this package: \begin{enumerate} \item a first one related to partial derivatives estimation, \item this one describing interpolation related stuff \item and a third one dealing with triangulations and Voronoi mosaics. \end{enumerate} \section[Introduction]{Introduction} \label{sec:intro} The main aim of this \proglang{R} package is to provide interpolation algorithms for both regular and irregular data grids $$ \{((x_{i},y_{i})^{\intercal},z_{i})|x_{i},y_{i},z_{i}\in\R \quad i=1,\ldots,n\} $$ From the early days of \proglang{S} and \proglang{S-Plus} there was a function \cmd{interp()} which solved this task. It used Akima's spline interpolation algorithms available at \cmd{netlib}\footnote{\url{https://netlib.org/toms/526.gz}} twice: Once to determine a triangulation of the data which is needed for a piecewise linear interpolation. This is the default application case of this function and as shown in \citet{bivand:17} the most common use of it, especially in other R packages depending on it. Second to get the spline interpolation based an the same triangulation. These algorithms have been available since 1998 in \proglang{R} via the package \cmd{akima}. Unfortunately this package inherits a non-free license from the underlying \proglang{Fortran} code. So the need to rewrite the algorithms under a free license, encouraged by the CRAN team, appeared convincing to the authors of this package. This is now mostly done and package \cmd{interp} provides plugin capable replacement functions for the interpolations delivered in package \pkg{akima}. For both of these interpolations to work it has to be ensured that no duplicate points $(x_{i},y_{i})$ may exist in the given point set $\{(x_{i},y_{i})|i=1,\ldots,n\}$. This is reached via the argument \cmd{duplicate} of \cmd{interp::interp()}. It offers three options: \begin{itemize} \item \cmd{"error"}: Stop with an error, this is the default. \item \cmd{"strip"}: Completely remove points with duplicates, or \item \cmd{"mean"},\cmd{"median"},\cmd{"user"}: apply some function to them. The Implementation provides \cmd{mean()}, \cmd{median()} or a user supplied function (\cmd{"dupfun"}). \end{itemize} \section{Bivariate Linear Interpolation} \label{sec:linear} The default behaviour of the \cmd{interp::interp()} function is to produce a piecewise linear interpolation. This interpolation takes the triangles of the Delaunay triangulation as also returned by \cmd{tri.mesh()} and simply fits a plane to the three vertices $(x_{i},y_{i},z_{i}), i=1,2,3$ of those triangles. As a natural consequence it is not possible to extrapolate this interpolation beyond the convex hull of the given point set. First load the data set used by Akima in his initial work on irregular gridded data \citep{akima:78}, see figure \ref{fig:akima}. <>= data(akima) library(scatterplot3d) scatterplot3d(akima, type="h", angle=60, asp=0.2, lab=c(4,4,0)) @ \begin{figure}[htb] \centering <>= <> @ \caption{Akimas test data in \cite{akima:78}} \label{fig:akima} \end{figure} The next plot in figure \ref{fig:lininterp} shows the linear nature of the isolines of the interpolation generated within all triangles: <>= li <- interp(akima$x, akima$y, akima$z, nx=150, ny=150) MASS::eqscplot(akima$x, akima$y) contour(li, nlevels=30, add=TRUE) plot(tri.mesh(akima), add=TRUE) @ \begin{figure}[htb] \centering <>= <> @ \caption{Piecewise linear interpolation} \label{fig:lininterp} \end{figure} In case the point data set resembles a regular rectangular grid it should be noted that no unique solution to the triangulation task exists. For each rectangle of this grid there are two possibilities to form triangles compatible with the main condition of a Delaunay triangulation: The interior of the circumcircle of each triangle does not contain any other point of the data set. Generally, as long as the data set contains more then 3 points on a common circumcircle which is otherwise empty of remaining points, it will lead to non uniqueness of the triangulation. This in turn means that a piece wise linear interpolation of rectangular gridded data is not unique. Nevertheless \cmd{interp::interp()} will always produce the same result as long as no jitter is applied to the data set. This can be done by explicitly via the argument \cmd{jitter} or it is applied automatically during the underlying triangulation, which applies this in some cases of collinear points to avoid error conditions. \section{Bivariate Spline Interpolation} \label{sec:spline} Akimas spline interpolator 'with the accuracy of a bicubic polynomial' \citep{akima:78a} for irregular gridded data is given by the following polynomial in $x$ and $y$: \begin{equation} \label{eq:akima} p(x,y)=\sum_{i=0}^{5}\sum_{j=0}^{5-i}p_{i,j}x^{i}y^{j} \end{equation} with 21 coefficients $p_{i,j}$, $0\le i\le j\le 5$. This polynomial is determined within each triangle $(v_{1},v_{2},v_{3})$ with vertexes $v_{i}\in\R^{2}, i=1,2,3$ of the Delaunay triangulation. The solution has to fulfill the following restrictions: \begin{enumerate} \item The interpolation itself (condition $(i)$ in \citep{akima:74}) results in 3 conditions. \item First and second order partial derivatives of $p(x,y)$ have to match estimated derivatives at the triangle vertices (Akima denotes them as condition $(ii)$). This makes up for 15 conditions. \item Finally the last three equations (condition $(iii)$) involve the directional derivatives along the normal vectors of the triangle sides. As the spline polynomial is of degree 5 these derivatives generally will be polynomials of degree 4. Now the condition demands that they are polynomials of degree 3 in that variable that is describing the position of that normal vector along the triangle side (later denoted as $s$ in a $(s,t)$ coordinate system), thus setting its highest degree coefficient to zero. This can be expressed by setting the appropriate 4th derivative of this directional derivative to zero. \end{enumerate} The same conditions are also used in an improved algorithm described in \citep{akima:96}, but e.g. the estimation of the partial derivatives is different to the old algorithm and a better triangulation based on the \cmd{TRIPACK} Fortran package has been used \citep{renka:96}. Next we will formulate the conditions at the triangle vertices $\Vect{v_{i}}=(x_i,y_i)^{\intercal}, i=1,2,3$ and for the normal vectors $\Vect{n}_{ij}= \begin{bmatrix} 0&1\\-1&0 \end{bmatrix} \Vect{t}_{ij} $ of the triangle sides $\Vect{t}_{ij}=(x_j,y_j)^{\intercal}-(x_i,y_i)^{\intercal}$ $(i,j)\in\{(1,3),(3,2),(2,1)\}$. \begin{equation} \label{eq:iiiiii} \begin{array}{lrclrclrcl} (i) & p(x_i,y_i)&=&z_i,&\multicolumn{6}{l}{i=1,2,3}\\ (ii)&\frac{\partial}{\partial x}p(x_i,y_i)&=&z_{x,i},& \frac{\partial}{\partial y}p(x_i,y_i)&=&z_{y,i},&\multicolumn{3}{l}{i=1,2,3}\\ &\frac{\partial^2}{\partial x\partial y}p(x_i,y_i)&=&z_{xy,i},& \frac{\partial^2}{\partial x^2}p(x_i,y_i)&=&z_{xx,i},& \frac{\partial^2}{\partial y^2}p(x_i,y_i)&=&z_{yy,i}\\ (iii)&\frac{\partial^{4}}{\partial s^{4}} \Vect{n}_{ij}\nabla p(x,y)&=&0&\multicolumn{6}{l}{(i,j)\in\{(1,3),(3,2),(2,1)\}} \end{array} \end{equation} where $z_{i}$ are the values to interpolate in $\Vect{v}_{i}=(x_{i},y_{i})^{\intercal}, i=1,2,3$ and $z_{x,i}=\frac{\partial}{\partial x}p(x_{i},y_{i})$, $z_{y,i}=\frac{\partial}{\partial y}p(x_{i},y_{i})$, $z_{xx,i}=\frac{\partial^{2}}{\partial x^{2}}p(x_{i},y_{i})$, $z_{xy,i}=\frac{\partial^{2}}{\partial x\partial y}p(x_{i},y_{i})$ and $z_{yy,i}=\frac{\partial^{2}}{\partial^{2} y}p(x_{i},y_{i})$ denote the estimates for partial derivatives at $\Vect{v}_{i}$. Note that the scalar product $\Vect{n}_{ij}\nabla p(x,y)$ represents the directional derivative mentioned above expressed in coordinates $s$ and $t$. All these conditions together ensure that the resulting spline interpolates the given data and the interpolating function is continuous and differentiable across the borders of all triangles. We now illustrate this with the same data set as above in figure \ref{fig:splinterp}. <>= si <- interp(akima$x, akima$y, akima$z, method="akima", nx=150, ny=150) MASS::eqscplot(akima$x, akima$y) contour(si, nlevels=30, add=TRUE) plot(tri.mesh(akima), add=TRUE) @ \begin{figure}[htb] \centering <>= <> @ \caption{Bivariate Spline Interpolation} \label{fig:splinterp} \end{figure} \section{Implementation details} \label{sec:impl} The call to \cmd{interp::interp()} follows this form: \begin{Schunk} \begin{Sinput} interp(x, y = NULL, z, xo = seq(min(x), max(x), length = nx), yo = seq(min(y), max(y), length = ny), linear = (method == "linear"), extrap = FALSE, duplicate = "error", dupfun = NULL, nx = 40, ny = 40, input="points", output = "grid", method = "linear", deltri = "shull", h=0, kernel="gaussian", solver="QR", degree=3, baryweight=TRUE, autodegree=FALSE, adtol=0.1, smoothpde=FALSE, akimaweight=TRUE, nweight=25) \end{Sinput} \end{Schunk} The arguments \cmd{duplicate} and \cmd{dupfun} have been introduced above, as well as \cmd{method} with its currently two available options \cmd{"linear"} and \cmd{"akima"}. Generally the input will be given as three vectors \cmd{x}, \cmd{y} and \cmd{z} of equal length. Omitting \cmd{y} implicates that \cmd{x} consist of a two column matrix or dataframe containing $x$ and $y$ entries. Additionally the argument \cmd{input} has to be set to \cmd{"points"} (which it is by default). If \cmd{input="grid"} is given, \cmd{z} is treated as a matrix of $z$ values containing $z_{i,j}$ for the $x$ and $y$ values given in the argument vectors \cmd{x} and \cmd{y} both of a length matching the dimensions of \cmd{z}. A similar scheme is applied to the output: If \cmd{output="grid"} is set (default) a matrix with rows and columns according to the output defining vectors \cmd{xo} and \cmd{yo} is returned. The output grid can also be specified by setting its dimension to \cmd{nx} times \cmd{ny}, it will then be chosen to cover the range of the input data. With \cmd{output="points"} \cmd{xo} and \cmd{yo} have to be of equal length and only a vector of $z$ values of the same length is returned. Extrapolation (\cmd{extrap=TRUE}) is only possible for spline interpolation but is disabled by default. The remaining parameters control several aspects of the algorithm and are at least partially explained later. Both methods are implemented via the \cmd{Rcpp} interface \citep{rcpp}. As mentioned before, step 1 of these interpolation methods is the Delaunay triangulation, described in another vignette (\cmd{vignette("tri")}) which is based on the sweep hull algorithm described in \citep{sinclair:16}. The access to the triangulation code is done internally via \proglang{C++}, not via the R function \cmd{interp::tri.mesh()}. In the second step the needed estimates for the partial derivatives up to degree 2 in all data points are determined. This is based on a local polynomial regression approach implemented in \proglang{C++}. These intermediate results are also available via \cmd{interp::locpoly()} described in a separate vignette (\cmd{vignette("partDeriv")}). All options of the related \cmd{interp::locpoly()} function are also available in \cmd{interp::interp()}, e.g. argument \cmd{kernel} specifies the kernel used. In contrast to Akima's interpolation we use a gaussian kernel by default and not a uniform one. Argument \cmd{h} contains the bandwidth, either as a scaler, or a vector of length 2. The first setting gives a percentage of the data set used for a local nearest neigbour bendwidth approach. If two bandwidths as a vector are given then two global bandwidths for $x$ and $y$ are chosen as the given percentage of their data range. If \cmd{h=0} then a minimum local bandwidth resulting in 10 nearest neigbours are choosen to be able to determine the 10 parameters of a \cmd{degree=3} polynomial. It is possible to choose different numerical solutions of the weighted least squares method behind the local regression via the argument \cmd{solver} (default is \cmd{"QR"}, but also \cmd{"LLT"}, \cmd{"SVD"}, \cmd{"Eigen"} and \cmd{"CPivQR"} are available) to be used in the local regression step, compare \cmd{fastLm()} in \citep{rcppeigen}. The third step performs the real interpolation. First the estimated derivatives are (optionally) smoothed according to the smoothing scheme detailed in \citep{akima:78}. Then the system of equations (\ref{eq:iiiiii}) is solved per triangle and the results are determined via \begin{dmath} p(x,y)=y\,\left(y\,\left(y\,\left(y\,\left(p_{0,5}\,y+p_{1,4}\,x+p_{0,4}\right)+x\,\left(p_{2,3}\,x+p_{1,3}\right)+p_{0,3}\right)+x\,\left(x\,\left(p_{3,2}\,x+p_{2,2}\right)+p_{1,2}\right)+p_{0,2}\right)+x\,\left(x\,\left(x\,\left(p_{4,1}\,x+p_{3,1}\right)+p_{2,1}\right)+p_{1,1}\right)+p_{0,1}\right)+x\,\left(x\,\left(x\,\left(x\,\left(p_{5,0}\,x+p_{4,0}\right)+p_{3,0}\right)+p_{2,0}\right)+p_{1,0}\right)+p_{0,0} \label{eq:poly} \end{dmath} which is equivalent to (\ref{eq:akima}) but numerically more stable. Optionally some methods to improve the results can be applied. They are choosen via the following arguments: \begin{itemize} \item \cmd{akimaweight}: As mentioned above, this sort of averaging is also done in Akimas original algorithms. It takes by default 25 (parameter \cmd{nweight}) estimates of that specific partial derivative and builds a weighted sum of them with the weights beeing constructed out of normal densities with mean and standard deviations of the according estmation errors. \item \cmd{baryweight}: The system of equations (\ref{eq:iiiiii}) is solved after transforming each triangle into a standardized triangle with vertices $(0,0)^{\intercal}, (1,0)^{\intercal}, (0,1)^{\intercal}$. So one of the three vertices of a triangle gets transformed into $(0,0)^{\intercal}$. During the development of the code it became apperent that the numerical errors for points near to this vertices are minimal and increase for the two other vertices. This weighting scheme repeats the interpolation for all three possibilities to transform a vertex into $(0,0)^{\intercal}$ and then merges the results using the barycentric coordinates (see \ref{sec:baryc-coord}) of the prediction points. That way results generated from a vertex mapped to $(0,0)^{\intercal}$ always dominate and all three vertices can benefit from the reduced numerical errors near $(0,0)^{\intercal}$ after transformation. Clearly this triples the computing time. But nevertheless this option is used by default. As motivation a result with barycentric weighting turned off is given below in figure \ref{fig:splinterpnobw}. <>= si.nobw <- interp(akima$x, akima$y, akima$z, method="akima", nx=150, ny=150, baryweight=FALSE) MASS::eqscplot(akima$x, akima$y) contour(si.nobw, nlevels=30, add=TRUE) plot(tri.mesh(akima), add=TRUE) @ The plot clearly shows (e.g. in the center of the upper left quadrant) the numerical problems of disconnected isolines across the triangle borders. Note, that these errors occur only on one triangle edge. It turned out this is opposite to the vertex mapped internally by the algorithm to $(0,0)^{\intercal}$. So we encourage to use this option even dispite the tripled computing time. Only if acurracy does not really matter one could reduce the computing time by turning it off. \begin{figure}[htb] \centering <>= <> @ \caption{Bivariate Spline Interpolation (Without barycentric weighting)} \label{fig:splinterpnobw} \end{figure} \item \cmd{smoothpde}: If \cmd{TRUE} smoothing of partial derivative estimates, if \cmd{akimaweight==TRUE} then Akimas weighting scheme is applied, otherwise a simple arithmetic mean is returned. Note that it is disabled by default which in turn means that also no Akima weighting is applied. If it is enabled then Akima weighting is used by default and a simple arithmetic mean if \cmd{akimaweight=FALSE} is given. \item \cmd{autodegree}: If the variability of the interpolates is above \cmd{adtol} then reduce the degree of the polynomial to get a smoother result. This is also disabled by default. \end{itemize} If \cmd{interp::interp()} is called with regular gridded data as input, it uses the same irregular grid based algorithm. This is in contrast to the old package \cmd{akima}, this also contained Akimas code for regular gridded data, based on \citep{akima:74} and \citep{akima:96a}. Maybe a future version of package \cmd{interp} will also contain a re-implementation of this old code. This package also implements bilinear interpolation for rectangular grids. Given a rectangle $\{(x_{1},y_{1})^{\intercal},(x_{2},y_{2})^{\intercal},(x_{3},y_{3})^{\intercal},(x_{4},y_{4})^{\intercal}\}$ and $y_{1}=y_{2}$, $y_{3}=y_{4}$, $x_{1}=x_{4}$ and $x_{2}=x_{3}$ (this makes it axis parallel) with counter clockwise indexed vertexes and according $z$ values $z_{1},z_{2},z_{3},z_{4},$ this algorithm can be described as follows: For a location $(x_{0},y_{0})^{\intercal}$ contained in this rectangle the interpolation is determined via: \begin{enumerate} % \item Calculate intermediate vertexes % $$(x_{12},y_{12})^{\intercal}=\frac{x_{0}-x_{1}}{x_{2}-x_{1}}((x_{2},y_{2})^{\intercal}+(x_{2},y_{2})^{\intercal}) % \quad\mbox{and}\quad % (x_{34},y_{34})^{\intercal}=\frac{x_{0}-x_{1}}{x_{2}-x_{1}}((x_{3},y_{3})^{\intercal}+(x_{4},y_{4})^{\intercal}).$$ \item Determine intermediate $z$ values for $(x_{0},y_{1})^{\intercal}$ and $(x_{0},y_{3})^{\intercal}$ as $$z_{01}=\frac{x_{0}-x_{1}}{x_{2}-x_{1}}(z_{1}+z_{2}) \quad\mbox{and}\quad z_{03}=\frac{x_{0}-x_{1}}{x_{2}-x_{1}}(z_{3}+z_{4}).$$ \item Now get $$z_{0}=\frac{y_{0}-y_{1}}{y_{4}-y_{1}}(z_{01}+z_{03}).$$ \end{enumerate} This results in a polynomial of degree 2 which is continuous but not differentiable at the borders of the rectangle. We use Franke function 1 \citep{franke:82} on a regular grid for the demonstration, see figure \ref{fig:bilinear}. <>= nx <- 8; ny <- 8 xg<-seq(0,1,length=nx) yg<-seq(0,1,length=ny) xyg<-expand.grid(xg,yg) fg <- outer(xg,yg,function(x,y)franke.fn(x,y,1)) # not yet implemented this way: # bil <- interp(xg,yg,fg,input="grid",output="grid",method="bilinear") bil <- bilinear.grid(xg, yg, fg, dx=0.01, dy=0.01) MASS::eqscplot(xyg[,1], xyg[,2]) contour(bil, add=TRUE) @ \begin{figure}[htb] \centering <>= <> @ \caption{Bilinear interpolation of regularly gridded data} \label{fig:bilinear} \end{figure} % FIXME: index bug in \cmd{BiLinear}: % <>= % bil <- BiLinear.grid(xg, yg, fg, dx=0.01, dy=0.01) % MASS::eqscplot(xyg[,1], xyg[,2]) % contour(bil, add=TRUE) % @ \section{One-Dimensional Data} \label{sec:1d} Akima also implemented algorithms for one-dimensional spline interpolation, see \citep{akima:72}. So it was a natural choice to include these algorithms also in the package \pkg{akima}. The functions \cmd{aspline()} and \cmd{aSpline()} are freely licensed re-implementations of this algorithm in \proglang{Fortran} and \proglang{C++}. It comes in two versions, one as described in \citep{akima:72} and an improved version as described in \citep{akima:91}, the newer algorithm also allows for higher degrees of the polynomial, not only degree 3, compare figure \ref{fig:aspline} <>= x <- c(-3, -2, -1, 0, 1, 2, 2.5, 3) y <- c( 0, 0, 0, 0, -1, -1, 0, 2) MASS::eqscplot(x, y, ylim=c(-2, 3)) lines(aspline(x, y, n=200, method="original"), col="red") lines(aspline(x, y, n=200, method="improved"), col="black", lty="dotted") lines(aspline(x, y, n=200, method="improved", degree=10), col="green", lty="dashed") @ \begin{figure}[htb] \centering <>= <> @ \caption{Spline interpolation of onedimensional data} \label{fig:aspline} \end{figure} \section{Appendix} \label{sec:appendix} \subsection{Barycentric Coordinates} \label{sec:baryc-coord} Points within a triangle can be expressed in barycentric coordinates as follows: Given a triangle with vertices $\Vect{v}_{i}=(x_i,y_i)^{\intercal}, i=1,2,3$ any interior point $\Vect{v}_{0}=(x_0,y_0)^{\intercal}$ of this triangle can be expressed as a convex linear combination $$ \Vect{v}_{0}=a\cdot\Vect{v}_{1}+b\cdot\Vect{v}_{2}+c\cdot\Vect{v}_{3} $$ with $a,b,c\in [0,1]$ and $a+b+c=1$ (notation: $[a:b:c]$). The vertices itself carry the representation $[1:0:0]$ , $[0:1:0]$ and $[0:0:1]$. In section \ref{sec:impl} we used these coordinates to build a weighted sum of three interpolation results. Component $a$ of the barycentric coordinates of a point near vertex $\Vect{v_1}$ will be close to 1 and so the interpolation result with the lowest numerical error (where vertex $\Vect{v}_{1}$ had been transformed to $(0,0)^{\intercal}$) will dominate the barycentric weighted sum mentioned above. Using this approach we cherry pick the numerically best portions of these three interpolation results. \bibliography{lit} %\addcontentsline{toc}{section}{Tables} %\listoftables \addcontentsline{toc}{section}{Figures} \listoffigures \end{document}