Parametric Cubic Spline Tutorial

[VBAhack]

Parametric equations are powerful and flexible. Perhaps the most familiar example is the equation of a circle in the form x = r*cos(θ), y = r*sin(θ). In this case, the parameter θ is the independent variable and increases monotonically (i.e. each successive value is larger than the previous one, which is a requirement of parametric equations in general), and x and y are dependent variables. The resulting plot of y vs x can wrap on itself, which a circle indeed does. There are many forms of polynomial parametric equations (e.g. Bezier, B-Spline, NURBS, Catmull-Rom) that have entire books devoted to their understanding and use. In this tutorial, we offer a brief introduction by way of examples to provide a basic understanding of the relative simplicity of parametric equations. Specifically, we will build upon the tools learned in the cubic spline tutorial (http://www.vbforums.com/showthread.php?t=480806) to extend our knowledge to parametric cubic splines.

Recall that a cubic spline is nothing more than a sequence of 3^rd order polynomials joined at the endpoints with enforced 1^st and 2^nd derivative compatibility at interior points and specified end conditions at the free ends. Extending this concept to parametric splines just means formulating two sets of equations instead of one using the exact same methodology as a standard (non parametric) cubic spline. In the case of parametric cubic splines, each spline segment is represented by 2 equations in an independent variable s:

x = f₁(s) = a_x(s-s₀)³ + b_x(s-s₀)² + c_x(s-s₀) + d_x
y = f₂(s) = a_y(s-s₀)³ + b_y(s-s₀)² + c_y(s-s₀) + d_y

where s₀ represents the value of the independent variable s at the beginning of the segment. Though not required, it is convenient to have s vary from 0 to 1. To get an idea of how powerful these equations are, consider the following simple equations for x & y:

x = 26s³ - 40s² + 15s - 1
y = -4s² + 3s

The resulting plot of y vs x (for s = 0 to 1) is truly amazing. It almost seems magical that a single cubic polynomial (well, two actually) can generate such a complicated shape. That’s the power of parametric equations.

Example Problem #1 (Connecting Spline Segments)

Let’s focus on joining several parametric cubic spline segments to form a shape. Consider the following 4 points x = [(1,0), (0,1), (-1,0), (0,-1)]. There are several ways the points can be connected to form a shape. We wish to connect them in a crossing fashion as illustrated in the picture.
.

[VBAhack]

Parametric Cubic Spline continued...

The methodology for determining the coefficients of the parametric cubic polynomials is exactly the same as before. The trick is to determine the values of the independent variable, s. One method is to use distance along the curve (arc length), or an approximation thereof. The values of the parameter s will be formed by summing the distances between the points, (Δx² + Δy²)^1/2. The key is to sum the distances in the sequence that the curve is to follow. In this case:

p₁ = [1,0]
p₂ = [-1, 0]
p₃ = [0, 1]
p₄ = [0, -1]

Calculating the distances between points:

d₁₂ = sqrt[(-1 - 1)² + (0 – 0)²] = sqrt[4 + 0] = 2 ==> s₂ = 2
d₂₃ = sqrt[(-1 – 0)² + (0 – 1)²] = sqrt[1 + 1] = 1.414 ==> s3 = s2 + 1.414 = 3.414
d₃₄ = sqrt[(0 – 0)² + (-1 – 1)²] = sqrt[0 + 4] = 2 ==> s4 = s3 + 2 = 5.414

And the starting point is s₁ = 0. As stated before, though not required, it is convenient (mainly for plotting purposes) to scale the values of s to range from 0 to 1. In this case, we just divide each value by s₄ (i.e. 5.414):

s₁ = 0
s₂ = 0.369
s₃ = 0.631
s₄ = 1

This gives us two sets of points, each of which will be fit with a cubic spline:

[s, x] = [(0, 1), (0.369, -1), (0.631, 0), (1, 0)]
[s, y] = [(0, 0), (0.369, 0), (0.631, 1), (1, -1)]

Recalling the matrix equation of a cubic spline, we have the following for n = 4 (Note: . means zero):

┌				┐ ┌       ┐    ┌                                            ┐
│ ?      ?           ?      ?	│ │x1" y1"│    │             ?                  ?           │
│h1   2(h1+h2)      h2	    .	│ │x2" y2"│ = 6│(x3-x2)/h2-(x2-x1)/h1  (y3-y2)/h2-(y2-y1)/h1│
│ .      h2      2(h2+h3)  h3	│ │x3" y2"│    │(x4-x3)/h2-(x3-x2)/h1  (y4-y3)/h2-(y3-y2)/h1│
│ ?      ?           ?      ?  	│ │x4" y3"│    │             ?                  ?           │
└   				┘ └       ┘    └                                            ┘

This is identical to the non-parametric case, except now we have two sets of 2^nd derivatives to determine and two right hand sides. The base matrix is the same for both. Note in this case, the h_i are formed by subtracting successive values of s_i (e.g. h₁ = s₂ – s₁), the independent variable. We can calculate the needed quantities:

h₁ = s₂ – s₁ = 0.369 – 0 = 0.369
h₂ = s₃ – s₂ = 0.631 – 0.369 = 0.261
h₃ = s₄ – s₃ = 1 – 0.631 = 0.369

2(h₁ + h₂) = 2(0.369 + 0.261) = 1.261
2(h₂ + h₃) = 2(0.261 + 0.369) = 1.261

x₂ – x₁ = -1 – 1 = -2
x₃ – x₂ = 0 – (- 1) = 1
x₄ – x₃ = 0 – 0 = 0

y₂ – y₁ = 0 – 0 = 0
y₃ – y₂ = 1 – 0 = 1
y₄ – y₃ = (-1) – 1 = -2

(x₃ - x₂) / h₂ – (x₂ - x₁) / h₁ = 1 / 0.261 – (-2) / 0.369 = 9.243
(x₄ - x₃) / h₃ – (x₃ - x₂) / h₂ = 0 / 0.369 – 1 / 0.261 = -3.828

(y₃ - y₂) / h₂ – (y₂ - y₁) / h₁ = 1 / 0.261 – 0 / 0.369 = 3.828
(y₄ - y₃) / h₃ – (y₃ - y₂) / h₂ = (-2) / 0.369 – 1 / 0.261 = -9.243

Choosing natural boundary conditions for x" (i.e. x₁" = x₄" = 0) and specific values for y" (to obtain a more pleasing plot), we get (Note: . means zero):

┌				┐ ┌       ┐    ┌               ┐
│1       .       .       .    	│ │x1" y1"│    │ 0        1    │
│0.369   1.261   0.261   .	│ │x2" y2"│ = 6│ 9.243    3.828│
│.       0.261   1.261   0.369	│ │x3" y2"│    │-3.828   -9.243│
│.       .       .       1  	│ │x4" y3"│    │ 0        1    │
└   				┘ └       ┘    └               ┘

x" = [0, 49.882, -28.544, 0,]^T
y" = [6, 27.088, -51.338, 6]^T

from which we can calculate the spline segment coefficients (one set for x and another set for y) and plot the result (left hand picture). By varying the sequence in which the points are connected and playing with the end conditions, interesting results can be obtained.
.

[VBAhack]

Example Problem #2 (Closed Curve)

Let us connect the points in the same sequence as before, but close the loop. In this case, we need a 5^th point in the sequence that is identical to the 1^st:

p₁ = [1,0]
p₂ = [-1, 0]
p₃ = [0, 1]
p₄ = [0, -1]
p₅ = p₁ = [1,0]

For s, arc length approximation could be used as before, but an easier approach will be used so the simplified equations can be employed. We will choose equal spacing: s = [0, 0.25, 0.5, 0.75, 1].

Boundary conditions require special consideration. In the case of a closed loop, the 1^st and last points are joined, which makes them interior points subject to the same 1^st and 2^nd derivative compatibility constraints as the other interior points. Thus, with closed loops, there are no choices to be made for boundary conditions! To derive the 1^st derivative compatibility equation, simply equate 1^st derivatives of p₁ and p₅:

p₁' = p₅'
p₁' = (p₂ – p₁)/h₁ - p₂"h₁/6 - p₁"h₁/3 = p₅' = (p₅ – p₄)/h₄ + p₅"h₄/3 + p₄"h₄/6
6(p₂ – p₁)/h₁ - p₂"h₁ – 2p₁"h₁ = 6(p₅ – p₄)/h₄ + 2p₅"h₄ + p₄"h₄
6[(p₂ – p₁)/h₁ - (p₅ – p₄)/h₄] = p₂"h₁ + 2p₁"h₁ + 2p₅"h₄ + p₄"h₄

The 2^nd derivative boundary conditions are easy: p₁" - p₅" = 0. (Note: using p just means there are 2 sets of equations, one for x and another for y). Also, h₁ = h₂ = h₃ = h₄ = h = s₂ = 0.25, which allows the use of the simplified equations. The resulting matrix equation to be solved is (Note: . means zero):

┌		┐ ┌       ┐          ┌                        ┐      ┌      ┐
│2  1  .  1   2	│ │x1" y1"│          │x2-x1-x5+x4  y2-y1-y5+y4│      │-3  -1│
│1  4  1  .   .	│ │x2" y2"│ = 6/(h*h)│  x3-2x2+x1    y3-2y2+y1│      │ 3   1│
│.  1  4  1   .	│ │x3" y3"│          │  x4-2x3+x2    y4-xy3+y2│ = 96 │-1  -3│
│.  .  1  4   1	│ │x4" y4"│          │  x5-2x4+x3    y5-2y4+y3│      │ 1   3│
│1  .  .  .  -1	│ │x5" y5"│          │          0            0│      │ 0   0│
└   		┘ └       ┘          └                        ┘      └      ┘

The solution is x" = [-120, 120, -72, 72, -120]^T and y" = [-72, 72, -120, 120, -72]^T, from which we can calculate the coefficients of the spline segments and generate the plot.
.

[VBAhack]

Example Problem #3 (Bezier curve)

As a last example, we will show how to relate a cubic Bezier representation to a parametric cubic spline formulated using the methodology shown here. A cubic Bezier curve is represented by 4 control points, 1 each to define the beginning and end of the curve, and 2 additional control points, not necessarily on the curve, to define the boundary conditions at either end. A common representation of a cubic Bezier (http://en.wikipedia.org/wiki/B%C3%A9zier_curve) is:

B(t) = (1-t)³P₀ + 3t(1-t)²P₁ + 3t²(1-t)P₂ + t³P₃..........(0 ≤ t ≤ 1)

which is just a compact form of representing 2 equations:

x(t) = (1-t)³x₀ + 3t(1-t)²x₁ + 3t²(1-t)x₂ + t³x₃..........(0 ≤ t ≤ 1)
y(t) = (1-t)³y₀ + 3t(1-t)²y₁ + 3t²(1-t)y₂ + t³y₃..........(0 ≤ t ≤ 1)

P₀ = (x₀,y₀) and P₃ = (x₃,y₃) are control points defining the endpoints of the Bezier curve and P₁ = (x₁,y₁) , P₂ = (x₂,y₂) are additional control points the define the shape of the curve. With the above equations it is easy to verify that x = x₀, y = y₀ at t = 0 and x = x₃, y = y₃ at t = 1. By expanding the above equations and collecting terms, it can be shown that a cubic Bezier curve can be also represented in traditional polynomial form:

x(t) = a_xt³ + b_xt² + c_xt + d_x..........(0 ≤ t ≤ 1)
y(t) = a_yt³ + b_yt² + c_yt + d_y..........(0 ≤ t ≤ 1)

where the coefficients of the polynomials are related to the Bezier control points by the following matrix equation (Note: . means zero):

┌  	┐   ┌              ┐ ┌      ┐  
│ax  ay	│   │-1   3  -3   1│ │x0  y0│  
│bx  by	│ = │ 3  -6   3   .│ │x1  y1│ 
│cx  cy	│   │-3   3   .   .│ │x2  y2│  
│dx  dy	│   │ 1   .   .   .│ │x3  y3│  
└  	┘   └              ┘ └      ┘

This is essentially a transformation matrix that allows conversion between 2 different representations of a parametric cubic polynomial. Let’s find the Bezier control points for example problem #1. Before proceeding, however, there is a subtle but important point that needs to be discussed. The above relations are valid for 0 ≤ t ≤ 1. In terms of the cubic spline formulation shown here, by comparing equations it is clear that t corresponds to the quantity (s-s₀). Thus, in order to convert a parametric cubic polynomial to traditional Bezier form, the formulation must be done in such a way that (s-s₀) varies from 0 to 1 in each segment. Actually, this is easy. All that is needed is to increment s by 1 for each point and not scale the result. Thus, s = [0, 1, 2, 3]. This method of choosing s has the additional benefit of h₁ = h₂ = h₃ = h = 1, which allows the use of the simplified version of the cubic spline equation. Choosing values for free end boundary conditions (for a pleasing result), we get:

┌		┐ ┌       ┐     ┌                      ┐   ┌         ┐
│1   .   .   .	│ │x1" y1"│     │        1          0.6│   │   6  3.6│
│1   4   1   .	│ │x2" y2"│ = 6 │x3-2x2+x1    y3-2y2+y1│   │  18    6│
│.   1   4   1 	│ │x3" y3"│     │x4-2x3+x2    y4-xy3+y2│ = │  -6  -18│
│.   .   .   1	│ │x4" y4"│     │     -0.6           -1│   │-3.6   -6│
└   		┘ └       ┘     └                      ┘   └         ┘

The solution is x" = [6, 3.36, -1.44, -3.6]^T and y" = [3.6, 1.44, -3.36, -6]^T. From the cubic spline tutorial, we can conveniently represent the equations to determine the spline segment coefficients by the following matrix equation (example shown for 1st segment).

┌  	┐   ┌                 ┐ ┌        ┐  
│ax  ay	│   │ .  .  -1/6   1/6│ │x1   y1 │  
│bx  by	│ = │ .  .   1/2     .│ │x2   y2 │ 
│cx  cy	│   │-1  1  -1/3  -1/6│ │x1"  y1"│  
│dx  dy	│   │ 1  .     .     .│ │x2"  y2"│  
└  	┘   └                 ┘ └        ┘

Since all x_i, y_i, x_i", y_i" are known, the spline coefficients can be readily calculated. Once the spline coefficients are determined, the associated Bezier control points can be calculated using the relation previously shown.

The polynomial coefficients are thus:

┌  				┐   ┌                 ┐ ┌                            ┐  
│ax1   ay1  ax2  ay2  ax3   ay3	│   │ .  .  -1/6   1/6│ │x1   y1   x2   y2   x3   y3 │  
│bx1   by1  bx2  by2  bx3   by3	│ = │ .  .   1/2     .│ │x2   y2   x3   y3   x4   y4 │ 
│cx1   cy1  cx2  cy2  cx3   cy3	│   │-1  1  -1/3  -1/6│ │x1"  y1"  x2"  y2"  x3"  y3"│  
│dx1   dy1  dx2  dy2  dx3   dy3	│   │ 1  .     .     .│ │x2"  y2"  x3"  y3"  x4"  y4"│  
└  				┘   └                 ┘ └                            ┘  
  ┌                 	┐ ┌                                      ┐  
  │ .   .   -1/6    1/6	│ │   1     0     -1      0      0      1│  
= │ .   .    1/2      .	│ │  -1     0      0      1      0     -1│ 
  │-1   1   -1/3   -1/6	│ │   6   3.6   3.36   1.44  -1.44  -3.36│  
  │ 1   .      .      .	│ │3.36  1.44  -1.44  -3.36   -3.6     -6│  
  └                 	┘ └                                      ┘ 
  ┌                                      ┐  
  │-0.44  -0.36  -0.8  -0.8  -0.36  -0.44│  
= │    3    1.8  1.68  0.72  -0.72  -1.68│ = polynomial coefficients
  │-4.56  -1.44  0.12  1.08   1.08   0.12│  
  │    1      0    -1     0      0      1│  
  └                                      ┘

And the Bezier control points are (note the matrix inverse):

┌  				┐   ┌              ┐-1 ┌                              ┐  
│x01   y01  x02  y02  x03   y03	│   │-1   3  -3   1│   │ax1   ay1  ax2  ay2  ax3   ay3│  
│x11   y11  x12  y12  x13   y13	│ = │ 3  -6   3   .│   │bx1   by1  bx2  by2  bx3   by3│ 
│x21   y21  x22  y22  x23   y23	│   │-3   3   .   .│   │cx1   cy1  cx2  cy2  cx3   cy3│  
│x31   y31  x32  y32  x33   y33	│   │ 1   .   .   .│   │dx1   dy1  dx2  dy2  dx3   dy3│  
└  				┘   └              ┘   └                              ┘ 
  ┌          	┐ ┌                                      ┐  
  │.   .   .  1	│ │-0.44  -0.36  -0.8  -0.8  -0.36  -0.44│  
= │.   . 1/3  1	│ │    3    1.8  1.68  0.72  -0.72  -1.68│ 
  │. 1/3 2/3  1	│ │-4.56  -1.44  0.12  1.08   1.08   0.12│  
  │1   1   1  1	│ │    1      0    -1     0      0      1│  
  └       	┘ └                                      ┘ 
  ┌                                     ┐  
  │    1      0     -1     0     0     1│  
= │-0.52  -0.48  -0.96  0.36  0.36  1.04│= Bezier control points
  │-1.04  -0.36  -0.36  0.96  0.48  0.52│  
  │   -1      0      0     1     0    -1│  
  └                                     ┘

Conclusion:

We have extended our knowledge of cubic splines to include the powerful and flexible parametric version and have fit various shapes (including a closed loop) to a set of given points. We have also shown how to determine the cubic Bezier control points of a curve.

.

[frankly747]

In the base matrix, what do the question marks stand for ?

[frankly747]

In the base matrix, what do the question marks stand for ?

[VBAhack]

In the base matrix, what do the question marks stand for ?

It just means the boundary conditions at each end of the spline need to be chosen.

[frankly747]

It just means the boundary conditions at each end of the spline need to be chosen.

Hi. I've tried doing some camull-rom interpolation for the following data.
x y
0.2 1.221402861
0.3 1.349864712
0.4 1.491929555
0.5 1.649697833
0.6 1.828165418
0.7 2.042000232
0.8 2.332915111
0.9 2.808281551
1 3.718281828
1.1 5.597908484
1.2 9.511853345
1.3 17.45514585
1.4 32.98066546
1.5 62.14672813
1.6 114.9041952
1.7 207.0733374
1.8 363.0963701
1.9 619.7925202
2 1031.389056

I chose to interpolate for point 1.1 (i'm assuming this point is non-existent).
I am using the following method:

f(x) = [1, x, x^2, x^3] * M * [p0, p1, p2, p3]

And M is the matrix,

0.0 1.0 0.0 1.0
-0.5 0.0 0.5 0.0
1.0 -2.5 2.0 -0.5
-0.5 1.5 -1.5 0.5

The "tension" value (greek letter tau) is 0.5.

Thus, the 4 points chosen were 0.9, 1, (2 before 1.1), 1.2, 1.3 (the 2 after 1.1).

I did it in matlab and got wrong values (using an input of 1.1 and expecting to get a value close to 5.5979). I have since learnt that the input needs to be between 0 and 1. How do I scale this value of 1.1 to be between 0 and 1 ? And does anything else change , in the M matrix, matrix of points [p0, p1, p2, p3] ?

Thankyou.

[VBAhack]

Aeroengie,

You need to include both x & y in the matrix of p's. It will be two columns. Should work.

[frankly747]

Aeroengie,

You need to include both x & y in the matrix of p's. It will be two columns. Should work.

Hi. Thanks for replying. But, if the P matrices are two columns, then the final answer will be a 1 X 2, that is, one row two columns

( 1x4 * 4X4 * 4X2 ), where 1X4 is the matrix of the interpolant (cubed, squared, raised to one and constant), 4X4 is the 'M' matrix and the 4X2 is the matrix you are proposong).

Thus, if I have a 1 X 2, this doesn not sound right since the answer I seek is a single value (the value of y, given an x).

Holla.

[sshfl]

VBAhack, thanks for your tutorials on cubic splines, they are very informative and have been a great help. I was wondering if you have had any experience with constraining the length of the spline segments. As a simple example, lets say I have 5 equally spaced points on a string that is swaying and moving around in three dimensions on a windy day. Since I know the length of the string, I would like to constrain the length of the 4 spline segments so that each spline segment is 25% of the total length of the string. I was thinking that maybe moving up to a quartic spline would add another degree of freedom that would allow me to also constrain the length of each spline segment, but the arc length integral cannot be computed analytically for these functions. Any ideas or suggestions would be greatly appreciated!

Thanks

[jemidiah]

I'm sorry, I don't think VBAHack visits these forums anymore. This thread was made in 2007 and the last response before yours was in January of 2010.

The arc length integral of a quartic spline should be very susceptible to numeric integration, since in any interval its derivatives would be bounded in magnitude. There would still be problems, since the resulting condition would be (highly) nonlinear. The only way I can even think of that approach working out is through a sequence of progressive approximations using quadrature on the parameter space--which is horribly messy to say the least, and would probably be very computationally intensive. There would probably be uniqueness issues, too, where multiple solutions exist.

My first thought is to replace the polynomial interpolation with another type of interpolation, where the arc length condition becomes a linear one. I'm not sure what type of interpolation would have this property; using sum_n a_n cos(nt) for n=0, ..., k comes close at first mental glance. Perhaps some modification would work.

My second thought is to break up your line into a large number of stiff segments and apply Lagrangian physics to them with energy functions tuned to avoid kinks. Some of the segments could end in fixed points. You would then basically run a physics simulation and wait until variation nearly stops. This method pretty much simulates a rubber band of fixed length passing through a few fixed points.

Good luck.