30th January, 2011

Cubic spline interpolation

The following may make more sense when the latex-mathjax plugin is enabled. You’ll see the equations then, not the LaTeX source.

I’ve got a set of $(x,y)$ points. What I want to do is draw a reasonably smooth curve through all of them, rather than a set of straight lines connecting them. To do so, I need some higher order function such as a cubic. In particular the cubic Bezier function is a well known graphics primitive. Therefore, given a set of points and the definition of a cubic Bezier, I need to calculate the control points for each curve.

A cubic Bezier function interpolates between $K_0$ and $K_1$ using control points $C_{0,1}$ and $C_{0,2}$ in the following manner
$$
B(t)=(1-t)^3K_0+3t(1-t)^2C_{0,1}+3t^2(1-t)C_{0,2}+t^3K_1.
$$
We may extend the idea to a spline. The knots of the spline are $langle K_0, K_1,ldots,K_nrangle$. This leads to $n-1$ cubic functions $B_0, B_1, ldots, B_{n-1}$. Where $B_i$ interpolates between $K_i$ and $K_{i+1}$ with control points $C_{i,1}$ and $C_{i,2}$.

We want $C^2$ continuity so that our curve is nice and smooth. Therefore, for each $i$ we require
begin{eqnarray*}
B^prime_i(1)&=&B^prime_{i+1}(0)textrm{ and}\
B^{primeprime}_i(1)&=&B^{primeprime}_{i+1}(0).
end{eqnarray*}
Furthermore,
$$
B_i^prime(t)=3K_{i+1}t^2-3C_{i,2}t^2+6C_{i,2}left(1-tright)t-6C_{i,1}left(1-tright)t+3C_{i,1}left(1-tright)^2-3K_ileft(1-tright)^2
$$
and
$$
B_i^{primeprime}(t)=6 K_i (1 – t) + 6 C_{i,1} t – 12 C_{i,1} (1 – t) – 12 C_{i,2} t + 6 C_{i,2} (1 – t) + 6 K_{i+1} t.
$$

Evaluating $B^prime_i(1)$ gives
begin{eqnarray*}
B_i^prime(1)&=&3K_{i+1}1^2-3C_{i,2}1^2+6C_{i,2}left(1-1right)1-6C_{i,1}left(1-1right)1+3C_{i,1}left(1-1right)^2-3K_ileft(1-1right)^2\
&=&3K_{i+1}-3C_{i,2}+6C_{i,2}(0)1-6C_{i,1}(0)1+3C_{i,1}(0)^2-3K_i(0)^2\
&=&3K_{i+1}-3C_{i,2}.
end{eqnarray*}
and evaluating $B_{i+1}^prime(0)$ gives
begin{eqnarray*}
B_{i+1}^prime(0)&=&3K_{i+2}0^2-3C_{i+1,2}0^2+6C_{i+1,2}left(1-0right)0\
&=&-6C_{i+1,1}left(1-0right)0+3C_{i+1,1}left(1-0right)^2-3K_{i+1}left(1-0right)^2\
&=&3C_{i+1,1}-3K_{i+1}
end{eqnarray*}
Therefore, for each $i$
begin{eqnarray*}
B_i^prime(1)&=&B_{i+1}^prime(0)\
3K_{i+1}-3C_{i,2}&=&3C_{i+1,1}-3K_{i+1}\
6K_{i+1}&=&3C_{i,2}+3C_{i+1,1}\
K_{i+1}&=&0.5C_{i,2}+0.5C_{i+1,1}.
end{eqnarray*}

