Draw Closed Smooth Curve with Bezier Spline

• Download source - 17.53 KB

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.