[This article was first published on Saturn Elephant, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't. 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 gridx <- y <- seq(-3, 3, by = 0.1)Grid <- expand.grid(X = x, Y = y)# make surface z <- with(Grid, 30 * dcauchy(X) * dcauchy(Y))# plotlibrary(deldir)delxyz <- deldir(Grid[["X"]], Grid[["Y"]], z = z)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]toremove <- with(Grid, X > -0.5 & X < 1 & Y > -0.5 & Y < 1)GridWithHole <- Grid[!toremove, ]# plot this gridplot(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]Constraints <- data.frame( 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 holeGridWithHole <- rbind(Constraints, GridWithHole)# remove duplicated pointsGridWithHole <- GridWithHole[!duplicated(GridWithHole), ]# verticespoints <- cbind(GridWithHole[["X"]], GridWithHole[["Y"]])# constraint edges must be given by vertex indicesedges <- rbind(c(1L, 2L), c(2L, 3L), c(3L, 4L), c(4L, 1L)) # the holeedges <- rbind(edges, edges + 4L) # the outer square# perform constrained Delaunay triangulationlibrary(delaunay)del <- delaunay(points, constraints = edges)
```
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 plotz <- with(GridWithHole, 30 * dcauchy(X) * dcauchy(Y))vertices <- cbind(points, z)rmesh <- tmesh3d( vertices = t(vertices), indices = t(del[["faces"]]), homogeneous = FALSE)rmesh <- addNormals(rmesh)# we plot the front side in red and the back side in grayopen3d(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 functionlibrary(interpolation)fun <- interpfun( GridWithHole[["X"]], GridWithHole[["Y"]], z, method = "linear")# make a grid of the holexnew <- ynew <- seq(-0.5, 1, by = 0.04)griddedHole <- expand.grid(X = xnew, Y = ynew)# interpolationznew <- fun(griddedHole[["X"]], griddedHole[["Y"]])# new 3D datanewPoints <- cbind(griddedHole[["X"]], griddedHole[["Y"]], znew)# plotopen3d(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:
fun <- interpfun( GridWithHole[["X"]], GridWithHole[["Y"]], z, method = "sibson")
And we obtain:
This is not exactly the true curve, but nevertheless this is impressive.
To leave a comment for the author, please follow the link and comment on their blog: Saturn Elephant.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.Continue reading: Interpolating non-gridded data