贝塞尔曲线与曲面(Bezier Curve and Surface)的详细介绍与代码实现

本文由 简悦 SimpRead 转码, 原文地址 zhuanlan.zhihu.com

前言

本文实在是有点长了,但是如果你坚持看完的话应该会对贝塞尔曲线和曲面有一个比较不错的了解(我自己是这样的)。同时学习完你还可以得到下面这样看着蛮有意思的效果,所以如果你有兴趣的话,那就坚持看完吧。

Demo 工程:

luckyWjr/Demo

图形学为什么要曲线和曲面

因为在生活中存在着各种各样光滑的曲线或曲面,例如汽车的表面,钢珠球等。

​在建模的时候,我们通常使用很多的小三角形来逼近这样的曲面,因此放大了看,就会发现这些面其实是凹凸不平的。但是实际生活中,例如我们身边的杯子,无论怎么看它都是一个光滑的表面。

因此我们想要来表达出这样一种光滑的曲线或曲面,那么我们应该怎么表达这些光滑的曲线或曲面呢?我们先从曲线入手,看看它的几种表达方式。

曲线的显示表示(Explicit representation)

学过函数我们知道, 画出来的图其实就是一条曲线,如下:

![](data:image/svg+xml;utf8,http://www.w3.org/2000/svg’ width=’280’ height=’316’>)

​即在画曲线前先通过向量绘制一个多边形,代表该曲线的趋势和走向。这就好比画画的时候先画个大致轮廓再画细节,雕刻的时候也是先雕个大致形状再慢慢打磨。

并且贝塞尔曲线具有交互性,也就是说我们可以通过修改向量,来修改曲线。

贝塞尔提出了如下公式用于计算曲线,将曲线表达成向量和基函数的乘积。

其中 代表的就是向量, 代表的是一个基函数,其内容如下:

注:该基函数本质上就是一个 n-1 次的多项式。

有了这个公式后,也就是说你给我一个多边形,我就可以用它算出一个曲线,从公式中我们也可发现,贝塞尔曲线属于一种参数形式(有参数 t)表示曲线的方式。该公式简单了解下就行了,不需要去理解它,据说贝塞尔本人都说不清这个公式,很有意思哈哈。

伯恩斯坦多项式

1972 年,Forrest 在 Computer Aided Design 杂志上发表了他的著名论文,在论文里他指出贝塞尔曲线可以借助伯恩斯坦多项式被定义在点集上

有关伯恩斯坦多项式的介绍可以参考,其中包括很多性质的推导:

王江荣:贝塞尔曲线中的伯恩斯坦多项式(Bernstein Polynomial)

原本贝塞尔曲线说的是向量相连,我们现在把这向量都变成一个个控制点,即给定控制顶点 ,贝塞尔曲线上的任意一点 可以定义为:

其中的 就是第 i 个 n 阶的伯恩斯坦多项式:

其中 为组合数(排列组合的那个组合),它的值为:

因为组合数有很多的性质,所以贝塞尔曲线的很多算法和定理,都是在弄这个组合公式。

伯恩斯坦多项式代码实现

有了公式我们不就可以写出代码来了吗,下面是 Unity 的简单示例:

  1. //伯恩斯坦多项式方法求曲线上一点
  2. Vector3 GetCurvePointByBernstein(Vector3[] referencePoint, float t)
  3. {
  4. uint i = 0;
  5. int n = referencePoint.Length - 1;// n 阶是 n+1 个顶点
  6. Vector3 point = Vector3.zero;
  7. for (; i <= n; i++)
  8. point += BernsteinNum(i, (uint)n, t) * referencePoint[i];
  9. return point;
  10. }
  11. //伯恩斯坦多项式
  12. float BernsteinNum(uint i, uint n, float t)
  13. {
  14. return CombinatorialNum(i, n) * Mathf.Pow(t, i) * Mathf.Pow(1.0f - t, n - i);
  15. }
  16. //组合数
  17. ulong CombinatorialNum(uint i, uint n)
  18. {
  19. if (i > n) return 0;
  20. return Factorial(n) / (Factorial(i) * Factorial(n - i));
  21. }
  22. //阶乘
  23. ulong Factorial(uint n)
  24. {
  25. if (n == 0) return 1;
  26. ulong result = n;
  27. for (n--; n > 0; n--)
  28. result *= n;
  29. return result;
  30. }

注:这里要注意的是公式里的 n 并不是说 n 个顶点,而是 0 到 n 一共 n+1 个顶点。还有代码就临时简单写的,也没优化,要是写的不好看别见怪!懂公式的话,很容易就能看懂。

这样我们就可以求出 t 从 0 到 1 所有曲线上的点,这些点就可以连成一条贝塞尔曲线。

德卡斯特里奥算法(de Casteljau Algorithm)

但是实际上,我们并不会使用上面那个算法来求曲线上的任意点 ,而是使用一种名为 de Casteljau 的递归算法来求。该算法比之前的方法慢,但在数值上更为稳定

接下来我们来看看用 de Casteljau 算法怎么求出曲线上的一点的。如下图,我们先随便选取三个控制点,并将它们前后连线。

然后我们在 线段上通过线性插值在线段上找到点 ,使得 ,其他线段也如此,我们就可得到下图

​接着我们连接 ,得到新的线段,然后在该线段上再取一点使得该线段被分为 t 和 1-t,那么就会得到下图:

​此时已经不能再连线了,而我们得到的点 就是这三个控制点对应的贝塞尔曲线在 t 位置上的点。

这样我们就可以通过使 t 从 0 变到 1,得到所有曲线上的点,从而得到曲线。这样的求曲线上任意一点方式也就是前面所说的 de Casteljau 算法

如果有更多的控制点,我们也可以使用相同的方法来求出曲线上的一点,如下图是四个控制点求曲线上一点的过程:

​对于三个顶点控制的贝塞尔曲线我们称之为二阶贝塞尔曲线(Quadratic Bezier),那么四个顶点自然是三阶贝塞尔曲线(Cubic Bezier),因此 n 阶贝塞尔曲线有 n+1 个顶点

而 de Casteljau 算法则是把这 n+1 个点,通过参数 t 用线性插值的方法,先变成 n 个新的点,然后在变成 n-1 个新的点,递归下去,最终得到一个点,这个点就是贝塞尔曲线上的点。

如下动图,很直接明了的表示了贝塞尔曲线的绘制过程:

de Casteljau 算法代码实现

该算法看起来就和我们伯恩斯坦多项式没有任何关系了,我们只需要写个递归的方法,方法体内求线性插值即可:

  1. //de Casteljau 算法求曲线上一点
  2. Vector3 GetCurvePointDeCasteljau(Vector3[] referencePoint, float t)
  3. {
  4. int len = referencePoint.Length;
  5. if (len < 2)
  6. Debug.LogError("控制点至少需要两个");
  7. Vector3[] newPoint = new Vector3[len - 1];
  8. for (int i = 0; i < len - 1; i++)
  9. newPoint[i] = Vector3.Lerp(referencePoint[i], referencePoint[i + 1], t);
  10. if (len == 2)
  11. return newPoint[0];
  12. else
  13. return GetCurvePointDeCasteljau(newPoint, t);
  14. }

伯恩斯坦多项式与 de Casteljau 算法的关系

回到公式部分,我们来看看伯恩斯坦多项式与 de Casteljau 算法的关系。我们拿最简单的二阶贝塞尔曲线举例,如下图:

​图中蓝色的点为控制点,他们的坐标我们是知道的,那么通过线性插值,我们可以得到求出红色点的坐标,公式如下( 我们用 来代替):


红色点坐标求出后,我们自然可以再求出绿色点的坐标:

把上面两个式子带入到下面的式子,得到:

我们还可以用这个方法去算三阶的,四阶的,乃至 n 阶的贝塞尔曲线,得到的结果为曲线上任意一点 P(t) 是各个顶点的线性组合,即:

而我们每个顶点前面的系数 k,就是伯恩斯坦多项式。例如二阶贝塞尔曲线对应的伯恩斯坦多项式为 ,其中 ,正好对应前面的三个系数。

因此可以得出结论,对于 n 阶的贝塞尔曲线,曲线上 t 位置上的点 P(t) 的坐标是由 n+1 个顶点和伯恩斯坦多项式的乘积求和

上面式子也就是 Forrest 在 1972 年提出的结论。

前面我们是从线性插值计算,逆推到伯恩斯坦多项式。现在我们来看看怎么直接使用伯恩斯坦多项式得到递归的结果。

注:下面公式看着都很长,其实都很简单,既没有微积分也没啥奇怪的公式,就是简单的加减乘除与一个插值的概念,所以看到它们的时候不要怕!

在讲到伯恩斯坦多项式的时候,我们说它具有递归性,即可以把 n 阶的伯恩斯坦多项式写成 n-1 阶的多项式组合:

注:后面会写到很多伯恩斯坦多项式的性质公式,想了解的话还是建议看下之前的文章。

也就是说原本的方程式:

可以写成:

展开来,可得:

我们会发现里面有 ,我们提取一下不就变成了 ,而 是啥子?不就是两点的线性插值么。对于后面的几项也是一样的,上面式子就可以写成:

那么就可以把贝塞尔曲线的方程式写成:

也就是说我们把原本 n 个控制点,通过 变成了 n-1 个新的控制点。那么 n-1 个控制点又可以按这个方法变成 n-2 个控制点,一直递归下去,最终只剩一个控制点,也就是曲线上的点。这个方法也正是我们之前贝塞尔曲线的绘制过程。

贝塞尔曲线的代码实现

简单用 Unity 写了一下:

  1. public class BezierCurve : MonoBehaviour
  2. {
  3. public Transform[] referencePointTransforms;
  4. public Transform curvePointTransform;
  5. Vector3[] mReferencePoints;
  6. int mSign = 1;
  7. float t = 0;
  8. int mResolution = 100;
  9. void Update()
  10. {
  11. TransformsToVectors();
  12. t = Mathf.Clamp(t, 0, 1);
  13. // curvePointTransform.position = GetCurvePointDeCasteljau(mReferencePoints, t);
  14. curvePointTransform.position = GetCurvePointByBernstein(mReferencePoints, t);
  15. t += Time.deltaTime * 0.2f * mSign;
  16. if (t > 1 || t < 0)
  17. mSign *= -1;
  18. }
  19. void OnDrawGizmos()
  20. {
  21. TransformsToVectors();
  22. if (mReferencePoints.Length == 0)
  23. return;
  24. //画控制点连线
  25. Gizmos.color = Color.red;
  26. for (int i = 0; i < mReferencePoints.Length - 1; i++)
  27. Gizmos.DrawLine(mReferencePoints[i], mReferencePoints[i + 1]);
  28. //画曲线
  29. Gizmos.color = Color.green;
  30. Vector3 start = mReferencePoints[0];
  31. for (int i = 1; i <= mResolution; i++)
  32. {
  33. Vector3 end = GetCurvePointDeCasteljau(mReferencePoints, (float)i / mResolution);
  34. Gizmos.DrawLine(start, end);
  35. start = end;
  36. }
  37. }
  38. void TransformsToVectors()
  39. {
  40. int referencePointCount = referencePointTransforms.Length;
  41. mReferencePoints = new Vector3[referencePointCount];
  42. for (int i = 0; i < referencePointCount; i++)
  43. mReferencePoints[i] = referencePointTransforms[i].position;
  44. }
  45. }

效果图如下:

![](data:image/svg+xml;utf8,http://www.w3.org/2000/svg’ width=’397’ height=’310’>)

端点性质

伯恩斯坦多项式具有端点性质,贝塞尔曲线同样具有,其实从绘制过程中,我们就很容易看出来,曲线是从第一个控制点开始到最后一个控制点结束的,即:


切向量(Tangent Vector)

切向量也就是对贝塞尔曲线求导,我们知道伯恩斯坦多项式的求导公式如下:

那么贝塞尔曲线的求导结果即为:

我们先来看看 t=0 的情况,可得:

由于伯恩斯坦多项式的端点性质:

所以 都为 0,可得:

同理也可证明,当 t=1 时:

也就是说贝塞尔曲线在起点处的切向量方向就是 的方向,终点处的切向量方向就是 的方向,如下图:

​通过公式,也很好理解为什么三阶的贝塞尔曲线起点和终点的切向量前面有个系数 3 了。

二阶导数

同样的,我们还可以求贝塞尔曲线的二阶导数,公式如下(这里就不推导了,简单的记个笔记):

因此:


有了一阶和二阶导数,我们就可以算出曲率,曲率是微分几何里面描述曲线弯曲程度的概念。通过曲率公式可得:


k 阶导数

贝塞尔曲线的 k 阶导数的差分形式如下:

代表着 k 阶差分,定义如下:

对称性

因为伯恩斯坦多项式有对称性,因此贝塞尔曲线同样拥有对称性。也就是说我们把 这些顶点顺序倒过来,求得的曲线和之前的曲线一模一样,只是方向相反,如下图:

凸包性(Convex Hull)

根据伯恩斯坦多项式的归一性,我们可以得出贝塞尔曲线的凸包性,也就是说贝塞尔曲线必定在所有控制点的凸包内部。

所谓凸包即是可以包含所有顶点的最小凸多边形,如下图:

![](data:image/svg+xml;utf8,http://www.w3.org/2000/svg’ width=’353’ height=’93’>)

我们以 AB 两点作为控制点形成贝塞尔曲线,那么该曲线自然就是 AB 线段,同样的,我们以 BC 两点再形成一个贝塞尔曲线为 BC 线段,AB 和 BC 形成一条直线段。此时 B 点左右的导数分别为 (B-A) 和(C-B),这两个值明显是不相等的,按微积分里的定义就是不连续。可是它明明是一条直线,怎么会不连续呢?这个反例意味着传统的连续性不适用于描述 CAD 和图形学中形状的连续性,因此提出了几何连续性(Geometric continuity)的概念。

所谓几何连续性的概念,就是导数的方向一样,我们即可认为是连续的。对于例子中的直线,我们认为它是几何连续的,但是速度不一样,例如一个小蚂蚁在 ABC 上爬,在 AB 上较慢,在 BC 上变快了。

那么我们就有如下定义,假设我们有 形成的两条贝塞尔曲线,如果:

,我们称之为 0 阶连续,标记为
,且这四个点共线,我们称之为 1 阶几何连续,如果 k=1,那么就是 1 阶连续,标记为
如果二阶导数相等,即曲率连续,我们称之为 2 阶连续,标记为

连续性的效果如下(http://math.hws.edu/eck/cs424/notes2013/canvas/bezier.html):

贝塞尔曲面

我们接下来看看下面这种情况:

如图,我们有四条一模一样的贝塞尔曲线,图中的小红点是这些曲线在 t 从 0 到 1 变化时对应的点。

那么假如我们现在把这些小红点互相连线,会形成啥呢?想象一下,没错,就是一个曲面!如下图:

噢!这该死的摩尔纹!!!

但是这样形成的曲面其实有很大的问题,它仅仅只是在一个方向上弯曲,在与其垂直的方向上依旧还是笔直的,如果我们简单修改下控制点的位置,就会出现下面这样的问题:

那么我们就要想办法优化它,人们想到,如果用这个四个小红点再做控制点,不就可以在垂直方向上又形成一条光滑的曲线了么?

没错,这就是贝塞尔曲线的思路,我们先用一组贝塞尔曲线在 t1 位置得到一组对应在曲线上的点,然后将这些点再作为控制点来做曲线,所有的曲线组合起来,即可得到一个光滑的曲面。

效果图如下:

图中可以发现,不管我们怎么修改控制点,都可以得到在所有方向上都是光滑的曲线。

害,要不是因为曲面的示意图少,我也不用专门写代码来展示了!!!瞎鸡儿写一下的代码如下:

  1. public class BezierSurface : MonoBehaviour
  2. {
  3. public BezierCurve[] curves;
  4. int mResolution = 100;
  5. void OnDrawGizmos()
  6. {
  7. Gizmos.color = Color.green;
  8. for (int i = 1; i <= mResolution; i++)
  9. {
  10. //获取每个曲线上 t 时刻的点
  11. Vector3[] curvePoints = new Vector3[curves.Length];
  12. for (int j = 0; j < curvePoints.Length; j++)
  13. {
  14. Vector3[] referencePoint = GetCurveReferencePoint(curves[j]);
  15. curvePoints[j] = GetCurvePointDeCasteljau(referencePoint, (float)i / mResolution);
  16. }
  17. //根据上面获得的点再画曲线
  18. Vector3 start = curvePoints[0];
  19. for (int j = 1; j <= mResolution; j++)
  20. {
  21. Vector3 end = GetCurvePointDeCasteljau(curvePoints, (float)j / mResolution);
  22. Gizmos.DrawLine(start, end);
  23. start = end;
  24. }
  25. // 那些点直接连线的错误示范
  26. // Vector3 start = curvePoints[0];
  27. // for (int j = 1; j < curvePoints.Length; j++)
  28. // {
  29. // Vector3 end = curvePoints[j];
  30. // Gizmos.DrawLine(start, end);
  31. // start = end;
  32. // }
  33. }
  34. }
  35. Vector3[] GetCurveReferencePoint(BezierCurve curve)
  36. {
  37. Vector3[] array = new Vector3[curve.referencePointTransforms.Length];
  38. for (int i = 0; i < array.Length; i++)
  39. array[i] = curve.referencePointTransforms[i].position;
  40. return array;
  41. }
  42. }

![](data:image/svg+xml;utf8,http://www.w3.org/2000/svg’ width=’400’ height=’224’>)

原理知道了,接下来我们看看怎么求贝塞尔曲面上任意一点,以及贝塞尔曲面的一些性质。

又到了密密麻麻又枯燥的公式部分了(其实公式一点也不枯燥,专研起来很有趣的,就是看着确实很头疼)。

贝塞尔曲面上的任意一点

首先我们有一个 n 阶的贝塞尔曲线,然后我们给它复制 m+1 份,那么在每个贝塞尔曲线的 t1 位置,我们可以得到 m+1 个曲线上的点。我们再把这些点作为控制点,可以得到一个 m 阶的贝塞尔曲线,而这个贝塞尔曲线就是最终贝塞尔曲面上的一条线,因此该曲线 t2 位置上的点就是贝塞尔曲面上的点。其中 t1 和 t2 通常会使用 u 和 v 来代替,它们的取值范围自然都是 0 到 1,那么对应贝塞尔曲面上的任意一点,我们就可以用 P(u,v) 来表示。

对于上面这种描述得到的曲面,我们称之为 m*n 阶的贝塞尔曲面,即有 m+1 条 n 阶的贝塞尔曲线,一共会得到 (n+1)*(m+1) 个顶点,这里我们定义第一个 n 阶的贝塞尔曲线的第一个控制点为 ,该曲线往后的控制点为 ,那么第 i 条 n 阶的贝塞尔曲线的第 j 个顶点就是 ,最后一条的最后一个顶点自然是 了。

要求贝塞尔曲面上的一点,我们先来复习一下贝塞尔曲线上一点的公式(毕竟文章稍微有点长了,在往上翻去找有点费劲):

那么第 i 条 n 阶的贝塞尔曲线上的点 P(u) 即为:

那么我们就可以得到 m 个 P(u),我们可以标记为 ,这些 P(u) 会再形成贝塞尔曲线,也就是曲面上的曲线,那么该曲线上的点 P(v) 即为:

那么两个公式合起来,就是我们求贝塞尔曲面上任意点 P(u,v) 的公式了:

此外,我们还可以用矩阵的形式来表示,如下:

简单的解读下,首先是 1*(n+1) 的矩阵乘以 (n+1)(m+1) 的矩阵,得到的结果是 1(m+1) 的矩阵,这一步其实就是在求 m+1 条贝塞尔曲线 u 位置的顶点。然后把 1*(m+1) 的矩阵再乘以 (m+1)*1 的矩阵,得到的结果就是一个值,这一步就是求新的曲线在 v 位置上的点。

当然了,同样的,我们也可以用 de Casteljau 算法来算曲面上的任意一点的坐标

端点性质

控制网格的四个顶点同样是贝塞尔曲面的顶点,即:




切平面

以下三角形指定了贝塞尔曲面四个角的切平面:




此外,贝塞尔曲线带有的凸包性,对称性,几何不变性,在贝塞尔曲面中也都有这些性质。

连续性

首先要两个贝塞尔曲面连续,它们的控制顶点数必须相同,不然一个 9 个点的曲面和一个 16 个点的曲面,很难做到连续。因此对与连续性,我们只考虑两个曲面控制顶点相同的情况。

假设有两个 m*n 阶的贝塞尔曲面,它们的控制顶点分别是

那么 0 阶连续( 的表示如下:

i = 0,1,…,m

这个式子很好理解,就是我们第一个贝塞尔曲面的每条贝塞尔曲线的最后一个控制点要和另一个贝塞尔曲面的每条贝塞尔曲线的第一个控制点重叠。

每个控制点都一样,那么得到的曲线也都一样,所以曲线上对应 v 位置的点也都一样,所以也可以写成如下形式:

我们用之前的烂代码在 Unity 模拟一下,如下图,是两个 3*3 阶的贝塞尔曲面,我们使头尾控制点位置相同,图中圈起来的地方

然后我们不管怎么移动其他控制点,链接处都保持这一样的曲线,这种就是所谓的 0 阶连续

而对于 1 阶连续(,同样是对切向量进行限制。为了方便理解,我们举个例子,如下图:

红点是我们两个曲面连接的控制点,我们来看 这个点。首先我们知道这个点会有如下两个切线,一个是 方向的,一个是 方向的。其次由于小红点又形成了贝塞尔曲线,因此在小红点这个方向上 还有个切线,即 。因为四个小红点本身形成的就是一个光滑连续的曲线,所以在这个方向上,一个点不会存在左右两个不同的切线。所以在两个曲面的连接处上的任意一点,我们其实可以得到三个切线。而只要这三个切线共面,我们就称之为一阶几何连续。如果除了小红点方向的切线,另外两个切线的值相等方向相反的话,那么就是一阶连续。

三角域贝塞尔曲面

前面介绍的贝塞尔曲面我们称之为矩形域贝塞尔曲面,u 和 v 都是 0 到 1,定义在矩形上。接下来我们介绍一下定义在三角形上的贝塞尔曲面,也就是三角域贝塞尔曲面

没错,又要用到重心坐标了,再次发下自己的之前的链接,嘿嘿

王江荣:直线与三角形的重心坐标(Barycentric Coordinates)的推导与应用

我们知道三角形上任意一点,我们可以用 (u,w,w) 的坐标来表示,其中 u+v+w=1,这个表示方式就是重心坐标。

如下图,如果我们只有如下三个顶点,那么每条边自然都是直线,形成的也就是个三角形平面。

![](data:image/svg+xml;utf8,http://www.w3.org/2000/svg’ width=’346’ height=’304’>)

三条边变曲线好搞,我们要怎么让整个三角形面都变成光滑的曲面么?一样的思路,我们要求出三角形内任意一点的坐标,利用重心坐标,我们对三角形内任意点记作 P(u,v,w) 。

然后我们再定义一些规范,首先每条边的顶点数要相同,不能一条边 4 个顶点,一条边 10 个顶点,这不是为难人嘛。如果每条边都有 n+1 个顶点,我们称之为 n 阶的三角域贝塞尔曲面

为了方便理解,我们首先假设每条边被各自的顶点分成了等长的线段,那么我们就可以求出每个顶点的重心坐标了(至于哪个顶点对应 u,哪个对应 v,哪个对应 w,随你喜欢),如下:

得到重心坐标后,我们把这些值同时乘以阶数 n ,得到的值就作为控制点的下标。例如上图中三角形是 3 阶的,那么就把这些重心坐标的值都乘以 3,例如重心坐标为 (1/3, 0, 2/3) 的控制点,乘以 3 后为(1,0,2),那么该控制点就标记为 ,这样我们可以得到下图:

![](data:image/svg+xml;utf8,http://www.w3.org/2000/svg’ width=’364’ height=’313’>)

代码实现

有了公式我们就可以写代码了(依旧简单的写一下,写的不好别见怪)

  1. [Serializable]
  2. public struct TriangleReferencePoint
  3. {
  4. public string ijk;
  5. public Transform trans;
  6. public bool isEqual(uint i, uint j, uint k)
  7. {
  8. uint num = uint.Parse(ijk);
  9. if (100 * i + 10 * j + k == num)
  10. return true;
  11. return false;
  12. }
  13. }
  14. public class TriangleBezierSurface : MonoBehaviour
  15. {
  16. public uint n = 3;
  17. public TriangleReferencePoint[] referencePoints;
  18. int mResolution = 30;
  19. Vector3[,,] PointsToVector3s()
  20. {
  21. Vector3[,,] array = new Vector3[n + 1, n + 1, n + 1];
  22. for (uint i = 0; i <= n; i++)
  23. {
  24. for (uint j = 0; j <= n-i; j++)
  25. {
  26. uint k = n - i - j;
  27. array[i, j, k] = FindReferencePoint(i, j, k).trans.position;
  28. }
  29. }
  30. return array;
  31. }
  32. TriangleReferencePoint FindReferencePoint(uint i, uint j, uint k)
  33. {
  34. foreach (var point in referencePoints)
  35. {
  36. if (point.isEqual(i, j, k))
  37. return point;
  38. }
  39. Debug.LogError("没找到对应控制点,错误的控制点下标");
  40. return new TriangleReferencePoint();
  41. }
  42. void OnDrawGizmos()
  43. {
  44. if ((n + 1) * (n + 2) / 2 != referencePoints.Length)
  45. {
  46. Debug.LogError("顶点数不匹配");
  47. return;
  48. }
  49. Gizmos.color = Color.green;
  50. Vector3[,,] referencePointPositions = PointsToVector3s();
  51. float u, v, w;
  52. uint i, j, k;
  53. float step = 1.0f / mResolution;
  54. //取重心坐标
  55. for (u = 0; u <= 1.0f; u += step)
  56. {
  57. for (v = 0; v <= 1.0f - u; v += step)
  58. {
  59. w = 1.0f - u - v;
  60. //利用三角域的伯恩斯坦多项式,计算每个(u,v,w)对应的曲面上的点P(u,v,w)
  61. Vector3 p = Vector3.zero;
  62. for (i = 0; i <= n; i++)
  63. {
  64. for (j = 0; j <= n - i; j++)
  65. {
  66. k = n - i - j;
  67. p += TriangleBernsteinNum(i, j, k, n, u, v, w) * referencePointPositions[i, j, k];
  68. }
  69. }
  70. Gizmos.DrawSphere(p, 0.02f);
  71. }
  72. }
  73. }
  74. //三角域的伯恩斯坦多项式
  75. float TriangleBernsteinNum(uint i,uint j,uint k, uint n, float u, float v, float w)
  76. {
  77. return TriangleCombinatorialNum(i, j, k, n) * Mathf.Pow(u, i) * Mathf.Pow(v, j) * Mathf.Pow(w, k);
  78. }
  79. //三角域的组合数
  80. ulong TriangleCombinatorialNum(uint i,uint j,uint k, uint n)
  81. {
  82. if (i + j + k != n) return 0;
  83. return Factorial(n) / (Factorial(i) * Factorial(j) * Factorial(k));
  84. }
  85. //阶乘
  86. ulong Factorial(uint n)
  87. {
  88. if (n == 0) return 1;
  89. ulong result = n;
  90. for (n--; n > 0; n--)
  91. result *= n;
  92. return result;
  93. }
  94. }

](data:image/svg+xml;utf8,<svg xmlns='<a href=http://www.w3.org/2000/svg’ width=’159’ height=’191’></svg>)![” title=”” />

我们来看看得到的效果:

很骚有没有!!!

三角域贝塞尔曲面的 de Casteljau 算法

前面我们知道三角域的贝塞尔曲面上的任意一点定义如下:

我们把它展开来可以得到下面式子(注:为了简洁,就省略伯恩斯坦多项式里带的 (u,v,w) 了):

然后我们来看看三角域伯恩斯坦多项式的递归公式:

带入上面公式,得到:

我们会发现里面有这么三项相加:

我们提取一下不就变成了:

括号里的不就是 三角形内的一点嘛。

我们最终可以得出,如果 n 阶的三角形(其控制点为 )变成 n-1 阶的三角形(其控制点为 ),那么有如下关系:

这样一直递归下去,最终会得到只剩一个点,那么这个点就是这个三角域贝塞尔曲面上的点。

有点难懂?我们来看下示意图,前面我们有一个 3 阶的三角形,那么降阶后就会变成一个二阶的三角形,如下图,红色点就是降阶后得到的新的控制点:

![](data:image/svg+xml;utf8,http://www.w3.org/2000/svg’ width=’332’ height=’280’>)

目前个人对贝塞尔曲线和曲面了解的就这么多了~ 谢谢阅读。此外如今还有更加牛逼的 B - 样条曲线(B-spline curve),当然也更加的复杂。

贝塞尔曲线中的伯恩斯坦多项式(Bernstein Polynomial)
光栅化与深度缓存