Interpolating non-gridded data
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 grid x <- y <- seq(-3, 3, by = 0.1) Grid <- expand.grid(X = x, Y = y) # make surface z <- with(Grid, 30 * dcauchy(X) * dcauchy(Y)) # plot library(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 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] 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 hole GridWithHole <- rbind(Constraints, GridWithHole) # remove duplicated points GridWithHole <- GridWithHole[!duplicated(GridWithHole), ] # vertices points <- cbind(GridWithHole[["X"]], GridWithHole[["Y"]]) # constraint edges must be given by vertex indices edges <- rbind(c(1L, 2L), c(2L, 3L), c(3L, 4L), c(4L, 1L)) # the hole edges <- rbind(edges, edges + 4L) # the outer square # perform constrained Delaunay triangulation library(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 plot z <- 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 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) fun <- interpfun( GridWithHole[["X"]], GridWithHole[["Y"]], z, method = "linear" ) # make a grid of the hole xnew <- ynew <- seq(-0.5, 1, by = 0.04) griddedHole <- expand.grid(X = xnew, Y = ynew) # interpolation znew <- fun(griddedHole[["X"]], griddedHole[["Y"]]) # new 3D data newPoints <- cbind(griddedHole[["X"]], griddedHole[["Y"]], znew) # 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:
fun <- interpfun( GridWithHole[["X"]], GridWithHole[["Y"]], z, method = "sibson" )
And we obtain:
This is not exactly the true curve, but nevertheless this is impressive.
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.