14,025,186 members
alternative version

#### Stats

87.1K views
42 bookmarked
Posted 2 Mar 2009
Licenced CPOL

# Draw Closed Smooth Curve with Bezier Spline

How to draw a smooth closed curve through the set of 2D points with Bezier drawing primitives

## 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

1. `P<sub>0</sub>` – first knot point
2. `P<sub>1</sub>` – first control point (at first knot)
3. `P<sub>2</sub>` – second control point (at second knot)
4. `P<sub>3</sub>` – 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<sub>i</sub> (0,..,n-1)` and n subintervals `S<sub>i</sub> (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<sub>i</sub> (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<sub>n-1</sub>`:

`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<sub>i</sub> (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<sub>n-1</sub>`:

`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<sub>i</sub> (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<sub>n-1</sub>`:

`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<sub>0</sub>`:

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

at `P<sub>i</sub> (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<sub>n-1</sub>`:

`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:

```/// <summary>
/// Closed Bezier Spline Control Points calculation.
/// </summary>
public static class ClosedBezierSpline
{
/// <summary>
/// Get Closed Bezier Spline Control Points.
/// </summary>
/// <param name="knots">Input Knot Bezier spline points.</param>
/// <param name="firstControlPoints">
/// Output First Control points array of the same
/// length as the <paramref name="knots"> array.</param>
/// <param name="secondControlPoints">
/// Output Second Control points array of the same
/// length as the <paramref name="knots"> array.</param>
public static void GetCurveControlPoints(Point[] knots,
out Point[] firstControlPoints, out Point[] secondControlPoints)
{
int n = knots.Length;
if (n <= 2)
{ // There should be at least 3 knots to draw closed curve.
firstControlPoints = new Point[0];
secondControlPoints = new Point[0];
return;
}

// Calculate first Bezier control points

// The matrix.
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;
}

// Right hand side vector for points X coordinates.
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;
}
// Solve the system for X.
double[] x = Cyclic.Solve(a, b, c, 1, 1, rhs);

// Right hand side vector for points Y coordinates.
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;
}
// Solve the system for Y.
double[] y = Cyclic.Solve(a, b, c, 1, 1, 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.
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").

```/// <summary>
/// Solves the cyclic set of linear equations.
/// </summary>
/// <remarks>
/// The cyclic set of equations have the form
/// ---------------------------
/// b0 c0  0 · · · · · · ß
///	a1 b1 c1 · · · · · · ·
///  · · · · · · · · · · ·
///  · · · a[n-2] b[n-2] c[n-2]
/// a  · · · · 0  a[n-1] b[n-1]
/// ---------------------------
/// This is a tridiagonal system, except for the matrix elements
/// a and ß in the corners.
/// </remarks>
public static class Cyclic
{
/// <summary>
/// Solves the cyclic set of linear equations.
/// </summary>
/// <remarks>
/// All vectors have size of n although some elements are not used.
/// The input is not modified.
/// </remarks>
/// <param name="a">Lower diagonal vector of size n; a[0] not used.</param>
/// <param name="b">Main diagonal vector of size n.</param>
/// <param name="c">Upper diagonal vector of size n; c[n-1] not used.</param>
/// <param name="alpha">Bottom-left corner value.</param>
/// <param name="beta">Top-right corner value.</param>
/// <param name="rhs">Right hand side vector.</param>
/// <returns>The solution vector of size n.</returns>
public static double[] Solve(double[] a, double[] b,
double[] c, double alpha, double beta, double[] rhs)
{
// a, b, c and rhs vectors must have the same size.
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]; // Avoid subtraction error in forming bb[0].
// Set up the diagonal of the modified tridiagonal system.
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];
// Solve A · x = rhs.
double[] solution = Tridiagonal.Solve(a, bb, c, rhs);
double[] x = new Double[n];
for (int k = 0; k < n; ++k)
x[k] = solution[k];

// Set up the vector u.
double[] u = new Double[n];
u[0] = gamma;
u[n-1] = alpha;
for (int i = 1; i < n - 1; ++i)
u[i] = 0.0;
// Solve A · z = u.
solution = Tridiagonal.Solve(a, bb, c, u);
double[] z = new Double[n];
for (int k = 0; k < n; ++k)
z[k] = solution[k];

// Form v · x/(1 + v · z).
double fact = (x[0] + beta * x[n - 1] / gamma)
/ (1.0 + z[0] + beta * z[n - 1] / gamma);

// Now get the solution vector x.
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”).

```/// <summary>
/// Tridiagonal system solution.
/// </summary>
public static class Tridiagonal
{
/// <summary>
/// Solves a tridiagonal system.
/// </summary>
/// <remarks>
/// All vectors have size of n although some elements are not used.
/// </remarks>
/// <param name="a">Lower diagonal vector; a[0] not used.</param>
/// <param name="b">Main diagonal vector.</param>
/// <param name="c">Upper diagonal vector; c[n-1] not used.</param>
/// <param name="rhs">Right hand side vector</param>
/// <returns>system solution vector</returns>
public static double[] Solve(double[] a, double[] b, double[] c, double[] rhs)
{
// a, b, c and rhs vectors must have the same size.
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.");
// If this happens then you should rewrite your equations as a set of
// order N - 1, with u2 trivially eliminated.

ulong n = Convert.ToUInt64(rhs.Length);
double[] u = new Double[n];
double[] gam = new Double[n]; 	// One vector of workspace,
// gam is needed.

double bet = b[0];
u[0] = rhs[0] / bet;
for (ulong j = 1;j < n;j++) // Decomposition and forward substitution.
{
gam[j] = c[j-1] / bet;
bet = b[j] - a[j] * gam[j];
if (bet == 0.0)
// Algorithm fails.
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]; // Backsubstitution.

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.

## About the Author

 Team Leader Russian Federation
No Biography provided

 Pro

## Comments and Discussions

 First Prev Next
 Area of Closed Smooth Curve Matt Flinders13-Apr-11 5:28 Matt Flinders 13-Apr-11 5:28
 hi,there is something wrong with the source file grace_shuyan7-Dec-09 11:58 grace_shuyan 7-Dec-09 11:58
 Re: hi,there is something wrong with the source file Oleg V. Polikarpotchkin7-Dec-09 19:08 Oleg V. Polikarpotchkin 7-Dec-09 19:08
 Re: hi,there is something wrong with the source file grace_shuyan7-Dec-09 22:58 grace_shuyan 7-Dec-09 22:58
 Thanks for your reply.It is WinRAR unzip error. WinRAR:Diagnostic messages are as follows: ! C:\Documents and Settings\Junjun\Desktop\ClosedBezierSpline.zip: Cannot create ClosedBezierSpline - 颸ClosedBezierSpline.sln ! The filename, directory name, or volume label syntax is incorrect. ! C:\Documents and Settings\Junjun\Desktop\ClosedBezierSpline.zip: Cannot create folder ClosedBezierSpline - 颸DrawSampledCircle ! The filename, directory name, or volume label syntax is incorrect. I have tried to unzip it on different computers,but it all ends up with the same unzip error. Many thanks for your help. Regards Grace
 Re: hi,there is something wrong with the source file Oleg V. Polikarpotchkin8-Dec-09 20:00 Oleg V. Polikarpotchkin 8-Dec-09 20:00
 Re: hi,there is something wrong with the source file grace_shuyan8-Dec-09 23:49 grace_shuyan 8-Dec-09 23:49
 Re: hi,there is something wrong with the source file Oleg V. Polikarpotchkin9-Dec-09 18:36 Oleg V. Polikarpotchkin 9-Dec-09 18:36
 Re: hi,there is something wrong with the source file grace_shuyan15-Dec-09 3:50 grace_shuyan 15-Dec-09 3:50
 Tangled curves SteveAbbott24-Mar-09 6:10 SteveAbbott 24-Mar-09 6:10
 Re: Tangled curves Oleg V. Polikarpotchkin24-Mar-09 20:01 Oleg V. Polikarpotchkin 24-Mar-09 20:01
 Re: Tangled curves SteveAbbott25-Mar-09 7:25 SteveAbbott 25-Mar-09 7:25
 What is a "Tridiagonal System" ? ajstadlin2-Mar-09 17:24 ajstadlin 2-Mar-09 17:24
 Re: What is a "Tridiagonal System" ? Oleg V. Polikarpotchkin2-Mar-09 18:45 Oleg V. Polikarpotchkin 2-Mar-09 18:45
 Re: What is a "Tridiagonal System" ? ajstadlin3-Mar-09 3:13 ajstadlin 3-Mar-09 3:13
 formatting... Adrian Pasik2-Mar-09 10:08 Adrian Pasik 2-Mar-09 10:08
 Last Visit: 19-Apr-19 16:21     Last Update: 19-Apr-19 16:21 Refresh 1

General    News    Suggestion    Question    Bug    Answer    Joke    Praise    Rant    Admin

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.