## How to draw a smooth curve through a set of 2D points with Bezier methods?

From time to time I face with the question: how to draw a smooth curve through a set of 2D points? It seems strange but we haven’t out of the box primitives to do so. Yes, we can draw a polyline, Bezier polyline or piece-wise cardinal spline but it all isn’t what is desired: polyline isn’t smooth; with Bezier and cardinal spline we have headache with additional parameters like Bezier control points or tension. Very often we haven’t any data to evaluate these additional parameters. Very often all we need is to draw the curve which passes through the points given and is smooth enough. Let’s try to find the solution.

We have interpolation methods at hands. The cubic spline variations, for example, give satisfactory results in most cases. We could use it and draw the result of interpolation, but there are some nasty drawbacks:

1. Cubic spline is a cubic polynomial but neither of Win32, .NET Forms or WPF provide methods to draw piecewise cubic polynomial curve. So, to draw that spline we have to invent an effective algorithm to approximate a cubic polynomial with polyline, which isn’t trivial (word “effective” is a point!).
2. Cubic spline algorithms are usually designed to deal with functions in the Cartesian coordinates, `y=f(x)`. This implies that the function is single valued; it isn’t always appropriate.

On other hand, the platform provides us with the suite of Bezier curve drawing methods. Bezier curve is a special representation of a cubic polynomial expressed in the parametric form (so it isn’t subject of single valued function restriction). Bezier curve starts and end with two points often named “knots”; form of that curve is controlled by two more points known as “control points”.

Bezier spline is a sequence of individual Bezier curves joined to form a whole curve. The trick to make it a spline is to calculate control points in such a way that the whole spline curve has two continuous derivatives.

I spent some time Googling for a code in any C-like language for Bezier spline, but haven’t find any cool, ready-to-use code.

# Equations

One more time: to make a sequence of individual Bezier curves to be a spline we should calculate Bezier control points so that the spline curve has two continuous derivatives at knot points.

Here we’ll deal with open-ended curves but the same approach could be applied to the closed curves.

Bezier curve on single interval is expressed as

B(t)=(1-t)3P0+3(1-t)2tP1+3(1-t)t2P2+t3P3

where t is in [0,1] and

1. P0 – first knot point
2. P1 – first control point (at first knot)
3. P2 – second control point (at second knot)
4. P3 – second knot point

First derivative is:

B’(t) = -3 (1 – t)2 P0 + 3 (3 t2 – 4t + 1) P1 + 3 (2t – 3t2) P2 + 3t2 P3

Second derivative is:

B’’(t) = 6 (1 – t) P0 + 3 (6t – 4) P1 + 3 (2 – 6t) P2+ 6 t P3

Let’s consider piece-wise Bezier curve with n+1 points (0,..,n) and, accordingly, n subintervals.

At the i-th subinterval Bezier curve will be:

Bi(t) = (1 – t)3 Pi-1 + 3 (1-t)2 t P1i + 3 (1-t) t2 P2i + t3 Pi; (i=1,..,n)

First derivative at the i-th subinterval:

B’i(t) = -3 (1-t)2 Pi-1 + 3(3t2 – 4t + 1) P1i + 3(2t – 3t2) P2i+3t2 Pi; (i=1,..,n)

First derivative continuity condition B’i-1(1)= B’i(0) give:

P1i + P2i-1 = 2Pi-1; (i=2,..,n) (1)

Second derivative at the i-th subinterval:

B’’i (t) = 6 (1-t) Pi-1 +6 (3t – 2) P1i + 6 (1 – 3t) P2i + 6t Pi; (i=1,..,n)

Second derivative continuity condition B’’i-1 (1)= B’’i (0) give:

P1i-1 + 2P1i = P2i + 2P2i-1; (i=2,..,n) (2)

Then, as always with spline’s, we’ll add two more conditions at the ends of total interval. Lets it’ll be the “natural conditions” B’’1 (0) = 0 and B’’n (1) = 0

2P11 – P21 = P0 (3)

2P2n – P1n = Pn (4)

Now we have 2n conditions (1-4) for n control points P1 and n control points P2. Excluding P2 we’ll have n equations for n control points P1:

2P11 + P12 = P0 + 2P1 (i=1)

….

