## Introduction

On December 18, 2008, I posted an article on how to draw a smooth curve through the set of 2D points with Bezier drawing primitives. Yesterday I got a question on how to draw a closed curve in the same manner. Indeed, there is a difference to open-ended curve.

## Equations

Note: 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 the knot points.

Bezier curve at the single interval is expressed as:

B(t)=(1-t)<sup>3</sup>P<sub>0</sub>+3(1-t)<sup>2</sup>tP<sub>1</sub>+3(1-t)t<sup>2</sup>P<sub>2</sub>+t<sup>3</sup>P<sub>3</sub>

where `t`

is in [0,1] and

`P`_{0}

– first knot point
`P`_{1}

– first control point (at first knot)
`P`_{2}

– second control point (at second knot)
`P`_{3}

– second knot point

First derivative is:

B’(t) = -3(1-t)<sup>2</sup>P<sub>0</sub>+3(3t<sup>2</sup>–4t+1)P<sub>1 </sub>+3(2t–3t<sup>2</sup>)P<sub>2</sub>+3t<sup>2</sup>P<sub>3</sub>

Second derivative is:

B’’(t) = 6(1-t)P<sub>0</sub>+3(6t-4)P<sub>1</sub>+3(2–6t)P<sub>2</sub>+6tP<sub>3</sub>

Let’s consider piece-wise Bezier curve on the interval with `n (n > 2)`

points `P`_{i} (0,..,n-1)

and n subintervals `S`_{i} (0,..,n-1)

.

S<sub>i</sub> = (P<sub>i</sub>, P<sub>i+1</sub>); (i=0,..,n-2)
S<sub>n-1</sub> = (P<sub>n-1</sub>, P<sub>0</sub>)

Bezier curve at `S`_{i} (i=0,..,n-2)

will be:

B<sub>i</sub>(t) = (1-t)<sup>3</sup>P<sub>i</sub>+3(1-t)<sup>2</sup>tP1<sub>i</sub>+3(1-t)t<sup>2</sup>P2<sub>i+1</sub>+t<sup>3</sup>P<sub>i+1</sub>; (i=0,..,n-2)

and the closing Bezier curve at `S`_{n-1}

:

B<sub>i</sub>(t) = (1-t)<sup>3</sup>P<sub>n-1</sub>+3(1-t)<sup>2</sup>tP1<sub>n-1</sub>+3(1-t)t<sup>2</sup>P2<sub>0</sub>+t<sup>3</sup>P<sub>0</sub>

First derivative at `S`_{i} (i=0,..,n-2)

will be:

B’<sub>i</sub>(t) = -3(1-t)<sup>2</sup>P<sub>i</sub>+3(3t<sup>2</sup>-4t+1)P1<sub>i</sub>+3(2t-3t<sup>2</sup>)P2<sub>i+1</sub>+3t<sup>2 </sup>P<sub>i+1</sub>

and at `S`_{n-1}

:

B’<sub>n-1</sub>(t) = -3(1-t)<sup>2</sup>P<sub>n-1</sub>+3(3t<sup>2</sup>-4t+1)P1<sub>n-1</sub>+3(2t-3t<sup>2</sup>)P2<sub>0</sub>+3t<sup>2</sup>P<sub>0</sub>

First derivative continuity condition gives:

P1<sub>i</sub>+P2<sub>i</sub> = 2P<sub>i</sub>; (i=0,..,n-1) (1)

Second derivative at `S`_{i} (i=0,..,n-2)

will be:

B’’<sub>i</sub>(t) = 6(1-t)P<sub>i</sub>+6(3t-2)P1<sub>i</sub>+6(1-3t)P2<sub>i+1</sub>+6tP<sub>i+1</sub>

and at `S`_{n-1}

:

B’’<sub>n-1</sub>(t) = 6(1-t)P<sub>n-1</sub>+6(3t-2)P1<sub>n-1</sub>+6(1-3t)P2<sub>0</sub>+6tP<sub>0</sub>

Second derivative continuity condition gives:

at `P`_{0}

:

2P1<sub>0</sub>+P1<sub>n-1</sub> = 2P2<sub>0</sub>+P2<sub>1</sub> (2.1)

at `P`_{i} (i=1,..,n-2)

:

2P1<sub>i</sub>+P1<sub>i-1</sub> = 2P2<sub>i</sub>+P2<sub>i+1</sub> (2.2)

and at `P`_{n-1}

:

2P1<sub>n-1</sub>+P1<sub>n-2</sub> = 2P2<sub>n-1</sub>+P2<sub>0</sub> (2.3)

Excluding `P2`

form (2.1-3) with (1), we get the set of equations for `P1`

:

4P1<sub>0</sub>+P1<sub>1</sub>+P1<sub>n-1</sub> = 4P<sub>0</sub>+2P<sub>1</sub> for P<sub>0</sub>
P1<sub>i-1</sub>+4P1<sub>i</sub>+P1<sub>i+1</sub> = 4 P<sub>i</sub>+2P<sub>i+1</sub> for P<sub>i</sub> (i=1,..,n-2)
P1<sub>0</sub>+P1<sub>n-2</sub>+4P1<sub>n-1</sub> = 4P<sub>n-1</sub>+ 2P<sub>0</sub> for P<sub>n-1</sub>

We got the system with the matrix which looks like:

4 1 0 0 ... 0 1
1 4 1 0 ... 0 0
0 1 4 1 ... 0 0
...............
1 0 0 0 ... 1 4

