使用贝塞尔曲线绘制闭合平滑曲线
4.31/5 (13投票s)
如何使用贝塞尔绘图图元在二维点集上绘制平滑的闭合曲线。
引言
2008年12月18日，我发布了一篇关于如何使用贝塞尔绘图原语绘制通过一组 2D 点的平滑曲线的 文章。 昨天我收到了一个关于如何以同样方式绘制闭合曲线的问题。 确实，它与开放式曲线有所不同。
方程
注意：为了使一系列单独的贝塞尔曲线成为样条线，我们应该计算贝塞尔控制点，以便样条曲线在结点处具有两个连续的导数。
单区间上的贝塞尔曲线表示为
B(t)=(1-t)3P0+3(1-t)2tP1+3(1-t)t2P2+t3P3
其中 t 位于 [0,1] 且
P0– 第一个结点P1– 第一个控制点（在第一个结点）P2– 第二个控制点（在第二个结点）P3– 第二个结点
一阶导数为
B’(t) = -3(1-t)2P0+3(3t2–4t+1)P1 +3(2t–3t2)P2+3t2P3
二阶导数为
B’’(t) = 6(1-t)P0+3(6t-4)P1+3(2–6t)P2+6tP3
让我们考虑在具有 n (n > 2) 个点 Pi (0,..,n-1) 和 n 个子区间 Si (0,..,n-1) 的区间上的分段贝塞尔曲线。
Si = (Pi, Pi+1); (i=0,..,n-2)
Sn-1 = (Pn-1, P0)
在 Si (i=0,..,n-2) 处的贝塞尔曲线将是
Bi(t) = (1-t)3Pi+3(1-t)2tP1i+3(1-t)t2P2i+1+t3Pi+1; (i=0,..,n-2)
以及在 Sn-1 处的闭合贝塞尔曲线
Bi(t) = (1-t)3Pn-1+3(1-t)2tP1n-1+3(1-t)t2P20+t3P0
在 Si (i=0,..,n-2) 处的一阶导数为
B’i(t) = -3(1-t)2Pi+3(3t2-4t+1)P1i+3(2t-3t2)P2i+1+3t2 Pi+1
以及在 Sn-1 处
B’n-1(t) = -3(1-t)2Pn-1+3(3t2-4t+1)P1n-1+3(2t-3t2)P20+3t2P0
一阶导数连续性条件给出
P1i+P2i = 2Pi; (i=0,..,n-1) (1)
在 Si (i=0,..,n-2) 处的二阶导数为
B’’i(t) = 6(1-t)Pi+6(3t-2)P1i+6(1-3t)P2i+1+6tPi+1
以及在 Sn-1 处
B’’n-1(t) = 6(1-t)Pn-1+6(3t-2)P1n-1+6(1-3t)P20+6tP0
二阶导数连续性条件给出
在 P0 处
2P10+P1n-1 = 2P20+P21 (2.1)
在 Pi (i=1,..,n-2) 处
2P1i+P1i-1 = 2P2i+P2i+1 (2.2)
以及在 Pn-1 处
2P1n-1+P1n-2 = 2P2n-1+P20 (2.3)
从 (2.1-3) 中排除 P2 和 (1)，我们得到 P1 的方程组
4P10+P11+P1n-1 = 4P0+2P1 for P0
P1i-1+4P1i+P1i+1 = 4 Pi+2Pi+1 for Pi (i=1,..,n-2)
P10+P1n-2+4P1n-1 = 4Pn-1+ 2P0 for Pn-1
我们得到了一个系统，其矩阵看起来像
4 1 0 0 ... 0 1
1 4 1 0 ... 0 0
0 1 4 1 ... 0 0
...............
1 0 0 0 ... 1 4
它是一个所谓的 "循环" 系统，可以使用来自书籍 "C 中的数值食谱" 的技巧来有效地解决，第 2.7 章 "稀疏线性系统"，主题 "循环三对角系统"。
在找到 P1 之后，很容易从 (1) 获得 P2。
代码
所有代码都附加在本文章的顶部。 这是一个 Visual Studio 2008 解决方案，针对 .NET 3.5，包含计算贝塞尔样条控制点的代码和演示其作用的示例应用程序。
目标是通过一组点（至少需要三个点）绘制闭合平滑曲线。 在我们有了这些点之后，我们需要使用 ClosedBezierSpline.GetCurveControlPoints 方法计算贝塞尔线段的控制点，连接相邻的点
/// <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]);
}
}
}
此方法非常简单：它获取结点，为 X 和 Y 点坐标组合上述循环系统，然后使用 Cyclic.Solve 方法解决它。 Cyclic.Solve 方法的代码是从 C 移植到 C# 的（参见书籍 "C 中的数值食谱"，第 2.7 章 "稀疏线性系统"，主题 "循环三对角系统"）。
/// <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;
}
}
此方法将循环系统转换为三对角系统，调用 Tridiagonal.Solve 来解决三对角系统，并将结果的逆变换应用于原始循环系统的解。 Tridiagonal.Solve 方法的代码是从 C 移植到 C# 的（参见书籍 "C 中的数值食谱"，第 2.4 章 “三对角和带状对角方程组”）。
/// <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;
}
}
此方法返回三对角系统的解。 它永远不会失败，因为我们的系统具有对角占优性。
该示例（WPF Windows 应用程序）演示了使用上述贝塞尔样条绘制采样圆。 您可以调整圆的半径、结点数量，并且可以显示/隐藏控制点，以查看它们与结点和贝塞尔线段的关系。
我用 C# 3.0 编译了这段代码，但我确信如果从类声明中删除关键字“static”，它可以在 C# 2.0 甚至 C# 1.0 中使用，而无需任何修改。