Evaluating $B_i^{primeprime}(1)$ gives
begin{eqnarray*}
B_i^{primeprime}(1)&=&6 K_i (1 – t) + 6 C_{i,1} t – 12 C_{i,1} (1 – t) – 12 C_{i,2} t + 6 C_{i,2} (1 – t) + 6 K_{i+1} t\
&=&6K_i(0)+6C_{i,1}-12C_{i,1}(0)-12C_{i,2}+6C_{i,2}(0)+6K_{i+1}\
&=&6C_{i,1}-12C_{i,2}+6K_{i+1}\
end{eqnarray*}
and evaluating $B_{i+1}^{primeprime}(0)$ gives
begin{eqnarray*}
B_{i+1}^{primeprime}(1)&=&6 K_{i+1} (1 – t) + 6 C_{i+1,1} t – 12 C_{i+1,1} (1 – t) – 12 C_{i+1,2} t + 6 C_{i+1,2} (1 – t) + 6 K_{i+2} t\
&=&6 K_{i+1} (1 – 0) + 6 C_{i+1,1} 0 – 12 C_{i+1,1} (1 – 0) – 12 C_{i+1,2} 0 + 6 C_{i+1,2} (1 – 0) + 6 K_{i+2} 0\
&=&6 K_{i+1} – 12 C_{i+1,1} + 6 C_{i+1,2}.
end{eqnarray*}
Therefore, for each $i$
begin{eqnarray*}
B_i^{primeprime}(1)&=&B_{i+1}^{primeprime}(0)\
6C_{i,1}-12C_{i,2}+6K_{i+1}&=&6 K_{i+1} – 12 C_{i+1,1} + 6 C_{i+1,2}\
0&=&-6C_{i,1}+12C_{i,2}-12 C_{i+1,1} + 6 C_{i+1,2}\
0&=&-C_{i,1}+2C_{i,2}-2 C_{i+1,1} + C_{i+1,2}\
end{eqnarray*}

We now have $(2times n-1)-2$ equations and $(2times n-1)$ unknown control points. We construct the linear system
$$
left(
begin{array}{ccccccc}
0&0.5&0.5&0&0&0&ldots\
-1&2&-2&1&0&0&ldots\
0&0&0&0.5&0.5&0&ldots\
0&0&-1&2&-2&1&ldots\
%0&0&0&0.5&0.5&0&ldots\
%0&0&1&-2&2&-1&ldots\
vdots&&&&&&
end{array}
right)
times
left(
begin{array}{c}
C_{0, 1}\
C_{0, 2}\
C_{1, 1}\
C_{1, 2}\
ldots
end{array}
right)
=
left(
begin{array}{c}
K_1\
0\
K_2\
0\
ldots
end{array}
right).
$$

We also have the natural boundary conditions that $B^{primeprime}_0(0)=0$ and $B_{n-1}^{primeprime}(1)=0$. The linear system then becomes
$$
left(
begin{array}{ccccccccc}
2&-1&0&0&0&0&ldots&0&0\
0&0.5&0.5&0&0&0&ldots&0&0\
-1&2&-2&1&0&0&ldots&0&0\
0&0&0&0.5&0.5&0&ldots&0&0\
0&0&-1&2&-2&1&ldots&0&0\
%0&0&0&0.5&0.5&0&ldots&0&0\
%0&0&1&-2&2&-1&ldots&0&0\
vdots&&&&&&&&\
0&0&0&0&0&0&ldots&-1&2
end{array}
right)
times
left(
begin{array}{c}
C_{0, 1}\
C_{0, 2}\
C_{1, 1}\
C_{1, 2}\
ldots
end{array}
right)
=
left(
begin{array}{c}
0\
K_1\
0\
K_2\
0\
ldots\
0
end{array}
right).
$$

Posted at 6:39 pm | Comments Off

7th January, 2011

More hills near Petra

What can I say? I like hills.

Posted at 12:00 pm | Comments Off

5th January, 2011

“The Treasury” III (Son of Treasury)

Yet another picture of the treasury.

Posted at 12:00 pm | Comments Off

3rd January, 2011

Tombs in Petra

There were so many tombs like this carved into the hillside in Petra that after a while you almost didn’t notice them and started to look instead for signs of the hidden alien spacecraft that brought the Nabataeans to our planet.

Posted at 12:00 pm | Comments Off