Smooth Bézier Spline Through Prescribed Points
In this article we show you how to generate a smooth spline through a set of prescribed points. The article also contains an interactive SVG demo that allows you to try this in real time. This demo may not work if you are reading this post through email or on a mobile device. In that case, please visit the article from your computer browser.
Introduction
In scientific computing, we often need to construct a line that passes through a set of prescribed controlled points, or knots. A common example is mesh generation. Even if we utilize a mesh-free method, we need to somehow prescribe the outer extent for the problem. This is typically accomplished by specifying boundary splines (the counterpart in 3D are patches).
Bézier Splines
Spline is a collection of polygonal segments. The segments can be linear, quadratic, cubic, or even higher order polynomials. In this article we derive the equations needed to draw a smooth curve through a set of control points using the cubic Bézier polynomial. Wikipedia has a very nice article on Bézier curves that includes animations that show how these polynomials work.
The cubic Bézier curve is given by
$$ \mathbf{B}(t)=(1-t)^3\mathbf{P}_0 + 3(1-t)^2t\mathbf{P}_1 + 3(1-t)t^2\mathbf{P}_2 + t^3\mathbf{P}_3, \quad t\in [0,1]$$
this can be rewritten as
$$ \mathbf{B}(t)=(1-t)^3\mathbf{P}_0 + 3(t-2t^2+t^3)\mathbf{P}_1 + 3(t^2-t^3)\mathbf{P}_2 + t^3\mathbf{P}_3$$
In this definition, the points 0 and 3 correspond to the end points (the knots). The other two points are control points that determine the shape of the curve. The curve does not in general pass through these points. Following the steps outlined on codeproject.com, we want the first and second derivatives to be continuous across the spline boundary. This will give us two equations (one for each derivative) at each spline interface which we use to find the control points.
The first derivative is given by
$$ \mathbf{B}'(t) = -3(1-t)^2 \mathbf{P}_0 + 3(1-4t+3t^2)\mathbf{P}_1 + 3(2t-3t^2)\mathbf{P}_2 +3t^2\mathbf{P}_3$$
At the left boundary of an “i-th” segment we can write
$$ \mathbf{B}’_i(0)=\mathbf{B}’_{i-1}(1)$$
or
$$ -3\mathbf{P}_{0,i}+3\mathbf{P}_{1,i}=-3\mathbf{P}_{2,i-1}+3\mathbf{P}_{3,i-1}$$
Now, since the curve is continuous, \(\mathbf{P}_{0,i}=\mathbf{P}_{3,i-1}=\mathbf{K}_i\), the “i-th” knot point, and we can simplify this expression as (1)
$$ 2\mathbf{K}_i=\mathbf{P}_{1,i}+\mathbf{P}_{2,i-1} \quad\quad\quad(1)$$
We also want the second derivative to be continuous. The second derivative is given by
$$ \mathbf{B}”(t)= 6(1-t)\mathbf{P}_0+3(-4+6t)\mathbf{P}_1+3(2-6t)\mathbf{P}_2+6t\mathbf{P}_3$$
at the boundary we have
$$ 6 \mathbf{P}_{0,i}-12\mathbf{P}_{1,i}+6\mathbf{P}_{2,i}=6\mathbf{P}_{1,i-1}-12\mathbf{P}_{2,i-1}+6\mathbf{P}_{3,i-1}$$
Simplifying and taking into account the shared knot point, we get equation (2),
$$ -2\mathbf{P}_{1,i}+\mathbf{P}_{2,i}=\mathbf{P}_{1,i-1}-2\mathbf{P}_{2,i-1} \quad\quad\quad(2)$$
The equations (1) and (2) are defined only at the internal knots – places where two segments come together. Mathematically, this means that we have \(2(n-1)\) equations for \(2n\) unknowns. In order to close the system, we prescribe two more natural boundary conditions, \(\mathbf{B}”_0(0)=0\) and \(\mathbf{B}”_{n-1}(1)=0\). In other words, the spline becomes linear at the end points. These two remaining equations (3) and (4) are
$$ \mathbf{K}_0-2\mathbf{P}_{1,0}+\mathbf{P}_{2,0}=0 \quad\quad\quad(3)$$
and
$$ \mathbf{P}_{1,n-1}-2\mathbf{P}_{2,n-1}+\mathbf{K}_n=0 \quad\quad\quad(4)$$
We can further simplify this system by substituting (1) into (2) to obtain
$$ \mathbf{P}_{1,i-1}+4\mathbf{P}_{1,i}+\mathbf{P}_{1,i+1}=4\mathbf{K}_i+2\mathbf{K}_{i+1} \quad\quad i\in[1,n-2]$$.
On the boundary nodes we have from (3) and (4)
$$ 2\mathbf{P}_{1,0}+\mathbf{P}_{1,1}=\mathbf{K}_0+2\mathbf{K}_1$$
and
$$ 2\mathbf{P}_{1,n-2}+7\mathbf{P}_{1,n-1} = 8\mathbf{K}_{n-1}+\mathbf{K}_n$$
This is a tri-diagonal system for \(P_1\) which we can solve using the Thomas algorithm. Once we determine \(\mathbf{P}_1\), we obtain \( \mathbf{P}_2\) from equations (1) and (4),
$$ \mathbf{P}_{2,i}=2\mathbf{K}_{i}-\mathbf{P}_{1,i} \quad\quad i\in [0,n-2]$$
and
$$ \mathbf{P}_{2,n-1}= (1/2)(\mathbf{K}_n+\mathbf{P}_{1,n-1})$$
Interactive Demo
Below you will find an interactive demo that implements this algorithm. If everything loaded fine, you should see four yellow circles connected by a smooth line composed of three differently colored segments. You can move any of the circles with the mouse. If you don’t see the circles, cannot move them, or the spline is not smooth, first try reloading the page. If that still fails, try a different browser. The demo worked for me with Chrome, Firefox, and Internet Explorer 9 under Windows. I couldn’t get it working on my Android phone (no graphics rendered) and a friend reported it didn’t work on iPad either. If you know how to get the demo working on these mobile devices, please let me know by leaving a comment below.
This demo uses SVG and Javascript (You can find a TON of SVG examples on Prof. Dailey’s website. Another good resource is the SVG tutorial on “Peter’s Website”). SVG (Scalable Vector Graphics) is one of two ways you can get interactive graphics in your website (the other being HTML5’s <canvas> element). Of these two, SVG seems to be much more versatile. There are even free vector graphics programs out there, such as Inkscape, which save the images in the SVG format. The three lines connecting the circles are defined using SVG’s Bézier cubic path syntax. Whenever you drag a circle, the code recomputes the control points and updates the path definition. Your browser then automatically repaints the view. Pretty nifty! You will find the complete code below.
Source Code
You can download the source code here: circles.svg. The image is embedded using the <embed> tag. You can also see just the Javascript: bezier-spline.js (this code is embedded inside the SVG).
As a side note, it took me way too long to get the code running. Javascript is different from Java in that, among other things, it does not define variable types. Initially, the code was using
for (i=0;i<4;i++) { x[i]=V[i].getAttributeNS(null,"cx") y[i]=V[i].getAttributeNS(null,"cy") }
to obtain the coordinates of the four control points. These were then passed to another function which performed the actual computation. The problem was, this actually resulted in x[i] being a string and not an integer. Javascript did not complain about this. Instead, where required it automatically converted the string to an int. However, the plus operator resulted in string concatenation instead of addition of two numbers. For example, for x[i]="50", x[i]+30 resulted in 5030 instead of the expected value of 80. For the longest time I was getting really funky results and could not figure out why despite rechecking my math multiple times. I finally found the bug by stepping through the code with the Javascript debugger that comes with Chrome. The simple fix was to force the conversion to integer,
for (i=0;i<4;i++) { /*use parseInt to convert string to int*/ x[i]=parseInt(V[i].getAttributeNS(null,"cx")) y[i]=parseInt(V[i].getAttributeNS(null,"cy")) }
Also make sure to checkout the interactive online mesh generator follow up article and an article on computing intersection between a cubic Bezier curve and a line.
You will also find a lot of information on Bezier splines, including interactive HTML5 canvas demos at http://processingjs.nihongoresources.com/bezierinfo/
I also just found out this doesn’t work in Firefox 10 due to a bug that got fixed in FF11. In 10, the spline doesn’t update since the code to svg.getElementById(‘id’) returns an exception, “Component returned failure code: 0x80004001 (NS_ERROR_NOT_IMPLEMENTED) [nsIDOMSVGSVGElement.getElementById]”
Update 6/20/2012: I modified the code to create the SVG elements dynamically. This also eliminated the need for the getElementById call so the new version should work in FF10 as well.
The Bezier curve renders on my Asus Prime tablet with Android ICS (it will probably work on my iPad3 as well), and I did get some touch-related functionality working (but only partially).
Btw, you might enjoy some of this SVG eye-candy:
https://code.google.com/p/svg-filters-graphics/
It would be nice to see and example of a closed path
I made an example of a closed path:
http://www.jacos.nl/jacos_html/spline/circular/index.html
Excellent! I never imagined that this one article will generate so much interest.
Hi
Great, very useful.
But is there a typo in calculating P2 after solving the TDMA?
Should it be P2(i) = 2K(i+1) – P1(i+1)
This thread is pretty old and it is unlikely the author will fix this typo, but i’m pretty happy other people had discovered the same issue i found when i was trying to repeat the same steps in this demonstration.
I have found the same expression, to compute P2(i)
P2(i) = 2K(i+1) – P1(i+1)
This way, to obtain P2(n-1) it is impossible with this equation because P1(n) doesn’t exist, so, we can use equation (4) to compute it:
P2(n-1) = 0.5*(K(n)+P1(n-1))
Thanks, but I don’t find this useful at all. I am not a math whiz and it has been ages I’ve been in high school. So tracking through half a page of math is a chore in itself. The worst part, you don’t provide a simple, clear explanation that the average person can understand.
Look at this example: http://www.benknowscode.com/2012/09/path-interpolation-using-cubic-bezier_9742.html
That is clear and provides an easy to understand explanation of what – generally – I have to do to calculate the control points for quadratic curves. The explanation is so well written that I can come up with the math to do it without relying on actual code sample.
I’m looking for a similar, clear explanation to calculate control points for cubic beziers.
Hi Ted, I am glad you were able to find that other page more useful. In my case, I prefer just to see the equations, since in the end, that is what needs to be coded up.
The math here is not difficult at all. Bezier cubic is a (duh!) a cubic polynomial, evaluated from t=0 to t=1 between the left and right end point. Two other “knot” points control the shape of it in between. The whole point of finding the smooth spline is satisfying two requirements:
The first requirement basically says that the end point of spline 1 must equal to starting point of spline 2. The second condition tells us that the derivatives need to be equal. The more derivatives that equal, the smoother the transition. Here we take first and second derivative, which is the maximum you can take for a cubic before getting a constant. These two sets of derivatives are made equal to each other on spline 1 and 2. This gives you two equations for two unknowns, solve these, and you are done.
I didn’t read the other article in detail, but my guess is that the author is doing something similar, except comparing only the first derivative, since he is dealing with quadratics.
Thank you for this.
When you say “we have 2(n-1) equations with 2n unknowns”, you mean segments. Then later you say i belongs to an interval of [1, n-2], where n means number of original points.
Also, what’s the meaning of n for the boundary conditions?
The switch from P3,n-1 to Kn is confusing, especially since P0,0 is replaced with K0.
Could you please clarify whan n means in each case?
Oh, never mind, I got it. It’s all consistent, I just got confused in indexes.
I’ve implemented this yesterday and the results are far from ideal in certain cases.
This image says it all: https://dl.dropboxusercontent.com/u/2776279/SVG/weird_smoothing.png
And this is the actual curve in SVG: https://dl.dropboxusercontent.com/u/2776279/SVG/weird_smoothing.svg
And here is the code that produces it: https://gist.github.com/yay/99055816806340d2f1b4421e14d0191e
Any ideas how to restrict control points so that lie inside the [P0x, P3x] interval? After all the equation is parametric…
Also, there must be some condition that prevents spikes inside the curve segment itself, even if transitions are perfect.
Thank you very much for this explanation. I have adapted the code to provide a smooth-line functionality to a graph in SVG, which works for the most part. However in particular circumstances, where the consecutive points fluctuate a lot in the Y direction, the resulting curve can often go past the points, thus indicating a minimum/maximum that may be false.
I was wondering if you had any suggestions on how to use cubic bezier functions to smooth a line graph, without this issue? My initial guess would be to modify the linear equations for the points where the previous and next points are both less than (for a maximum) or greater than (for a minimum) that point in the Y direction, so that the first derivative would equal 0, although I’m not entirely sure how this would affect the computation for the x direction.
Hi Brad, one option may be to fit a standard cubic instead of the Bezier. That’s in fact what I ended up doing in one my own codes for data fitting. There is a good write up on fitting standard cubics to points in Burden’ Numerical Analaysis (and probably any other book like that).
I’m not sure if I’m understanding you correctly, but I’m limited to what kind of curves SVG is able to draw (quadratic bezier, cubic bezier and elliptic arcs) for the particular application I’m working on. I’m not sure if it can draw a standard cubic, unless I can adapt the cubic bezier curve to do it, or is that what you are suggesting?
Using the equations you have specified, I’m looking to replace the equations at points where either (K(i) > K(i-1) and K(i) > K(i+1)) or (K(i) < K(i-1) and K(i) < K(i+1)) so that the first derivative at those points is equal to zero (signifying that the graph is at a minimum or maximum). This is for the y-axis, the points are always increasing on the x-axis, so no points would trigger this condition for that axis
This means that K(i) = P(1,i) and K(i) = P(2,i-1). I'm not sure if what to do with the 2nd derivative equations, or whether I need to do anything at all with them. I think I should be able to use the above relationship directly into the matrix that is used to solve the system of equations, correct?
I tried the K(i) = P(1,i) and K(i) = P(2,i-1) in my code, and it produces close to what I would like, in that the graph doesn’t go past the points that are the local minimum and local maximum. However, at some of the points in the graphs I’ve tested, the curvature is not smooth, so I think I need to take the 2nd derivative into account. With substituting the above K(i) into your equation 2, I end up with P(2,i) = P(1,i-1), which I’m not sure how to put into the code. I think I may have done something wrong here.
Very useful – thanks
{The formulas get garbled after a small delay after loading into the browser (Chrome). Reloading the page restores the formulas, which then get garbled (Squished to the center and overlapping) again}
I was wondering if you could possibly explain which elements of the arrays are used for the x and y coordinates of the control points.
Possibly if you could modify the code so that it shows the final control points, and the “intermediate” control points of the left and right sections, it would be more comprehensible to me what is going on.
I’m trying to adapt your js code to “C++” (POV-Ray SDL) so I can better understand the algorithm and hopefully generate Bezier patches with 16 points and have my data be on the surface of the patch.
When I run the adapted code in POV-Ray, with the initial coordinates, I get
p1[0] = 100.0, 122.5, 0.0 p1[1] = 280.0, 350.1, 0.0 p1[2] = 500.0, 277.1, 0.0
p2[0] = 160.0, 249.9, 0.0 p2[1] = 340.0, 322.9, 0.0 p2[2] = 500.0, 277.1, 0.0
but I need to put this information into an array to return from my macro.
What does “return {p1:p1, p2:p2};” mean?
Why are p1[2] and p2[2] the same?
Thank you for any assistance you might offer. This has been the most helpful example I’ve found after days of searching how to determine the control points from the data! Great job!
Bill,
return {s1:p1, s2:p2}; returns an array that has 2 indices (“s1” and “s2”) that point to variables p1 and p2.
return {p1:p1, p2:p2}; means that computeControlPoints(..) returns array that has indices “p1” and “p2”.
These indices point to the arrays p1 and p2 that are defined in computeControlPoints.
This might be confusing indeed.
>>> I was wondering if you could possibly explain which elements of the arrays are used for the x and y coordinates of the control points.
>>> Why are p1[2] and p2[2] the same?
I suppose that these questions are related?
Essentially, the calculations are performed independently for all x and all y coordinates.
In updateSplines(), px.p1[i] then contains the x-coordinate for the first control point of a Bézier curve, px.p2[i] of the second control points.
The start and end of the curve are x[i] and x[i+1].
Hi Bill, thanks for pointing out the issue with the equation rendering. I confirmed this on my end but not sure what is the solution. I am using library called MathJax to show the equations. You can right click on the equations to get MathJax settings. If you go Math Settings -> Math Renderer and set it to Common HTML, then the equations should look fine. It seems that the default renderer is causing trouble in Chrome. I’ll investigate how to change this.
About your other question, the return statement returns an associative object which does not have an equivalent in C++. I am return “p1” as a named entity also confusingly called p1. Meaning, the
px = computeControlPoints(x);
call results in there being px.p1 and px.p2 arrays. If I am not mistaken, these arrays then define the x coordinate of the p1 point for each spline. So px.p1[0] is the x coordinate of the p1 point in the first segment, px.p1[1] is the x coordinate of p1 in the second segment, and so on.
Dr. Brieda,
Thanks for your reply, and I wish you luck with getting MathJax to behave.
I’ve spent the last several days looking over the javascript code, and rewrote it to function in POV-Ray to do x, y, and z coordinates. I better understand the output now, although the meaning of the algorithm’s stepwise functions remains unclear.
(I was an organic chemist – not an engineer or a mathematician, so my linear algebra is self taught, and toddler-level 😉 )
The code I adapted “works”, in the sense that it will return control points that give me a smooth Bezier spline, however I’ve noticed some odd behaviour when my data points have negative numbers, or if the components of the data points being interpolated are arranged backwards (in descending order) or are close together and linear. I don’t notice this when I manipulate the data in the SVG.
The different results are curious.
http://news.povray.org/web.56ded7545a67670a5e7df57c0%40news.povray.org
There are values that are hard-coded for a, b, c, and r – I’ve seen some codings of this algorithm that don’t declare these static values for the matrix – are these coefficients derived from some sort of parametric equation for the whole Bezier curve?
Thank you again for putting this together – I’ve learned more from this over the last few days than I have all year.
thank you so much for this example and explanation. extremely helpful!
This page helped me enormously!
(both for the maths and the programming!)
When I checked the algorithm, I found a small error that becomes apparent when the spline has only 2 nodes.
The problem is with
2 P[1,n-2] + 7 P[1,n-1] = 8 K[n-1] + K [n]
I corrected it as below, and implemented it on http://www.jacos.nl/jacos_html/spline/
(note that I use a bigger value for n, so I calculate as many times P1 as there are nodes, although P1[n] is only used to calculate P2[n-1] and its coordinates don’t figure in the SVG).
/*computes control points given knots K, this is the brain of the operation*/
function computeControlPoints(K)
{
var p1, p2, n
var a,b,c,r
p2=new Array();
n = K.length; //different from Lubos Brieda’s version
/*rhs vector*/
a=new Array();
b=new Array();
c=new Array();
r=new Array();
/*left most segment*/
a[0]=0; // outside the matrix
b[0]=2;
c[0]=1;
r[0] = K[0]+2*K[1];
/*internal segments*/
for (i = 1; i < n – 1; i++)
{
a[i]=1;
b[i]=4;
c[i]=1;
r[i] = 4 * K[i] + 2 * K[i+1];
}
/*right segment*/
a[n-1]=1;//different from Lubos Brieda's version
b[n-1]=2;
c[n-1]=0; // outside the matrix
r[n-1] = 3*K[n-1]; //different from Lubos Brieda's version
/*solves Ax=b with the Thomas algorithm (from Wikipedia)*/
for (i = 1; i = 0; –i)
p1[i] = (r[i] – c[i] * p1[i+1]) / b[i];
/*we have p1, now compute p2*/
for (i=0;i<n-1;i++)
p2[i]=2*K[i+1]-p1[i+1];
/* the last element of p1 is only used to calculate p2 */
p1.splice(n-1,1) // remove the last element
return {p1:p1, p2:p2};
}
My current version sets the control points nearer to the knots, if knots are near each other.
In this way, I avoid the “overshoot” when dots are near each other.
The first and second derivatives are still continuous.
Very cool article. Thanks for putting it together.
I believe there is a typo in the second to last equation before the section heading “Interactive Demo.” That equation should read
P_{2,i} = 2K_{i+1} – P_{1,i+1}
rather than 2K_i – P_{1,i}.
Thanks again!
It’d be nice to put a real license on the code, like BSD or X11. The ad hoc text in there unfortunately isn’t clear enough to use it (eg. it says I can use the “algorithm”, but it doesn’t actually say I can distribute or modify the source code, which is what a license needs to do)…
I recommend: https://opensource.org/licenses/BSD-2-Clause
Is there a way to kinda control the “radius” (I mean control the “smoothness”) at every point ?
Where to put the variable ?
Hi, I am not familiar with JS and trying to convert this code into Java. Where is the actual plotting of the curved line happening? I think it’s in this section:
for (i=0;i<3;i++)
S[i].setAttributeNS(null,"d", path(x[i],y[i],px.p1[i],py.p1[i],px.p2[i],py.p2[i],x[i+1],y[i+1]));
Could you spell out exactly how I can plot the line from the initial points, px, and py?
Nevermind, figured it out
Hi, the article is awesome and concise.
I’m not good at math stuff, the domain of some equations confuses me a lot, like getting a successive equation by adding one to i, while the end of domain decreased by 1, the start is untouched.
This can lead to a even bigger confusion to me: what’s the purpose and meaning of P_{1,0}? The zero segment confuses me a lot.
Would you please expand the explanation a little more?thank you.