Curve Expedition
The world of mathematical curves is rich and vast, with scenic mountain tops, and refreshing meadows and brooks, but also with hazardous swamps and deserts that will ensnare the naïve traveler. In these notes we navigate through this world to the B-Spline curves, which are fundamental in many ways and which have all sorts of beautiful properties. We'll learn what these curves are and what makes them special.
B-Spline curves are relatively inaccessible, and so most people only ever see them from afar. You can see B-Splines from 30,000 feet, for instance by using them in a canned CAD application, however it's an entirely different and immeasurably richer experience to visit them on foot- to get into their roots and to understand their many details.
Our trek through the world of curves begins with the following preliminaries.
- What Is A Curve?
- Unit Vectors And Curvature
- Curvature Comb
- Optical Reflections
- General Parameterization
- Motion Along A Curve
- Continuity And Smoothness
- Order And Degree
At this point we've prepared our hiking gear- our boots, our food, our compass, etc. On our first actual foray into the world itself, we encounter cubic Bezier curves, which are ubiquitous in practical graphics work. If you've drawn a curve in a program like Adobe Illustrator, chances are it was a cubic Bezier curve.
The Bezier construction process is our foundation, and from this we build a new thing which is more general and more powerful.
We're deep into our trek now, and in a position to make some important observations.
We end with some of the more esoteric notes from our excursion, as well as references that other explorers may find helpful.
What Is A Curve?
When we talk about a curve, we mean a path traced out in the plane. Many of the ideas we explore extend easily to paths traced out in more exotic spaces, however we stick to the plane to minimize our bookkeeping. In addition, as much as we love space filling curves and fractals, we restrict our focus to curves that might reasonably be made into physical objects, such as the edges of a computer.
Workers in design who create the physical objects around us usually do so with curves that have components that are polynomial functions of some parameter. For instance, $$\mathbf{x}(t) = \left[\! \begin{array}{c}t\notag\\ 2t-t^2\end{array} \!\right]$$
As \(t\) increases from \(0\) to \(2\), this curve traces out the arc shown below in blue. Appending another segment to this curve makes it piecewise polynomial. We do this by defining $$\mathbf{x}(t)= \left\{ \begin{array}{l} \left[\! \begin{array}{c}t\notag\\2t+t^2\end{array} \!\right] \text{ when } t\in[-2,0)\\ \\ \left[\! \begin{array}{c}t\notag\\2t-t^2\end{array} \!\right] \text{ when } t\in[0,2) \end{array} \right.$$
Here is what this curve looks like, with pieces illustrated in red and blue.
When \(||d\mathbf{x}/dt||=1\), we say that \(t\) is an arc-length parameter. We typically denote an arc-length parameter by the letter \(s\).
Unit Vectors And Curvature
If \(\mathbf{x}(s)\) is a parametric curve, with \(s\) an arc-length parameter, then the unit tangent vector \(\mathbf{t}\) to the curve is given by $$\mathbf{t}=\frac{d\mathbf{x}}{ds}.$$
Restricting our attention to planar curves, we define the inclination angle as the angle \(\theta\) that \(\mathbf{t}\) makes with the horizontal. We also define \(\mathbf{t}^{\perp}\) as the \(\theta=90^{\circ}\) rotation of \(\mathbf{t}\).
We define the curvature \(\kappa\) as
$$\kappa=\frac{d\theta}{ds}$$As we move along a curve at constant speed, the curvature is simply the speed of rotation of the forward pointing unit vector \(\mathbf{t}\).
Another way of thinking about the curvature \(\kappa\) at a point \(\mathbf{x}(s)\) comes from the construction of a circle. Draw two lines perpendicular to the curve, one passing through \(\mathbf{x}(s)\), and the other passing through the nearby point \(\mathbf{x}(s+\epsilon)\). As \(\epsilon\) shrinks to zero, the point where these lines intersect approaches a point \(\mathbf{a}\). Our circle is centered at \(\mathbf{a}\) and passes through \(\mathbf{x}(s)\). We call it the osculating circle at \(\mathbf{x}(s)\) because it locally matches the curve so well there (the term derives from the Latin osculari which means to kiss). The radius of the osculating circle is equal to \(1/|\kappa(s)|\). The smaller the osculating circle, the larger the curvature.
Curvature Comb
A curvature comb is additional artwork which makes it easy to understand the curvature of a curve. The curvature comb associated with the curve \(\mathbf{x}(t)\) is a new curve \(\mathbf{u}(t)\) given by $$ \mathbf{u}(t) = \mathbf{x}(t)-\gamma\kappa\mathbf{t}^{\perp}, $$ where \(\gamma\) is a scaling factor that is chosen case by case. As an example, consider the following curve \(\mathbf{x}(t)\).
Two of the corners of \(\mathbf{x}(t)\) are perfect circular arcs that transition abruptly to straight lines, while the other two corners are B-Spline curves (we discuss B-Splines later in these notes). Subtle differences such as these can have a profound effect on the appearance of an object, especially when optical reflections are involved, as we discuss in the next section.
Optical Reflections
Consider an object viewed by reflection off a shiny surface.
The viewing angle \(\theta\) depends on the curvature \(\kappa\) of the reflecting surface. In the given scene, \(\kappa = 0.5\) and \(\theta = 1.01^{\circ}\), marked by the purple dot in the graph. The viewing angle corresponds to the reflected image size, and so an abrupt change in curvature causes an abrupt change in the reflected image.
As an example, we extrude our squircle from the previous page, creating a solid object with a mirror finish. We surround it with dots and see how they are reflected.
As expected, the size of the dots (in the horizontal direction) depends on the surface curvature. The discontinuity in curvature on the right causes an abrupt change in the amount of horizontal squishing. It's even easier to see these effects when the object is in motion.
Turning a shiny object in your hands and looking at how it reflects your surroundings is a great way to assess its smoothness. Reflected light is an unforgiving judge of surface geometry!
General Parameterization
Curves are often parameterized by a non arc-length parameter \(t\). In terms of \(t\), the unit tangent vector is given by $$ \mathbf{t}=\frac{d\mathbf{x}}{ds}=\frac{d\mathbf{x}}{dt}\frac{dt}{ds}=\frac{\dot{\mathbf{x}}}{\sqrt{x_{11}}}, $$ where a superposed dot indicates differentiation with respect to \(t\), and where \(x_{11}=\dot{\mathbf{x}}\!\cdot\!\dot{\mathbf{x}}\). Noting that \(\kappa\mathbf{t}^{\perp}=d\mathbf{t}/ds\), we differentiate our expression for \(\mathbf{t}\) and obtain $$ \kappa x_{11} = \ddot{\mathbf{x}}\cdot\mathbf{t}^{\perp}. \tag{1} $$ Differentiating again, we obtain $$ \dot{\kappa}x_{11}+3\kappa x_{12}=\dddot{\mathbf{x}}\cdot\mathbf{t}^{\perp}, \tag{2} $$ where \(x_{12}=\dot{\mathbf{x}}\!\cdot\!\ddot{\mathbf{x}}\). We can use \((1)\) to compute the curvature of any parameterized curve, and we can use \((2)\) to compute the rate of change of curvature. We can keep going in this way, differentiating again to obtain
$$ \ddot{\kappa}x_{11}-5\dot{\kappa}x_{12}+3\kappa(x_{22}+x_{13})=(\dddot{\mathbf{x}}+\ddddot{\mathbf{x}})\cdot\mathbf{t}^{\perp}. $$Obtaining derivatives with respect to the arc length \(s\) requires a bit more work. We note that
$$ \frac{d\kappa}{ds} = \dot{\kappa}\frac{1}{\sqrt{x_{11}}}, $$and that
$$ \frac{d^2\kappa}{ds^2} = \ddot{\kappa}\frac{1}{x_{11}} - \dot{\kappa}\frac{x_{12}}{x^2_{11}}. $$
example
Consider the parabola \(f(t)=t^2\). This curve and its derivatives are given by
$$
\mathbf{x} = \left[\!\begin{array}{c}t\notag\\t^2\end{array}\!\right],\;\;
\dot{\mathbf{x}} = \left[\!\begin{array}{c}1\notag\\2t\end{array}\!\right],\;\;
\ddot{\mathbf{x}} = \left[\!\begin{array}{c}0\notag\\2\end{array}\!\right],\;\;
\dddot{\mathbf{x}} = \left[\!\begin{array}{c}0\notag\\0\end{array}\!\right].
$$
The inner products \(x_{11}\) and \(x_{12}\) are given by
$$
\begin{align}
x_{11}&=1+4t^2,\notag\\
x_{12}&=4t.\notag
\end{align}
$$
We use \(\dot{\mathbf{x}}\) and \(x_{11}\) to obtain
$$
\mathbf{t} = \dot{\mathbf{x}}\frac{1}{\sqrt{x_{11}}}=
\left[\!\begin{array}{c}1\notag\\2t\end{array}\!\right]\frac{1}{\sqrt{1+4t^2}},
$$
and so
$$
\mathbf{t}^{\perp} =
\left[\!\begin{array}{c}-2t\notag\\1\end{array}\!\right]\frac{1}{\sqrt{1+4t^2}}.
$$
We now use \((1)\) to obtain
$$
\kappa = \frac{2}{(1+4t^2)^{3/2}}.
$$
The following graph shows \(f(t)=t^2\) in blue and \(\kappa=\kappa(t)\) in orange.
Motion Along A Curve
As the parameter \(t\) progresses through its values, the point \(\mathbf{x}(t)\) traces out a curve.
Sometimes we care about motion along the curve. For instance, imagine being a passenger in a vehicle located by \(\mathbf{x}(t)\).
Other times, motion along the curve is not important. For instance if \(\mathbf{x}(t)\) is the section curve of a satellite antenna such as the Stanford Dish shown below, then the only thing that matters is the totality of points \(\mathbf{x}(t)\) considered all at once.
Whether or not we care about motion along the curve affects how we think about smoothness.
Continuity And Smoothness
Continuity allows for a systematic and rigorous treatment of smoothness. Different notions of continuity are appropriate depending on whether we care about motion along a curve, or just the overall shape of the curve.
Basic Continuity
A function \(f(x)\) is continuous if having input values \(x\) and \(x_0\) close together causes the corresponding output values \(f(x)\) and \(f(x_0)\) to be close together. We say that such a function is \(C^0\). If the first derivative of a function is continuous we say the function is \(C^1\), and if the \(n\)th derivative of a function is continuous we say the function is \(C^n\). These ideas underly everything that follows.
Non Zero Derivative
We stick to curves \(\mathbf{x}(t)\) that have a derivative \(\dot{\mathbf{x}}(t)\) that is never zero. This is because if \(\dot{\mathbf{x}}=\mathbf{0}\) we can see things like a sharp cusp in a curve that is \(C^{\infty}\) (i.e., infinitely smooth). For instance, consider \(\mathbf{x}(t) = [t^3\;t^2]\).
Smoothness
If motion along the curve is important, then our notion of smoothness simply depends on \(\mathbf{x}(t)\) and its derivatives. If the first \(n\) derivatives of \(\mathbf{x}(t)\) are continuous, the curve is said to be parameterically continuous to order \(n\), and this is abbreviated by saying that the curve is \(C^n\). The higher the value of \(n\), the smoother the curve.
If we don't care about motion along the curve, then we need to step back and consider other possible motions besides the one encoded in \(\mathbf{x}(t)\). We do this by looking at reparameterizations of \(\mathbf{x}(t)\). If we can find a reparameterization continuous to order \(n\) (previous paragraph), then \(\mathbf{x}(t)\) is said to be geometrically continuous to order \(n\), and this is abbreviated by saying that this curve is \(G^n\).
In other words, if a function \(f(t)\) exists such that \(\mathbf{x}(f(t))\) is \(C^n\), then \(\mathbf{x}(t)\) is \(G^n\).
Examples
As an illustration of the idea of geometric continuity, consider a road being surveyed by a bad driver. We record position \(\mathbf{x}\) as a function of time \(t\), and because the driver is bad, the motion corresponding to \(\mathbf{x}(t)\) is filled with uncomfortable stops and starts. The question is, would it be possible for a good driver to do any better? Or is there a sharp turn in the road that would be uncomfortable even if the driver was really good? If we can reparameterize \(\mathbf{x}(t)\) to get rid of all the stops and starts, it means the road is smooth. If we can't, it means the road has that sharp turn.
Provided we restrict ourselves to curves with a non-zero derivative, the initial levels of geometric continuity have tangible interpretations as follows:
- If the unit tangent vector \(\mathbf{t}(t)\) of a curve is continuous, then the curve is \(G^1\).
- If the curvature \(\kappa\) of a curve is continuous, then the curve is \(G^2\).
- If the curvature of a curve has no cusps- that is, if \(\dot{\kappa}\) is continuous, then the curve is \(G^3\).
The challenge with continuity is usually at the junctions of piecewise polynomial curves. Certain piecewise polynomial curves (e.g., cubic splines) are intuitive and easy to work with, however it can be difficult to achieve \(G^2\) let alone \(G^3\) continuity at their junction points. Other piecewise polynomial curves (e.g., high enough order B-Splines) are \(G^2\) and \(G^3\) automatically, however they are more involved to work with.
Order And Degree
The order of a polynomial is the number of coefficients, while the degree is the highest power. For instance, the following polynomial has order \(k+1\) and degree \(k\). $$ f(x) = a_0+a_1x+a_2x^2+\cdots+a_kx^k $$ Order is always one more than degree.
Bezier Curves
The Bezier construction process is at the foundation of much of our later work. The process is simple, and yet it generates beautiful curves that are easy to control.
A Bezier curve is parameterized by a chain of vectors \(\mathbf{u}_1\), ..., \(\mathbf{u}_n\) that in this context are referred to as control points. From this we build a new chain of one fewer points \(\mathbf{v}_1\), ..., \(\mathbf{v}_{n-1}\), with $$ \mathbf{v}_i=(1-t)\mathbf{u}_i+t\mathbf{u}_{i+1}, $$ where \(t\) is a parameter between \(0\) and \(1\). This step is repeated over and over, always using the same value of \(t\), each time producing a new shorter chain. Ultimately we get a chain of one point which traces out a curve as \(t\) goes from \(0\) to \(1\).
The curve traced by the final vector begins at the first control point \(\mathbf{u}_1\) (when \(t=0\)) and ends at the last control point \(\mathbf{u}_n\) (when \(t=1\)). The curve follows the control points in an intuitive way- manipulating the curve is done by moving around the control points.
When we call a curve a Bezier curve, we mean that the curve has been parameterized as described above. Any polynomial curve can be parameterized in this way, and so the terms Bezier curve and polynomial curve can be used interchangeably.
Blossoms & Subdivision
By making the interpolation values at each step of the Bezier process independent, we obtain a vector valued function known as a blossom.
In the above example, the blossom is the function \(\mathbf{X}(t_1,t_2,t_3)\). Of course, when the \(t_i\)'s are equal, we obtain the original cubic polynomial (aka Bezier) curve which we denote by \(\mathbf{x}(t)\). $$\mathbf{x}(t)= \mathbf{X}(t,t,t)$$ Every cubic polynomial \(\mathbf{x}:\mathbb{R}\longrightarrow V\) has a uniquely associated blossom \(\mathbf{X}:\mathbb{R}^3\longrightarrow V\). So far, we've worked with chains of control points defined with respect to the parameter interval \([0,1]\). Generally, working with respect to the parameter interval \([a,b]\), our construction of the blossom for a cubic polynomial is as follows.
The blossom \(\mathbf{X}\) is a weighted sum of the \(\mathbf{u}_i\)'s, with each weight consisting of the sum of products of one term from each row of the diagram. It follows that \(\mathbf{X}\) is linear in each of its arguments. Another property of \(\mathbf{X}\) comes from the observation that the following diagrams are equivalent (i.e., \(\mathbf{f}=\mathbf{g}\) for any collection of inputs \(\mathbf{u}_i\)).
From this it follows (immediately) that the outputs \(\mathbf{v}_i\) of a longer diagram such as the following are unaffected by swapping \(\alpha\) and \(\beta\).
And from this it is clear that the rows in any diagram like \((\bigstar)\) can be swapped without affecting the output. That is, \(\mathbf{X}(t_1,t_2,t_3)\) and generally \(\mathbf{X}(t_1,t_2,...,t_n)\) are unaffected by a reordering of their arguments.
We have seen that a chain of control points \(\mathbf{u}_i\) and an associated parameter interval \([a,b]\) define a polynomial (aka Bezier) curve as well as its associated blossom. In the other direction, we can start with a blossom \(\mathbf{X}(t_1,t_2,...,t_n)\) and an interval \([a,b]\), and find the associated chain of control points. We get to these control points from the vertices of the hypercube \([a,b]^n\), as shown below in the case \(n=3\). Here, the function \(\mathbf{X}\) maps from its domain \(\mathbb{R}^3\) on the left, to the vector space \(V\) on the right.
Generally, given the blossom \(\mathbf{X}(t_1,t_2,...,t_n)\), the \(k\)th control point is obtained by setting \(k-1\) of the inputs equal to \(b\), and the rest equal to \(a\). We remember this fact as \(\clubsuit\). In fact, all the construction points that arise during a Bezier construction can be expressed by the blossom. For instance, for a cubic Bezier curve:
Subdivision
Consider the edge of the above diagram circled in green, and recall \(\clubsuit\). These points are the control points for this blossom with respect to the interval \([a,t]\). That is, the points that arise during the Bezier construction of \(\mathbf{X}(t,t,t)=\mathbf{x}(t)\) are the control points associated with \([a,t]\). This process, of finding an alternative set of control points corresponding to a different interval for the same polynomial is called subdivision. Each interval \([a,b]\) corresponds to a hypercube in the blossom domain, the vertices of which map to the chain of control points corresponding to the interval.
We have seen how to construct control points from a blossom. It is also straightforward to construct these points from the corresponding polynomial \(\mathbf{x}(t)\). This involves working with the derivatives of \(\mathbf{x}(t)\). For instance, if \(\mathbf{x}(t)\) is a cubic polynomial, then the Bezier construction diagram over the interval \([a,b]\) is as follows:
We find that
| \(\mathbf{u}_1 = \mathbf{x}(a)\) |
|---|
| \(\mathbf{u}_2 = \mathbf{x}(a)+\frac{b-a}{3}\dot{\mathbf{x}}(a)\) |
| \(\mathbf{u}_3 = \mathbf{x}(b)-\frac{b-a}{3}\dot{\mathbf{x}}(b)\) |
| \(\mathbf{u}_4 = \mathbf{x}(b)\) |
This can be generalized to larger diagrams, and there's also lots to explore on the topic of blossom differentiation. For additional material on blossoms I recommend Ron Goldman's notes.
As an alternative to working with blossoms, subdivision can be established purely in terms of the construction diagram for the polynomial curve.
The key to this is the magical equivalence between the following two diagrams.
This equivalence is obvious when there are only two inputs, corresponding to linear interpolation between two points. Moving some fraction \(\alpha\) of the distance from one point \(u_1\) to another point \(u_2\) gives us a new point \(u_3\). Moving some fraction \(\beta\) of the distance from \(u_1\) to \(u_3\) gives us a final point \(u_4\), which is the fraction \(\alpha\beta\) along the line from \(u_1\) to \(u_2\). In terms of our diagram, this base case is represented as follows.
We prove the general case by induction. Assuming the equivalence holds for a certain number of inputs say \(k\), we feed our base case into each of these, obtaining a diagram with \(k+1\) inputs.
From our assumption and our base case, we immediately see that the above two diagrams are equal. It only remains to transform the diagram on the left into the proper form. We do this by dividing it into two cases.
Every path through the left diagram encounters a single \ before passing through the main network of blue lines. There is no change to having this encounter occur after the main networkk of blue lines as shown below. We reconfigure the right diagram similarly.
Next we add some channels and a new input to the left diagram, which have no effect whatsoever on the output value.
We combine the left and right diagrams together to obtain what we were after, establishing the induction step.
Cubic Curves
We do so much work with cubic curves that it is worth examining them in detail. Here is the Bezier construction process for building a cubic curve.
In addition to the Bezier control points \(\mathbf{u}_i\), we find it helpful to work with an alternative set of parameters. These are the end points \(\mathbf{a}=\mathbf{u}_1\) and \(\mathbf{b}=\mathbf{u}_4\), the unit tangent vectors \(\mathbf{e}_a\) and \(\mathbf{e}_b\) at the end points, and scalars \(\alpha\) and \(\beta\) which satisfy \(\alpha\mathbf{e}_a = \mathbf{u}_2-\mathbf{u}_1\) and \(\beta\mathbf{e}_b = \mathbf{u}_4-\mathbf{u}_3\).
In terms of these alternative parameters, the curve can be written as $$ \mathbf{x}(t) = \mathbf{a}f(t)+\mathbf{b}g(t)+3\alpha\mathbf{e}_au(t)+3\beta\mathbf{e}_bv(t), $$ where \(f(t)\), \(g(t)\), \(u(t)\), and \(v(t)\) are a basis for the cubic polynomials, defined as follows by their values and derivative values at \(t=0\) and \(t=1\).
| \(f(0)=1\), | \(g(0)=0\), | \(u(0)=0\), | \(v(0)=0\) |
|---|---|---|---|
| \(f(1)=0\), | \(g(1)=1\), | \(u(1)=0\), | \(v(1)=0\) |
| \(f'(0)=0\), | \(g'(0)=0\), | \(u'(0)=1\), | \(v'(0)=0\) |
| \(f'(1)=0\), | \(g'(1)=0\), | \(u'(1)=0\), | \(v'(1)=1\) |
The functions that satisfy these equations are as follows.
These are simple functions that are easy to work with.
Cubic Curvature
There is a straightforward expression for the curvature at the endpoints of a cubic curve. If the Bezier control points defining the curve are \(\mathbf{u}_1\), \(\mathbf{u}_2\), \(\mathbf{u}_3\), and \(\mathbf{u}_4\) as shown below, with \(\alpha\) the distance from \(\mathbf{u}_1\) to \(\mathbf{u}_2\), and with \(h\) the displacement from the line \(l\) through \(\mathbf{u}_1\) and \(\mathbf{u}_2\) to the point \(\mathbf{u}_3\), then the curvature at \(\mathbf{u}_1\) is given by $$ \kappa_1=\frac{2h}{3\alpha^2}. $$
As \(\mathbf{u}_3\) crosses from one side of \(l\) to the other, the sign of the displacement \(h\) changes, and this causes the sign of \(\kappa_1\) to flip. If \(\alpha\) is finite, the only way to get zero curvature at \(\mathbf{u}_1\) is to have \(\mathbf{u}_1\), \(\mathbf{u}_2\), and \(\mathbf{u}_3\) be co-linear.
For two cubic curves to have matching curvature where they meet, the following condition must hold.
Having aligned control arms of equal length (i.e., \(\alpha=\beta\) ) is necessary for \(C^1\) continuity at a junction point, but is unrelated to curvature continuity. In the above image for instance, we have curvature continuity at the junction point even though \(\alpha\neq\beta\).
B-Spline Curves
With only a small refinement, the Bezier construction process yields a whole family of so-called B-Spline curves. Each B-Spline curve is a chain of polynomial segments. As an example, from a list of 10 input points, we can obtain:
The original data is shown as small black dots connected by gray lines, while the output B-Spline curve segments are colored blue and light orange. The segments in the first B-Spline curve are discrete points that coincide with the original data. Successive B-Spline curves have fewer segments but these are higher order. The last curve in the progression is simply the Bezier curve associated with the original data, and so B-Spline curves include Bezier curves.
We can have \(C^{k-2}\) continuity at the junctions between order \(k\) segments. As an example, the cubic B-Spline shown above (with order 4 segments) is \(C^2\) at its junction points.
B-Spline curves with lower order segments are attractive for several reasons. Unlike a single Bezier curve through the same data, a low order B-Spline curve conforms closely to the data. Also, the segments of a B-Spline curve are only affected by nearby control points. This allows us to make adjustments to some segments without affecting others, and allows us to extend the original data with any number of new points without affecting anything except the last segment of the curve that has already been constructed. Lower order curves are numerically preferable as well.
B-Spline Construction
B-Spline curves are built by repeated linear combinations. We keep track of these using the dot diagram shown below. Dots in the top row represent the original data. Neighbors from this row are linearly combined to create the dots in the next row, and so on. The B-Spline curve is represented by the dots in the final row.
Each dot represents a vector valued function, mapping from an interval on the real line to a vector space (in our case the plane). The dots in the top row correspond to functions that are constant, simply equal to one of the original input data points. The dots in the bottom row correspond to polynomial functions which comprise the output B-Spline. We start off knowing the functions in the top row, and from these we construct the functions in the lower rows.
Intervals And Knots
The domains of the functions we build are intervals on the real line. Each dot in our diagram has an associated interval.
We specify the intervals for dots in the bottom row with a list of end points \(t_i\) that are traditionally known as knots. The intervals for dots in the next row up are the union of the intervals of their children.
Continuing in this way, we construct all the remaining intervals.
Interpolation
In the previous section we constructed the domain for each dot, going from the bottom of the dot diagram to the top. We now construct the function for each dot, going from top to bottom. The following figure shows the domains for functions \(\mathbf{u}_1\) and \(\mathbf{u}_2\), and their child function \(\mathbf{u}_3\). Note that \(\alpha_2=\alpha_3\) and \(\beta_1=\beta_3\).
Given \(\mathbf{u}_1\!:\![\alpha_1,\beta_1]\longrightarrow\mathbb{R}^2\) and \(\mathbf{u}_2\!:\![\alpha_2,\beta_2]\longrightarrow\mathbb{R}^2\), we construct the new function \(\mathbf{u}_3\!:\![\alpha_3,\beta_3]\longrightarrow\mathbb{R}^2\) according to $$ \mathbf{u}_3(t) = \mathbf{u}_1(t)+\frac{t-\alpha_3}{\beta_3-\alpha_3}(\mathbf{u}_2(t)-\mathbf{u}_1(t)). $$ Here's an example of functions generated by this process, starting with six input data points in the first row, and ending with three order 4 segments in the last row.
Segments \(\mathbf{u}_i\) in the final row comprise a piecewise polynomial function \(\mathbf{x}\), which is the output B-Spline. The domains of these \(\mathbf{u}_i\)'s are adjacent intervals \(U_i\) on the real line, and so we define \(\mathbf{x}\) over the union of the \(U_i\)'s by \(\mathbf{x}(t)=\mathbf{u}_i(t)\), where \(i\) is the index of the \(U_i\) that contains \(t\). In the appendix, we show that \(\mathbf{x}\) can be \(C^{n-2}\), where \(n\) is the number of rows.
Thoughts On Notation
A notation is a convention for conveying information by positioning symbols in two dimensional space (e.g., on a page or a screen). As an example, the notation \(a\mathbf{x}\) usually represents the product of a scalar \(a\) and a vector \(\mathbf{x}\). Consider the following alternative notation for this product, in which the action of scalar multiplication is represented by an arrow labeled with the scalar.
This is subtly different than the common maps to notion \(\mathbf{x}\mapsto\alpha \mathbf{x}\). Our labeled arrow actually represents the multiplication operation. The arrow transforms the vector \(\mathbf{x}\) assigned to the base into the vector \(\alpha\mathbf{x}\) assigned to the tip. Developing this further, let's agree that arrows coming together at a point have their output values summed and assigned to the point, and that arrows leading away from a point all use the point as the same base. In this way, by considering labeled arrows as a proper notation, we turn our dot diagram for B-Spline construction into a formal representation of the sums and products needed to build a B-Spline. As an example, consider the following equation from the Interpolation Section. $$ \mathbf{u}_3(t) = \mathbf{u}_1(t)+\frac{t-\alpha_3}{\beta_3-\alpha_3}(\mathbf{u}_2(t)-\mathbf{u}_1(t)) $$ In terms of our alternative notation, this expression can be represented as
Influence Of A Control Point
The influence of a control point \(\mathbf{x}_i\) on a B-Spline segment \(\mathbf{u}_j(t)\) arises from the paths that can be drawn between the dots corresponding to \(\mathbf{x}_i\) and \(\mathbf{u}_j(t)\).
Taking products along the paths, and then summing up the outputs, we find that \(\mathbf{x}_i\) contributes to \(\mathbf{u}_j(t)\) as \(\mathbf{x}_iB_{ij}(t)\), where \(B_{ij}(t)\) is a polynomial in \(t\) that depends on the knot values. If there are no paths between \(\mathbf{x}_i\) and \(\mathbf{u}_j(t)\), then \(B_{ij}(t)=0\). The segment \(\mathbf{u}_j(t)\) consists entirely of similar contributions from all of the input points, and so, if there are \(m\) input points, $$ \mathbf{u}_j(t) = \sum_{i=1}^m\mathbf{x}_iB_{ij}(t). $$ As long as it isn't identically zero, \(B_{ij}(t)\) is an order \(n\) polynomial in \(t\), where \(n\) is the number of rows in the dot diagram. We keep track of \(n\) by appending it as a superscript, writing \(B^n_{ij}(t)\). The following (distorted and truncated) dot diagram represents all the paths contributing to \(B_{ij}^n\) and \(B_{i+1,j}^n\) as clusters of red and blue lines respectively. We add a row to the top, and consider all the paths from the new dot marked \(i+1\) to \(\mathbf{u}_j(t)\).
This gives us the following recurrence relation. $$ B_{i+1,j}^{n+1}(t)=\frac{t-\alpha_i}{\beta_i-\alpha_i}B_{ij}^n(t)+\frac{\beta_{i+1}-t}{\beta_{i+1}-\alpha_{i+1}}B_{i+1,j}^n(t) $$ When the dot diagram consists of a single row, the B-Spline is piecewise constant, and the \(B_{ij}^1\)'s are either identically zero or one. Building up from these, we can use this recurrence relation to construct \(B_{ij}^n(t)\) for any \(n\).
Practical Construction
In the previous section we introduced the functions \(B_{ij}(t)\). The function \(B_{ij}(t)\) determines the influence of the control point \(\mathbf{x}_i\) on the jth segment \(\mathbf{u}_j(t)\) of the B-Spline.
$$ \mathbf{u}_j(t) = \sum_{i=1}^m\mathbf{x}_iB_{ij}(t) $$For each \(t\) let \(b_i(t)\) denote the output of the \(B_{ij}(t)\) which maps from the knot interval \([t_j,t_{j+1}]\) containing \(t\). The functions \(b_i(t)\) comprise a basis with respect to the given set of knots, and allow us to express the entire B-Spline \(\mathbf{x}(t)\) as
$$ \mathbf{x}(t) = \sum_{i=1}^m\mathbf{x}_ib_i(t) $$The following diagram shows the \(b_i\)'s and the grid of underlying \(B_{ij}\)'s, each represented by a dot. Unfilled dots (i.e., circles) represent the \(B_{ij}\)'s that are identically zero.
We have skewed the grid so that the non-trival \(B_{ij}\)'s are in the middle.
We refer to this array as \(B^k\) where \(k\) is the number of rows (and also the order of the polynomials \(B_{ij}\)). Our iterative construction of the \(B^k\)'s begins with \(k=1\), and in this base case all of the \(B_{ij}\)'s are simply equal to 1.
We increase the value of \(k\) as follows:
where \(U\) and \(V\) are arrays of second order (degree 1) polynomials in \(s\).
$$ U = \frac{(s+T)-L}{R-L},\;\;\;\;\; V = \frac{R-(s+T)}{R-L} $$The array \(L\) is the same size as \(B^k\), with each entry the left end point (knot value) of the support for the \(b_i\) passing through the entry.
Likewise, the array \(R\) is the same size as \(B^k\), with each entry the right end point (knot value) of the support for the \(b_i\) passing through the entry.
Finally, the array \(T\) is the same size as \(B^k\), with each entry the left end point (knot value) of the domain of the \(B_{ik}\) corresponding to the given entry.
Here is Matlab pseudo code for building these arrays.
t = [t1 t2 t3 ... tm];
L = repmat(t(1:m-1),k,1);
T = L;
R = repmat(t(2:m),k,1);
for i = 2:m-1
L(1:k-1,i) = L(2:k,i-1);
R(2:k,m-i) = R(1:k-1,m-i+1);
end
Polynomials on the interval \([t_i,t_{i+1}]\) are expressed in terms of the local variable \(s = t-t_i\). Working in terms of \(s\) (instead of \(t\)) keeps our coefficients from getting out of hand.
Knot Derivatives
Let \(f = f(x_1,x_2)\), and let \(g = g(y_1,y_2) = f(y_1+y_2,y_2)\). Note that
$$ \frac{\partial g}{\partial y_1} = \frac{\partial f}{\partial x_1},\;\;\;\;\;\; \frac{\partial g}{\partial y_2} = \frac{\partial f}{\partial x_1}+\frac{\partial f}{\partial x_2}, $$and so it follows that
$$ \frac{\partial f}{\partial x_2} = \frac{\partial g}{\partial y_2} - \frac{\partial g}{\partial y_1}. $$We care about this because in order to keep our coefficients from exploding, we express our B-Spline basis functions \(b_i\) in terms of local coordinates such as \(s=t-t_l\), where \(t_l\) is a knot value. That is, we express \(b_i(t,t_l)\) as \(b_i(t_l+s,t_l)\) which we call \(g_i(s,t_l)\). The array \(B^k\) from the previous section contains these \(g_i\)'s, and so the functions \(\frac{\partial}{\partial t_l}b_i(t,t_l)\) are given by
This seems simple now that I've written it up, however an error with these derivatives caused a subtle bug in some code I wrote that took me quite a while to track down.
A Basis For Each Segment
A segment from an order \(n\) B-Spline is the sum of \(n\) terms \(\mathbf{x}_iB_{ij}^n(t)\). The \(\mathbf{x}_i\)'s are from a vector space \(V\), while the scalar functions \(B_{ij}^n(t)\) arise by summing dot diagram paths from an input point to the given segment.
The \(n\) functions \(B_{ij}^n(t)\) that contribute to the given B-Spline segment are all order \(n\) polynomials in \(t\). Here we show that these functions comprise a basis for the order \(n\) polynomials in \(t\). We do this by using induction to show that these functions are independent.
Base Case
A single (nonzero) function is trivially independent.
Induction Step
A set of functions is independent if and only if whenever a linear combination of them is equal to zero, it follows that all the linear combination multipliers are zero. We suppose this is the case for the \(k\) functions associated with the possible paths from any set of \(k\) adjacent input points to the single descendant point that they and only they have in common. We call this single descendant point the apex of the \(k\) input points.
Pick any \(k+1\) adjacent input points, and for each of these choose a multiplier that will scale the functions emanating from the point. Let \(y_a(t)\) be the sum of the functions associated with the apex of the first \(k\) points, and let \(y_b(t)\) be the sum of the functions associated with the apex of the last \(k\) points.
The function \(f(t)\) associated with the apex of the \(k+1\) points is $$ f(t) = y_a(t)\cdot g(t)+y_b(t)\cdot(1-g(t)), $$ where \(g(t)\) is some order 2 (i.e., linear) function in \(t\). If \(f(t)=0\), then $$ (y_b(t)-y_a(t))\cdot g(t) = y_b(t). $$ Because \(g(t)\) is an order 2 polynomial, and \(y_a(t)\) and \(y_b(t)\) are order \(k\) (if they are not identically zero), the only way for this equation to be satisfied is for \(y_a(t)\) and \(y_b(t)\) to both be identically zero. It follows that the functions going from the selected \(k+1\) input points to their apex are independent as desired.
A Bigger Vector Space
The B-Splines generated by our dot diagram comprise their own vector space which we call \(S\). If \(V\) is the vector space of the dot diagram input points, then \(S\) is the set of piecewise polynomial functions from \(\mathbb{R}\) to \(V\), defined with respect to the knot values \(t_i\).
The vector space \(S\) is isomorphic to \(V\times \dots \times V\), which is great news because it is easy to deal with lists of vectors. We can construct any B-Spline \(\mathbf{x}(t)\) in \(S\) as
$$ \mathbf{x}(t) = \sum_{i=1}^m\mathbf{x}_ib_i(t) $$where the \(b_i\)'s are the scalar functions introduced in the section on Practical Construction. We refer to the vectors \(\mathbf{x}_i\in V\) as control points. Working with control points is an easy way to navigate and maneuver within \(S\). We like to think of these vectors as the controls of our very own (vector) space ship.
B-Spline Continuity
Our B-Spline construction process terminates with a collection of segments comprising a single piecewise polynomial function (the B-Spline) mapping from \(\mathbb{R}\) to a vector space such as \(\mathbb{R}^2\). This function is \(C^{\infty}\) on each of its segments, however not at the junctions between segments (i.e., at the knot points). We show here that the function can be \(C^{n-2}\) at these junctions, where \(n\) is the order of the function segments.
The following diagram shows the domains of segments from our B-Spline construction process that include the same knot. The knot is at the origin, and each domain is given by \([-a_i,b_i]\), with \(a_i,b_i\geq0\).
Each interval is the domain of a vector valued function \(\mathbf{f}_{ij}(t)\). To simplify our notation, we use \(\mathbf{f}_{ij}\) to denote \(\mathbf{f}_{ij}(0)\), that is, the value of the function at the knot value of interest. Doing this causes our interpolation scheme to become $$ \mathbf{f}_{k+1,l} = \mathbf{f}_{kl}\frac{b_l}{a_{k+l}+b_l}-\mathbf{f}_{k,l+1}\frac{a_{k+l}}{a_{k+l}+b_l} \tag{1} $$ These successive multiplications and additions can be considered all at once in a kind of elaborate Pascal's triangle, as follows.
Because \(a_7\) and \(b_1\) are zero, it's clear that the output value ? is zero. This means that \(\mathbf{f}_{61}=\mathbf{f}_{62}\), that is, that the B-Spline is at least \(C^0\) at the knot value in question. The right most hexagon which we isolate is an example of a B-Hive Sum. Establishing higher levels of continuity requires us to work with larger and larger B-Hive sums. Differentiating \((1)\), we obtain $$ \begin{align} \dot{\mathbf{f}}_{k+1,l} &= \dot{\mathbf{f}}_{kl}\frac{b_l}{a_{k+l}+b_l}-\dot{\mathbf{f}}_{k,l+1}\frac{a_{k+l}}{a_{k+l}+b_l} \nonumber \\ &+\frac{1}{a_{k+l}+b_l}(\mathbf{f}_{k,l+1}-\mathbf{f}_{kl}).\nonumber \end{align} $$ We use this to decompose \(\mathbf{f}^{[1]}_{61}\) as the sum \(\mathbf{f}^{[1]}_{61}.\!A+\mathbf{f}^{[1]}_{61}.\!B\). The \(A\) term consists of functions that have been differentiated one more time than those in the \(B\) term. We decompose \(\mathbf{f}^{[1]}_{62}\) similarly, and collect like terms (\(A\) and \(B\)) into their own B-Hive sums.
All B-Hive sums map to zero (see B-Hive Sum), and so it must be that \(\mathbf{A}\) and \(\mathbf{B}\) are both zero. From this it follows that \(\mathbf{f}^{[1]}_{61}=\mathbf{f}^{[1]}_{62}\), and so the B-Spline is \(C^1\) at the given knot.
We perform a similar decomposition for the second derivatives- the associated B-Hive sums are shown below, establishing that the B-Spline is \(C^2\) at the given knot point.
This pattern continues up to \(C^{n-2}\), with \(n=6\) in the current case. Generally, to establish \(C^k\), we have \(k+1\) B-Hive sums.
B-Hive Sums
A B-Hive sum is a linear map from \(\mathbb{R}^n\) to \(\mathbb{R}\) defined by \(n\) nonzero values \(a_i\) and \(n\) nonzero values \(b_i\) as shown below. The \(n\) input values \(x_i\) generate new values as they move from left to right along the tracks shown in gray. Encountering a value on a track means that you are multiplied by that value. Merging tracks correspond to a sum, and splitting tracks cause the same value to be sent in both of the new directions. The output is the value \(f(x_i)\) on the right.
The colored bands highlight a helpful pattern- the \(a_i\)'s are the same within each orange band, and the \(b_i\)'s are the same within each blue band.
The main result about B-Hive sums is that they map to zero. Regardless of the input values \(x_i\), the output is always zero. We establish this by induction.
base case
The property clearly holds when \(n=1\). If \(a_1\) and \(b_1\) are nonzero, then any \(x_1\) gets mapped to zero.
induction step
We suppose the property holds when \(n=3\), and we show that as a result, it must hold when \(n=4\). Start with four nonzero \(a_i\)'s, four nonzero \(b_i\)'s, and four numbers \(x_i\). Form a beehive sum with the first three \(a_i\)'s, the first three \(b_i\)'s, and the first three \(x_i\)'s.
We've labeled four points \(u_i\). Because the beehive sum is zero by hypothesis, it immediately follows that $$ u_0-\frac{1}{b_1}(u_1-\frac{1}{b_2}(u_2-\frac{1}{b_3}(b_3u_3)))=0. \tag{1} $$ Our next step is to form a B-Hive sum with all four of our \(a_i\)'s, \(b_i\)'s, and \(x_i\)'s. Note that the \(u_i\)'s in this new \(n=4\) B-Hive sum have the same values as in the \(n=3\) B-Hive sum.
The points marked \(v_i\) and \(w_i\) satisfy the following equations.
| \(w_0=\frac{1}{a_4}v_0-\frac{1}{b_1}w_1\), | ||
|---|---|---|
| \(v_0=u_0-\frac{1}{a_4+b_1}v_1\), | \(\;\;\;\;w_1=\frac{1}{a_4+b_1}v_1-\frac{1}{b_2}w_2\), | |
| \(v_1=u_1-\frac{1}{a_4+b_2}v_2\), | \(\;\;\;\;w_2=\frac{1}{a_4+b_2}v_2-\frac{1}{b_3}w_3\), | |
| \(v_2=u_2-\frac{1}{a_4+b_3}v_3\), | \(w_3=\frac{1}{a_4+b_3}v_3-x_4\), | |
| \(\!\!v_3=b_3u_3+a_4x_4\). |
Working through these equations from top to bottom, at each line substituting in the results from the next line and simplifying, we see that $$ w_0=\frac{1}{a_4}\left(u_0-\frac{1}{b_1}(u_1-\frac{1}{b_2}(u_2-\frac{1}{b_3}(b_3u_3)))\right).\notag $$ From \((1)\) and the fact that \(a_4\) is nonzero, it follows that \(w_0=0\), as desired. This same pattern gets us from any \(n\) to \(n+1\), and so it follows that all B-Hive sums map to zero.
References And Resources
Videos and interactive essays are a powerful way to share the kind of material we've covered in these notes. Some of my favorite examples of this are:
The videos made by 3b1b have fundamentally changed how many people think about math pedagogy. They are spectacular- beautifully executed, clear, and with an abundance of delightful ah-hah moments. They set the bar for explaining complex technical topics.
Bartosz Ciechanowski maintains a blog filled with beautiful and thought provoking interactive essays. Bartosz posts new essays all the time. One of my favorites is on Curves and Surfaces from 2021.
This is a free online book about all things Bezier. This book contains a thoughtful introduction to Bezier (and other) curves, derivations, and numerous interactive diagrams that have really helped me out.
Freya Holmer's videos are exquisitely produced and fun to watch. The two of her videos that most directly relate to curves are The Beauty of Bezier Curves, and The Continuity of Splines.
Distill was a scientific journal for machine learning research. These articles aren't related to curves, however they are examples of truly aspirational technical communication. Distill has been on an indefinite hiatus since 2021.
As much as I love videos and interactive websites, the foundation blocks of my own mathematical learning are still well written books. Some of my favorites that relate to B-Splines are as follows.
By Sheldon Axler
Linear Algebra Done Right is now an Open Access book that you can download here. This book is clear, sharp, and rigorous yet accessible. It was from this book that I really got some of the important ideas and theorems from linear algebra. I cannot recommend this book highly enough.
By Nick Trefethan and David Bao
A gentle introduction to applied linear algebra. I particularly like the development of the singular value decomposition. This book can be purchased here.
By Les A. Piegl and Wayne Tiller
Finding The NURBS Book was such a relief. Its development of NURBS is straightforward, clean, and easy to follow.
Additional resources:
By Patch KesslerThese notes of mine from a few years ago revolve around the representation of a perfect circle with NURBS. I play around with a simplex based extension of NURBS to higher dimensions, and prove that these higher dimensional NURBS can indeed perfectly represent higher dimensional circles.