# Interpolating non-gridded data

I still wrapped a part of the C++ library **CGAL** in a R package, namely **interpolation**.

The purpose of this package is to perform interpolation of bivariate functions. As compared to existing packages, it can do more: it can interpolate vector-valued functions (with dimension two or three), and it does not require that the given data are gridded. I will illustrate this second point here.

First, let’s plot a surface \(z = f(x, y)\).

```
# make a square grid
y <- seq(-3, 3, by = 0.1)
x <- expand.grid(X = x, Y = y)
Grid <-# make surface
with(Grid, 30 * dcauchy(X) * dcauchy(Y))
z <-# plot
library(deldir)
deldir(Grid[["X"]], Grid[["Y"]], z = z)
delxyz <-library(rgl)
open3d(windowRect = 50 + c(0, 0, 512, 512))
persp3d(delxyz, color = "red")
```

Now we will make a hole in this surface, and then we will interpolate it.

```
# make a hole in the grid: the square [-0.5, 1]x[-0.5, 1]
with(Grid, X > -0.5 & X < 1 & Y > -0.5 & Y < 1)
toremove <- Grid[!toremove, ]
GridWithHole <-# plot this grid
plot(GridWithHole[["X"]], GridWithHole[["Y"]], pch = 19, asp = 1)
```

Now, to plot the surface with the hole, I will use a constrained Delaunay triangulation. I didn’t find a more straightforward way.

```
# constraint edges: squares [-0.5, 1]x[-0.5, 1] and [-3, 3]x[-3, 3]
data.frame(
Constraints <-X = c(-0.5, 1, 1, -0.5, -3, 3, 3, -3),
Y = c(-0.5, -0.5, 1, 1, -3, -3, 3, 3)
)# add constraint edges to the grid with the hole
rbind(Constraints, GridWithHole)
GridWithHole <-# remove duplicated points
GridWithHole[!duplicated(GridWithHole), ]
GridWithHole <-# vertices
cbind(GridWithHole[["X"]], GridWithHole[["Y"]])
points <-# constraint edges must be given by vertex indices
rbind(c(1L, 2L), c(2L, 3L), c(3L, 4L), c(4L, 1L)) # the hole
edges <- rbind(edges, edges + 4L) # the outer square
edges <-# perform constrained Delaunay triangulation
library(delaunay)
delaunay(points, constraints = edges) del <-
```

Note that the **delaunay** package is also a wrapper of **CGAL**.

This Delaunay triangulation provides triangular faces that we can use to create a 3D **rgl** mesh.

```
# create surface plot
with(GridWithHole, 30 * dcauchy(X) * dcauchy(Y))
z <- cbind(points, z)
vertices <- tmesh3d(
rmesh <-vertices = t(vertices),
indices = t(del[["faces"]]),
homogeneous = FALSE
) addNormals(rmesh)
rmesh <-# we plot the front side in red and the back side in gray
open3d(windowRect = 50 + c(0, 0, 512, 512))
shade3d(rmesh, color = "red", back = "cull")
shade3d(rmesh, color = "gray", front = "cull")
persp3d(delxyz, add = TRUE, alpha = 0.2)
```

Good. Now let’s interpolate.

```
# make the interpolating function
library(interpolation)
interpfun(
fun <-"X"]], GridWithHole[["Y"]], z, method = "linear"
GridWithHole[[
)# make a grid of the hole
ynew <- seq(-0.5, 1, by = 0.04)
xnew <- expand.grid(X = xnew, Y = ynew)
griddedHole <-# interpolation
fun(griddedHole[["X"]], griddedHole[["Y"]])
znew <-# new 3D data
cbind(griddedHole[["X"]], griddedHole[["Y"]], znew)
newPoints <-# plot
open3d(windowRect = 50 + c(0, 0, 512, 512))
shade3d(rmesh, color = "red", back = "cull")
shade3d(rmesh, color = "gray", front = "cull")
points3d(newPoints)
```

Not very nice, you think? Right, but I used the linear method of interpolation here. The **interpolation** package also provides the *Sibson* method, this one is not linear. One just has to repeat the above code but starting with:

```
interpfun(
fun <-"X"]], GridWithHole[["Y"]], z, method = "sibson"
GridWithHole[[ )
```

And we obtain:

This is not exactly the true curve, but nevertheless this is impressive.