Fun with Catmull-Rom splines using GNUplot and Ruby

By Greg Goltsov

So recently we were studying all kinds of curves: starting with Bézier, then looking at Catmull-Rom and B-splines. One of the assignments we had was to actually plot Catmull-Rom spline, given a few control points. Now, the problem is quite simple, but I spent too much time thinking of what to learn in order to plot it. R? It's a bit ugly. Matlab? Can't afford a copy for home. Julia looked very nice, but I wasn't really in the mood for yet another language, especially with a bit of a shift if programming paradigm (very procedural). Python also looked great (oh, NumPy), but I've already heavily invested in Ruby; although I'll probably end up needing Python at some later time.

And then I realised - why not just use what I know and love? Unfortunately, there isn't much going on in scientific/technical scene in Ruby — SciRuby is really the only project that pops up in search.

So, really, I decided to go somewhat low-level and just compute everything by hand, numerically, and plot everything using the very sturdy GNUplot.

A bit of theory

The are a few reasons why Catmull-Rom splines are nice: they interpolate all the points, so the curve actually passes through them (as opposed to Bézier spline). Also they are quite easy to compute and preserve the tangents over multiple segments.

So, what exactly does it look like? Well, the parametric equation for a Catmull-Rom spline, \( \mathbf{S}(t) \), is

$$\mathbf{S}(t) = \mathbf{c}_{0} + \mathbf{c}_{1}t + \mathbf{c}_{2}t^2 + \mathbf{c}_{3}t^3$$

Since it's most usually used as a cubic interpolation, it takes in 4 points, meaning that any complex curve will have to be split into segments (we'll cover that later). Also, one of the peculiarities of CR spline is that it doesn't interpolate the first and the last points (see the figure at the beginning).

Now, the cool thing is that to get the coefficients for the parametric equation, we just use a simple matrix multiplication over all the segments,

$$\begin{pmatrix} {\mathbf{c}}_{0} \\ {\mathbf{c}}_{1} \\ {\mathbf{c}}_{2} \\ {\mathbf{c}}_{3} \end{pmatrix}=\begin{pmatrix} 0 & 1 & 0 & 0 \\ -\tau & 0 & \tau & 0 \\ 2\tau & \tau -3 & 3-2\tau & -\tau \\ -\tau & 2-\tau & \tau -2 & \tau \end{pmatrix}\begin{pmatrix} {\mathbf{p}}_{i-2} \\ {\mathbf{p}}_{i-1} \\ {\mathbf{p}}_{i} \\ {\mathbf{p}}_{i+1} \end{pmatrix}$$

And practice

ruby_gnuplot provides a nice interface to GNUplot's process and Plot and SPlot objects. So, on the most basic level,

For every curve we're plotting,

  1. Create a new DataSet object
  2. Split the given control points into segments
  3. Generate a parametric expression for each segment
  4. Sample the parametric equation with values from \(t=0\) to \(t=1\) with a given step \(t_{step}\)
  5. Plot the resulting datapoints

Using this simple procedure, we can get quite a decent-looking graph (this one uses an svg terminal, as pngcairo isn't working for me, for some strange reason; you can read more about output terminals here).

Here, the control points are barely visible grey crosses.

Going further

But you may ask, "why did you bother with choosing Ruby?". Well, having a simple plotter is great, but having a full-blown general-purpose programming language lets you do more cool things. For example, we know that CR splines are parametrised basically by only one parameter, \(\tau\), so what happens when we change it? The answer is a simple for loop!

(0..5).step(0.1) do |τ|
  p.data <<
    Gnuplot::DataSet.new( CR_curve(0.01, points, τ) ) do |ds|
      ds.notitle
      ds.with = 'lines'
      ds.linewidth = 2
    end
end

Here, we're generating a separate CR spline for 50 separate \(\tau\)-values. And here's the result:

But the \([0, 5]\) range is all nice and well, but what if we go to some crazy values, say, 50? Boom!

(0..50).step(1) do |τ|
  p.data <<
    Gnuplot::DataSet.new( CR_curve(0.01, points, τ) ) do |ds|
      ds.notitle
      ds.with = 'lines'
      ds.linewidth = 2
    end
end

You can just barely see the small red triangle in the middle - that's our initial \(\tau = 0\) plot.

Going further, what if we generate a lot of curves, and make the colour dependent on the value of \(\tau\)? Sure, can do that:

# This is a double plot, gives a black and white gradient plot
# with each curve's blackness depending on how close to 1 τ is.
(0..1).step(0.01) do |τ|
  p.data <<
    Gnuplot::DataSet.new( CR_curve(0.01, points, τ) ) do |ds|
      ds.notitle
      ds.with = 'lines'
      ds.linewidth = 2
      # 0.2 is just to offset the min value, so it's not pure white
      colour = ColorMath::HSL.new(0, 0, (1-τ+0.2)/1.0).hex
      ds.linecolor = "rgb '#{colour}'"
    end
end
(1..5).step(0.01) do |τ|
  p.data <<
    Gnuplot::DataSet.new( CR_curve(0.01, points, τ) ) do |ds|
      ds.notitle
      ds.with = 'lines'
      ds.linewidth = 2
      colour = ColorMath::HSL.new(0, 0, (τ-1)/(5-1)).hex
      ds.linecolor = "rgb '#{colour}'"
    end
end

And this gives us a very nice, smooth, gradient:

I know I could've used GNUplot's colorbox for some colour gradients, but I just wanted to keep it quite minimal.

So, that's all for now. This was a small hack that I wrote to see how CR splines behave, but it led to quite a lot of experimentation, so I thought I'd share. I'm planning on implementing the Bézier curve plotter soon.

How do I get this?

All the code is available on GitHub. Feel free to fork, edit and experiment!