P1i-1 + 4P1i + P1i+1 = 4Pi-1 + 2Pi (i=2,..,n-1) (5)

….

P1n-1 + 2P1n= 3Pn-1 (i=n)

System (5) is the tridiagonal one with diagonal dominance hence we can solve it by simple variable elimination without any tricks.

When P1 is found P2 could be calculated from (1) and (4).

# The code

```/// <summary>
/// Bezier Spline methods
/// </summary>
public static class BezierSpline
{
/// <summary>
/// Get open-ended Bezier Spline Control Points.
/// </summary>
/// <param name="knots">Input Knot Bezier spline points.</param>
/// <param name="firstControlPoints">Output First Control points array of knots.Length - 1 length.</param>
/// <param name="secondControlPoints">Output Second Control points array of knots.Length - 1 length.</param>
public static void GetCurveControlPoints(Point[] knots, out Point[] firstControlPoints, out Point[] secondControlPoints)
{
int n = knots.Length - 1;
if (n < 1)
{
firstControlPoints = new Point[0];
secondControlPoints = new Point[0];
return;
}

// Calculate first Bezier control points
// Right hand side vector
double[] rhs = new double[n];

// Set right hand side X values
for (int i = 1; i < n - 1; ++i)
rhs[i] = 4 * knots[i].X + 2 * knots[i + 1].X;
rhs[0] = knots[0].X + 2 * knots[1].X;
rhs[n - 1] = 3 * knots[n - 1].X;
// Get first control points X-values
double[] x = GetFirstControlPoints(rhs);

// Set right hand side Y values
for (int i = 1; i < n - 1; ++i)
rhs[i] = 4 * knots[i].Y + 2 * knots[i + 1].Y;
rhs[0] = knots[0].Y + 2 * knots[1].Y;
rhs[n - 1] = 3 * knots[n - 1].Y;
// Get first control points Y-values
double[] y = GetFirstControlPoints(rhs);

// Fill output arrays.
firstControlPoints = new Point[n];
secondControlPoints = new Point[n];
for (int i = 0; i < n; ++i)
{
// First control point
firstControlPoints[i] = new Point(x[i], y[i]);
// Second control point
if (i < n - 1)
secondControlPoints[i] = new Point(2 * knots[i + 1].X - x[i + 1], 2 * knots[i + 1].Y - y[i + 1]);
else
secondControlPoints[i] = new Point((knots[n].X + x[n - 1]) / 2, (knots[n].Y + y[n - 1]) / 2);
}
}

/// <summary>
/// Solves a tridiagonal system for one of coordinates (x or y) of first Bezier control points.
/// </summary>
/// <param name="rhs">Right hand side vector.</param>
/// <returns>Solution vector.</returns>
private static double[] GetFirstControlPoints(double[] rhs)
{
int n = rhs.Length;
double[] x = new double[n]; // Solution vector.
double[] tmp = new double[n]; // Temp workspace.

double b = 2.0;
x[0] = rhs[0] / b;
for (int i = 1; i < n; i++) // Decomposition and forward substitution.
{
tmp[i] = 1 / b;
b = (i < n - 1 ? 4.0 : 2.0) - tmp[i];
x[i] = (rhs[i] - x[i - 1]) / b;
}
for (int i = 1; i < n; i++)
x[n - i - 1] -= tmp[n - i] * x[n - i]; // Backsubstitution.
return x;
}
}
```

Although I compiled this code in C# 3.0 I don’t see why it can’t be used without any modification in C# 2.0 and even in C# 1.0 if you remove keyword “static” from the class declaration.

Note that the code uses System.Windows.Point structure so it is intended for use with Windows Presentation Foundation and using System.Windows directive is required. But all you should do to use it with Windows Forms is to replace System.Windows.Point type with System.Drawing.PointF. Equally I see straightforward to convert this code to C/C++ if desired.

# The sample

http://cid-39d56f0c7a08d703.skydrive.live.com/embedrowdetail.aspx/.Public/BezierSpline/BezierSpline.zip

It is Visual Studio 2008 solution targeted to .NET 3.5. It contains WPF Windows Application project designed to demonstrate some curves drawn with Bezier spline above. You can select one of the curves from Combo Box at the top of the Window, experiment with point counts and set appropriate XY Scales. You can even add you own curve, but this requires coding as follows: