65.9K
CodeProject 正在变化。 阅读更多。
Home

使用贝塞尔曲线绘制一组 2D 点的平滑曲线

starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.94/5 (55投票s)

2008 年 12 月 18 日

CPOL

6分钟阅读

viewsIcon

457651

downloadIcon

12895

计算分段贝塞尔曲线控制点以形成样条。

引言

我时不时会遇到这样一个问题:如何绘制通过一组二维点的平滑曲线?这似乎很奇怪,但我们没有现成的图元可以做到这一点。是的,我们可以绘制折线、贝塞尔折线或分段基样条曲线,但它们都不是我们想要的:折线不平滑;对于贝塞尔曲线和基样条曲线,我们需要处理额外的参数,如贝塞尔控制点或张力。很多时候,我们没有数据来评估这些额外的参数。而且,很多时候,我们需要的只是绘制通过给定点的足够平滑的曲线。让我们试着找到一个解决方案。

我们手头有插值方法。例如,三次样条曲线的变体在大多数情况下都能给出令人满意的结果。我们可以使用它并绘制插值结果,但有一些令人讨厌的缺点。

  1. 三次样条曲线是三次多项式,但 Win32、.NET Forms 或 WPF 不提供绘制分段三次多项式曲线的方法。因此,要绘制该样条曲线,我们必须发明一种有效的算法来用折线近似三次多项式,这并非易事(“有效”这个词是关键!)。
  2. 三次样条算法通常设计用于处理笛卡尔坐标系中的函数,y=f(x)。这意味着该函数是单值函数;这并非总是合适的。

另一方面,平台为我们提供了一套贝塞尔曲线绘制方法。贝塞尔曲线是三次多项式的一种特殊表示形式,以参数形式表示(因此它不受单值函数限制)。贝塞尔曲线以通常称为“节点”的两个点开始和结束;曲线的形状由另外两个点(称为“控制点”)控制。

贝塞尔样条曲线是由一系列独立的贝塞尔曲线连接而成的一条完整曲线。实现样条曲线的关键在于以这样一种方式计算控制点,使得整个样条曲线具有两个连续的导数。

我花了一些时间在 Google 上搜索任何类似 C 语言的贝塞尔样条曲线代码,但没有找到任何很酷、即用型的代码。

在这里,我们将处理开放式曲线,但相同的方法也可以应用于闭合曲线。

贝塞尔曲线表示

在单个区间上的贝塞尔曲线表示为

$B(t)=(1-t)^3 P_0 +3(1-t)^2t P_1 +3(1-t)t^2 P_2 +t^3 P_3 \ldots (1)$

其中 t[0,1] 中,并且

  1. P0– 第一个节点
  2. P1– 第一个控制点(靠近 P0
  3. P2– 第二个控制点(靠近 P3
  4. P3– 第二个节点

(1) 的一阶导数是

$B'(t)=-3(1-t)^2 P_0 +3(3t^2 - 4t+1) P_1 + 3(2t-3t^2)P_2 +3t^2 P_3$

(1) 的二阶导数是

$B''(t)=6(1 - t)P_0 +3(6t-4)P_1+3(2 - 6t)P_2 +6tP_3$

单段

如果我们只有两个节点,我们的“平滑”贝塞尔曲线应该是一条直线,即在 (1) 中,带有 2 和 3 次幂的项的系数应该为零。很容易推断出其控制点应计算为

$3P_1 = 2P_0 + P_3 P_2 = 2P_1 - P_0$

多段

这是我们有超过两个点的情况。再说一次:要使一系列独立的贝塞尔曲线成为样条曲线,我们应该计算贝塞尔控制点,使样条曲线在节点处具有两个连续的导数。

考虑一组具有 n+1 个点和 n 个子区间的 piecewise Bezier 曲线,第 (i-1) 条曲线应连接到第 i 条。现在我们将点表示如下:

  1. Pi– 第 ith 个节点 (i=1,..,n)
  2. P1i– 靠近第一个控制点Pi
  3. P2i– 靠近第二个控制点Pi

在第 ith 个子区间,贝塞尔曲线将是

$B_i (t)=(1-t)^3 P_{i-1} +3(1-t)^2 tP1_i +3(1-t)t^2 P2_i +t^3 P_i \ldots (i=1,..,n)$

在第 ith 个子区间的一阶导数是

$B'_i (t) = -3(1-t)^2 P_i-1 +3(3t^2 - 4t+1)P1_i +3(2t-3t^2)P2_i +3t^2 P_i; \ldots (i=1,..,n)$

一阶导数连续性条件 \(B'_{i-1}(1) = B'_i(0)\) 给出

$P1_i +P2_{i-1} = 2P{i-1}; \ldots (i=2,..,n) (2)$

在第 ith 个子区间的二阶导数是

$B''_i(t) = 6(1-t)P_{i-1} + 6(3t-2)P1_i + 6(1-3t)P2_i +6tP_i; \ldots (i=1,..,n)$

二阶导数连续性条件 B’’i-1(1)=B’’i(0) 给出

$P1_{i-1} + 2P1_i = P2_i +2P2_{i-1}; \ldots (i=2,..,n) (3)$
P1i-1+2P1i=P2i+2P2i-1; (i=2,..,n) (3)

然后,像往常一样使用样条曲线,我们将在总区间的两端添加另外两个条件。这些将是“自然条件” B’’<sub>1</sub>(0) = 0B’’<sub>n</sub>(1) = 0

$2P1_1 - P2_1 = P_0 (4) 2P2_n -P1_n = P_n \ldots (5)$
2P11-P21 = P0 (4) 2P2n-P1n = Pn (5)

现在,我们有 2n 个条件 (2-5) 用于 n 个控制点 P1n 个控制点 P2。排除 P2,我们将得到 n 个方程用于 n 个控制点 P1

$ \begin{align} 2P1_1 + P1_2 &= P_0 + 2P_1 P1_1 + 4P1_2 + P1_3 = 4P_1 + 2P_2 \\ P1_{i-1} + 4P1_i + P1_{i+1} &= 4P_{i-1} + 2P_i \ldots (6) \\ P1_{n-2} + 4P1_{n-1} + P1_n &= 4P_{n-2} + 2P_{n-1} 2P1_{n-1} + 7P1_n = 8P_{n-1} + P_n \end{align} $

系统 (6) 是对角占优的三对角系统,因此我们可以通过简单的变量消除来解决它,无需任何技巧。

当找到 P1 时,可以从 (2) 和 (5) 计算 P2

代码

/// <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>
	/// <exception cref="ArgumentNullException"><paramref name="knots"/>
	/// parameter must be not null.</exception>
	/// <exception cref="ArgumentException"><paramref name="knots"/>
	/// array must contain at least two points.</exception>
	public static void GetCurveControlPoints(Point[] knots,
		out Point[] firstControlPoints, out Point[] secondControlPoints)
	{
		if (knots == null)
			throw new ArgumentNullException("knots");
		int n = knots.Length - 1;
		if (n < 1)
			throw new ArgumentException
			("At least two knot points required", "knots");
		if (n == 1)
		{ // Special case: Bezier curve should be a straight line.
			firstControlPoints = new Point[1];
			// 3P1 = 2P0 + P3
			firstControlPoints[0].X = (2 * knots[0].X + knots[1].X) / 3;
			firstControlPoints[0].Y = (2 * knots[0].Y + knots[1].Y) / 3;

			secondControlPoints = new Point[1];
			// P2 = 2P1 – P0
			secondControlPoints[0].X = 2 *
				firstControlPoints[0].X - knots[0].X;
			secondControlPoints[0].Y = 2 *
				firstControlPoints[0].Y - knots[0].Y;
			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] = (8 * knots[n - 1].X + knots[n].X) / 2.0;
		// 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] = (8 * knots[n - 1].Y + knots[n].Y) / 2.0;
		// 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 : 3.5) - 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;
	}
}

虽然我用 C# 3.0 编译了这段代码,但我看不出为什么它不能在 C# 2.0 中不经任何修改地使用,甚至在 C# 1.0 中,如果您从类声明中移除关键字“static”也可以。

请注意,代码使用了 System.Windows.Point 结构,因此它旨在与 Windows Presentation Foundation 配合使用,并且需要 using System.Windows; 指令。但是,要将其与 Windows Forms 配合使用,您只需将 System.Windows.Point 类型替换为 System.Drawing.PointF。同样,如果需要,将此代码转换为 C/C++ 也很简单。

示例

本文提供的示例是一个面向 .NET 3.5 的 Visual Studio 2008 解决方案。它包含一个 WPF Windows 应用程序项目,旨在演示使用上述贝塞尔样条绘制的一些曲线。您可以从窗口顶部的组合框中选择一条曲线,尝试不同的点数,并设置适当的 XY 比例。您甚至可以添加自己的曲线,但这需要进行如下编码:

  1. 将您的曲线名称添加到 CurveNames 枚举中。
  2. 将您的曲线实现添加到 *Curves* 区域。
  3. OnRender 覆盖中添加对您的曲线的调用。

在示例中,我使用自定义 Canvas 上的 Path 元素来渲染曲线,但在实际应用程序中,您可能会使用更有效的方法,如视觉层渲染。

历史

  • 2008年12月18日:初次发布
  • 2009年3月23日:第二次文章修订,包含以下更正和补充
    1. 最重要的是 bug 修复。这个 bug 及其修正由 Peter Lee 发现。非常感谢他!这个 bug 在节点点数较少的情况下产生了视觉上可区分的行为。
    2. 如果传入无效参数,GetCurveControlPoints 现在会抛出异常。
    3. 添加了对节点数组只有两个点的情况的特殊处理。
    4. 添加了单元测试。
© . All rights reserved.