Calculate a 2D spline curve in R

Question

I'm trying to calculate a Bezier-like spline curve that passes through a sequence of x-y coordinates. An example would be like the following output from the cscvn function in Matlab (example link):

I believe the (no longer maintained) grid package used to do this (grid.xspline function?), but I haven't been able to install an archived version of the package, and don't find any examples exactly along the lines of what I would like.

The bezier package also looks promising, but it is very slow and I also can't get it quite right:

library(bezier)

set.seed(1)
n <- 10
x <- runif(n)
y <- runif(n)
p <- cbind(x,y)
xlim <- c(min(x) - 0.1*diff(range(x)), c(max(x) + 0.1*diff(range(x))))
ylim <- c(min(y) - 0.1*diff(range(y)), c(max(y) + 0.1*diff(range(y))))
plot(p, xlim=xlim, ylim=ylim)
text(p, labels=seq(n), pos=3)

bp <- pointsOnBezier(cbind(x,y), n=100)
lines(bp$points)
arrows(bp$points[nrow(bp$points)-1,1], bp$points[nrow(bp$points)-1,2],
  bp$points[nrow(bp$points),1], bp$points[nrow(bp$points),2]
)

As you can see, it doesn't pass through any points except the end values.

I would greatly appreciate some guidance here!

Unruly

On a technical note, what you're looking for is called a Catmull-Rom spline. It's in the same class of curves as Bezier curves but has different "which points define what" properties. In this case, each segment between points n and n+1 has a tangent parallel to the line (n-1)--(n+1) at point n, and n--(n+2) at point `(n+1), with the length of the tangent vector determined by the "tension" you choose to work with. At tension 1/2, it's an equivalence curve with cubic Beziers (and yes that means the tangents at the first and last point need to be "made up" =) – Brewery

Answer 1

There is no need to use grid really. You can access xspline from the graphics package.

Following from your code and the shape from @mrflick:

set.seed(1)
n <- 10
x <- runif(n)
y <- runif(n)
p <- cbind(x,y)
xlim <- c(min(x) - 0.1*diff(range(x)), c(max(x) + 0.1*diff(range(x))))
ylim <- c(min(y) - 0.1*diff(range(y)), c(max(y) + 0.1*diff(range(y))))
plot(p, xlim=xlim, ylim=ylim)
text(p, labels=seq(n), pos=3)

You just need one extra line:

xspline(x, y, shape = c(0,rep(-1, 10-2),0), border="red")

Comeau

Answer 2

It may not the be the best approach, bit grid certainly isn't inactive. It's included as a default package with the R installation. It's the underlying graphics engine for plotting libraries like lattice and ggplot. You shouldn't need to install it, you should just be able to load it. Here's how I might translate your code to use grid.xpline

set.seed(1)
n <- 10
x <- runif(n)
y <- runif(n)
xlim <- c(min(x) - 0.1*diff(range(x)), c(max(x) + 0.1*diff(range(x))))
ylim <- c(min(y) - 0.1*diff(range(y)), c(max(y) + 0.1*diff(range(y))))

library(grid)
grid.newpage()
pushViewport(viewport(xscale=xlim, yscale=ylim))
grid.points(x, y, pch=16, size=unit(2, "mm"), 
    default.units="native")
grid.text(seq(n), x,y, just=c("center","bottom"), 
    default.units="native")
grid.xspline(x, y, shape=c(0,rep(-1, 10-2),0), open=TRUE, 
    default.units="native")
popViewport()

which results in

note that grid is pretty low-level so it's not super easy to work with, but it does allow you far more control of what and where you plot.

And if you want to extract the points along the curve rather than draw it, look at the ?xsplinePoints help page.

Gothic

You can extend the range of a vector with extendrange. – Seleta

@VincentGuillemot - Thanks for the helpful function – Unruly

@Gothic - Great answer as usual - thanks so much. This will definitely get me started. I just need to get used to using those low-level functions. Cheers. – Unruly

there's something odd going on with that segment from point 6 to point 7; give the tangent at p7 parallel to the line p6--p8, there should not be any inflections in that segments, but it looks like there are. – Brewery

Answer 3

Thanks to all that helped with this. I'm summarizing the lessons learned plus a few other aspects.

Catmull-Rom spline vs. cubic B-spline

Negative shape values in the xspline function return a Catmull-Rom type spline, with spline passing through the x-y points. Positive values return a cubic B type spline. Zero values return a sharp corner. If a single shape value is given, this is used for all points. The shape of end points is always treated like a sharp corner (shape=0), and other values do not influence the resulting spline at the end points:

# Catmull-Rom spline vs. cubic B-spline
plot(p, xlim=extendrange(x, f=0.2), ylim=extendrange(y, f=0.2))
text(p, labels=seq(n), pos=3)
# Catmull-Rom spline (-1)
xspline(p, shape = -1, border="red", lwd=2) 
# Catmull-Rom spline (-0.5)
xspline(p, shape = -0.5, border="orange", lwd=2) 
# cubic B-spline (0.5)
xspline(p, shape = 0.5, border="green", lwd=2) 
# cubic B-spline (1)
xspline(p, shape = 1, border="blue", lwd=2)
legend("bottomright", ncol=2, legend=c(-1,-0.5), title="Catmull-Rom spline", col=c("red", "orange"), lty=1)
legend("topleft", ncol=2, legend=c(1, 0.5), title="cubic B-spline", col=c("blue", "green"), lty=1)

Extracting results from xspline for external plotting

This took some searching, but the trick is to apply the argument draw=FALSE to xspline.

# Extract xy values
plot(p, xlim=extendrange(x, f=0.1), ylim=extendrange(y, f=0.1))
text(p, labels=seq(n), pos=3)
spl <- xspline(x, y, shape = -0.5, draw=FALSE) 
lines(spl)
arrows(x0=(spl$x[length(spl$x)-0.01*length(spl$x)]), y0=(spl$y[length(spl$y)-0.01*length(spl$y)]),
       x1=(spl$x[length(spl$x)]), y1=(spl$y[length(spl$y)])
)

Unruly