It's so-called "Cyclic" system which can be solved as effectively as a tridiagonal system with the trick from the book "Numerical Recipes in C", Chapter 2.7 "Sparse Linear Systems", topic "Cyclic Tridiagonal Systems".

After `P1`

's found, it's easy to get `P2`

's from (1).

## The Code

All the code is attached at the top of this article. It is a Visual Studio 2008 solution targeted to .NET 3.5 and contains both the code to calculate Bezier Spline control points and the sample application to see it in action.

The goal is to draw the closed smooth curve through the set of points (at least three points required). After we have the points, we need to calculate control points of the Bezier segments, joining adjacent points, with `ClosedBezierSpline.GetCurveControlPoints`

method:

public static class ClosedBezierSpline
{
public static void GetCurveControlPoints(Point[] knots,
out Point[] firstControlPoints, out Point[] secondControlPoints)
{
int n = knots.Length;
if (n <= 2)
{ firstControlPoints = new Point[0];
secondControlPoints = new Point[0];
return;
}
double[] a = new double[n], b = new double[n], c = new double[n];
for (int i = 0; i < n; ++i)
{
a[i] = 1;
b[i] = 4;
c[i] = 1;
}
double[] rhs = new double[n];
for (int i = 0; i < n; ++i)
{
int j = (i == n - 1) ? 0 : i + 1;
rhs[i] = 4 * knots[i].X + 2 * knots[j].X;
}
double[] x = Cyclic.Solve(a, b, c, 1, 1, rhs);
for (int i = 0; i < n; ++i)
{
int j = (i == n - 1) ? 0 : i + 1;
rhs[i] = 4 * knots[i].Y + 2 * knots[j].Y;
}
double[] y = Cyclic.Solve(a, b, c, 1, 1, rhs);
firstControlPoints = new Point[n];
secondControlPoints = new Point[n];
for (int i = 0; i < n; ++i)
{
firstControlPoints[i] = new Point(x[i], y[i]);
secondControlPoints[i] = new Point
(2 * knots[i].X - x[i], 2 * knots[i].Y - y[i]);
}
}
}

This method is very straightforward: it takes the knots, composes the cyclic system above for both X and Y point coordinates and then solves it with the `Cyclic.Solve`

method. `Cyclic.Solve`

method code is ported to C# from C (see the book "Numerical Recipes in C", Chapter 2.7 "Sparse Linear Systems", topic "Cyclic Tridiagonal Systems").

public static class Cyclic
{
public static double[] Solve(double[] a, double[] b,
double[] c, double alpha, double beta, double[] rhs)
{
if (a.Length != b.Length || c.Length != b.Length ||
rhs.Length != b.Length)
throw new ArgumentException
("Diagonal and rhs vectors must have the same size.");
int n = b.Length;
if (n <= 2)
throw new ArgumentException
("n too small in Cyclic; must be greater than 2.");
double gamma = -b[0]; double[] bb = new Double[n];
bb[0] = b[0] - gamma;
bb[n-1] = b[n - 1] - alpha * beta / gamma;
for (int i = 1; i < n - 1; ++i)
bb[i] = b[i];
double[] solution = Tridiagonal.Solve(a, bb, c, rhs);
double[] x = new Double[n];
for (int k = 0; k < n; ++k)
x[k] = solution[k];
double[] u = new Double[n];
u[0] = gamma;
u[n-1] = alpha;
for (int i = 1; i < n - 1; ++i)
u[i] = 0.0;
solution = Tridiagonal.Solve(a, bb, c, u);
double[] z = new Double[n];
for (int k = 0; k < n; ++k)
z[k] = solution[k];
double fact = (x[0] + beta * x[n - 1] / gamma)
/ (1.0 + z[0] + beta * z[n - 1] / gamma);
for (int i = 0; i < n; ++i)
x[i] -= fact * z[i];
return x;
}
}

This method converts Cyclic system to Tridiagonal system, calls `Tridiagonal.Solve`

to solve the Tridiagonal system and applies the reverse transformation of the result to the solution of the original Cyclic system. `Tridiagonal.Solve`

method code is ported to C# from C (see the book "Numerical Recipes in C", Chapter 2.4 “Tridiagonal and Band Diagonal Systems of Equations”).

public static class Tridiagonal
{
public static double[] Solve(double[] a, double[] b, double[] c, double[] rhs)
{
if (a.Length != b.Length || c.Length != b.Length ||
rhs.Length != b.Length)
throw new ArgumentException
("Diagonal and rhs vectors must have the same size.");
if (b[0] == 0.0)
throw new InvalidOperationException("Singular matrix.");
ulong n = Convert.ToUInt64(rhs.Length);
double[] u = new Double[n];
double[] gam = new Double[n];
double bet = b[0];
u[0] = rhs[0] / bet;
for (ulong j = 1;j < n;j++) {
gam[j] = c[j-1] / bet;
bet = b[j] - a[j] * gam[j];
if (bet == 0.0)
throw new InvalidOperationException
("Singular matrix.");
u[j] = (rhs[j] - a[j] * u[j - 1]) / bet;
}
for (ulong j = 1;j < n;j++)
u[n - j - 1] -= gam[n - j] * u[n - j];
return u;
}
}

This method returns the solution of the Tridiagonal System. It never fails as our system has diagonal dominance.

The example (WPF Windows Application) demonstrates drawing of the sampled circle with Bezier spline above. You can play with the circle radius, knot points count and you can show/hide control points to see how they relate to the knot points and Bezier segments.

I compiled this code in C# 3.0 but I'm sure it can be used without any modification in C# 2.0 and even in C# 1.0 if you remove keyword “`static`

” from the class declarations.