计算机图形学十贝塞尔曲线

錶哥

已于 2023-08-08 02:21:08 修改

阅读量1.9k

点赞数 5

分类专栏：计算机图形学文章标签： c++

于 2022-12-20 16:14:15 首次发布

本文链接：https://blog.csdn.net/weixin_44683202/article/details/128379717

版权

计算机图形学专栏收录该内容

23 篇文章 6 订阅

订阅专栏

贝塞尔曲线

我们在建模时，一般使用小三角形来逼近这样的曲面，因此放大了看，这些表面很多都是凹凸不平的。但是生活中的很多面，比如汽车，比如玻璃球，都是一个光滑的表面。
因此我们想要来表达出一种光滑的曲线或者说曲面。

曲线的表现形式分为两种：显式表现形式和隐式表现形式。

显式表现例如在这里插入图片描述
通俗的说，显示表现就是，以自变量x的值来表达因变量y的值
曲线的隐式表示的表现形式，就是利用隐式方程来刻画一条曲线，在二维空间中一个隐式的曲线可以通过 f(x,y)=0 来表达。
其中，当f的值小于0时，该点一般在曲线围成的区域的内部，当大于0时，该点一般在曲线围成的区域的外部。

曲线在计算机图形学当中最常用的表示形式是参数形式表示（Parametric Form）。一个曲线的参数形式是通过一个自变量（参数）t来表达曲线上每个点的空间坐标，它也可算作是显示表示的一种。
在三维空间中，我们可以用如下的三个显示的函数，来表示一条空间曲线。
在这里插入图片描述
t=0，代表曲线的起点，t=1代表曲线的终点，因此可以求出t（0<=t<=1）代表曲线上的任意一点P（t），然后将t的值带入三个显式的方程，即可求出曲线上点P的坐标：

参数曲线并不是唯一的，一个给定的曲线或者曲面可以通过不同的形式来表达。给定一个表达式我们可以画出一个唯一的曲线，但是该曲线可能还有别的表达式可以来表示。
参数形式有一个好处。比如上面的式子表达的是三维空间的曲线，如果我们删除了z=z(t)，那么它就变成了一个二维空间的曲线，即降维很方便。并且同样的也容易推广到高维，新增新维度对应的显示方程即可。

我们可以对每个点进行求导，即
在这里插入图片描述
得到的结果被视为曲线的绘制速度，该导数指向曲线的切向。

举一个二维函数的例子，比方说我有一个在这里插入图片描述
其中x的取值范围为1-4，那么这个曲线如下图：
我们来看看怎么用参数的形式来对其进行表示，首先是x和t的关系。因为x是1-4，而t是0-1，因此我们可以得到方程x=3t+1。然后是y，因为我们知道这个曲线y和x的关系，因此我们只需要将x关于t的表示带入即可。
所以我们上述的曲线形式可以表示为：
在这里插入图片描述
通过带入不同的t值，即可得到该曲线的各个点的坐标，例如P(0)=(1，1.3)，P(1)=(4，5.8)。

参数形式与多项式

我们曲线的参数形式写成一个多项式的形式。例如我们前面的二维例子，我们可以把它写成一个二阶的多项式，即在这里插入图片描述
的形式。但是其中a,b,c并不是常数，而是代表着向量。其结果为：

这个成为样条曲线。

贝塞尔曲线（Bezier Curve）

相对之前的几种曲线表示方法，贝塞尔提出的是一种通过连接向量（connected vectors）来表示曲线的方法，如下：
在这里插入图片描述
即在画曲线前先通过向量绘制一个多边形，代表该曲线的趋势和走向。这就好比画画的时候先画个大致的轮廓再画细节，雕刻的时候也是先雕个大致形状再慢慢打磨。
并且贝塞尔曲线具有交互性，也就是我们可以通过修改向量，来修改曲线。
贝塞尔提出了如下公式来计算曲线，将曲线表达成向量和基函数的乘积。
在这里插入图片描述
其中
代表的是向量，
代表的是一个基函数，其内容如下：

注：该基函数的本质上就是一个n-1次的多项式。
有了这个公式后，也就是说你给我一个多边形，我就可以用它算出一个曲线，从公式中我们可以发现，贝塞尔曲线属于一种参数形式（有参数t）表示曲线的方式。
该公式简单了解一下即可。

伯恩斯坦多项式

1972年，Forrest在Computer Aided Design杂志上发表了他的著名论文，在论文里他指出贝塞尔曲线可以借助伯恩斯坦多项式被定义在点集上。
原本贝塞尔曲线说的是向量相连，现在我们把一个个向量当做控制点，即给定控制顶点P0，P1，P2、、、、、Pn，贝塞尔上任意一点P(t)可以定义为：在这里插入图片描述
其中的就是第i个n阶的伯恩斯坦多项式：

其中，就是高中数学中学过的排列组合中的组合数，值为

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

伯恩斯坦多项式代码实现

//伯恩斯坦多项式方法求曲线上一点
Vector3 GetCurvePointByBernstein(Vector3[] referencePoint, float t)
{
    uint i = 0;
    int n = referencePoint.Length - 1;// n 阶是 n+1 个顶点
    Vector3 point = Vector3.zero;
    for (; i <= n; i++)
        point += BernsteinNum(i, (uint)n, t) * referencePoint[i];
    return point;
}
//伯恩斯坦多项式
float BernsteinNum(uint i, uint n, float t)
{
    return CombinatorialNum(i, n) * Mathf.Pow(t, i) * Mathf.Pow(1.0f - t, n - i);
}
//组合数
ulong CombinatorialNum(uint i, uint n)
{
    if (i > n) return 0;
    return Factorial(n) / (Factorial(i) * Factorial(n - i));
}
//阶乘
ulong Factorial(uint n)
{
    if (n == 0) return 1;
    ulong result = n;
    for (n--; n > 0; n--)
        result *= n;
    return result;
}

德卡斯特里奥算法（de Casteljau Algorithm）

但是实际上，我们不会用上面的算法来求曲线上的任意点P(t)，而是使用一种名为deCasteljau的递归算法来进行求解。该算法比之前的方法慢，但是在数值上更为稳定。
接下来我们看看用deCasteljau算法怎么求出曲线上的一点的。如下图，我们先随便选取三个控制点，并将他们前后进行连线。
在这里插入图片描述
然后我们在P0P1线段上通过线性插值在线段上找到点P0,1，使得P0P0,1：P0,1P1 = t:1-t，其他线段也是如此，因此我们就可以得到下图：

在这里插入图片描述
接着我们连接P01P11，得到新的线段，然后在该线段上再取一点使得该线段被分为t和1-t，那么我们就会得到下图：

此时已经不能再进行连线了，而我们得到的点P0,2就是这三个控制点对应的贝塞尔曲线在t位置上的点。
这样我们就可以通过使t从0变到1，得到曲线上所有的点，从而得到曲线。这样的求曲线上任意一点的方式就是前面所说的de Casteljau算法。
如果有更多的控制点，我们也可以使用相同的方法来求出曲线上的一点，如下图是四个控制点求曲线上一点的过程：

在这里插入图片描述
对于三个顶点控制的贝塞尔曲线我们称之为二阶贝塞尔曲线，那么四个顶点自然是三阶贝塞尔曲线，因此n阶贝塞尔曲线有n+1个顶点。

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

在这里插入图片描述
该算法看起来就和伯恩斯坦多项式没有什么大关系了，因此我们只需要写个递归的方法，方法体内求线性插值即可：

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

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

回到公式部分，我们来看看伯恩斯坦多项式与de Casteljau算法的关系。我们拿最简单的二阶贝塞尔曲线举例，如下图：
在这里插入图片描述
图中蓝色点为控制点，他们的坐标我们是知道的，那么通过线性插值，我们可以得到求出红色点的坐标，公式如下：

红色点坐标求出来之后，我们自然可以再求出绿色点的坐标。

把上面两个式子带入到下面的式子，得到：
在这里插入图片描述
我们还可以用这个方法去算三阶的，四阶的，乃至n阶的贝塞尔曲线，得到的结果为曲线上任意一点P(t)是各个顶点的线性组合，即：

而我们每个顶点前面的系数k，就是伯恩斯坦多项式。例如二阶贝塞尔曲线对应的伯恩斯坦多项式为在这里插入图片描述，
其中，
正好对应前面的三个系数。

因此可以得出结论，对于n阶的贝塞尔曲线，曲线上t位置上的点P(t)的坐标是由n+1个顶点和伯恩斯坦多项式的乘积求和：
在这里插入图片描述
前面我们是从线性插值计算，逆推到伯恩斯坦多项式。现在我么来看看怎么使用伯恩斯坦多项式得到递推的结果。

在讲到伯恩斯坦多项式时，我们说它具有递归性，即可以把n阶伯恩斯坦多项式写成n-1阶的多项式组合：
在这里插入图片描述
原本的方程式：

可以写成：

展开来，可得：
可以化简把贝塞尔曲线写为：

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

端点性质：

伯恩斯坦多项式具有端点性质，贝塞尔曲线同样具有。曲线是从第一个控制点开始到最后一个控制点结束的，即：
在这里插入图片描述

切向量（tangent vector）

切向量也就是对贝塞尔曲线求导，我们知道伯恩斯多项式的求导公式如下：
在这里插入图片描述
那么贝塞尔曲线的求导结果即为：

我们先来看看t=0的情况，可得：

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

所以
而
都为0，可得：

同理也可证明，当 t=1 时：

也就是说贝塞尔曲线在起点处的切向量方向就是P0到P1的方向，终点处的切向量方向就是在这里插入图片描述到的方向，如下图：

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

二阶导数

同样的，我们还可以求贝塞尔曲线的二阶倒数，公式如下
在这里插入图片描述
因此：

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

k阶导数

贝塞尔曲线的k阶导数的差分形式如下：
在这里插入图片描述
代表着k阶差分，定义如下：

对称性

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

凸包性

根据伯恩斯坦多项式的归一性，我们可以得出贝塞尔曲线的凸包性，也就是说贝塞尔曲线必定在所有控制点的凸包内部。
所谓凸包即是可以包含所有顶点的最小凸多边形，如下图：
在这里插入图片描述
黑色的点是我们的控制点，那么蓝色的形状就是这些控制点的凸包。我们可以把这些控制点想象成木板上的钉子，蓝色的线是一根可拉到无限大的橡皮筋，当我们把橡皮筋拉到包含所有顶点时，一松手，橡皮筋收缩后形成的形成就是凸包。
在这里插入图片描述

几何不变性

即贝塞尔曲线的一些几何性质不随坐标系变换，因为贝塞尔曲线的位置和性质依赖于控制点，而不是坐标系。

仿射变换

对贝塞尔曲线做仿射变换，我们可以通过先对控制点做仿射变换，然后利用变换或的控制点来得到新的曲线，这个曲线和直接对曲线做仿射变换得到的结果是一样的。

但是需要注意，透视投影变换不能这么操作。

升阶操作

所谓升阶操作，就是比如我们原本由n+1个控制点得到一个n阶贝塞尔曲线，我们可以把它变成是n+2个新控制点并且得到n+1阶的贝塞尔曲线，这两个曲线一模一样。

升阶操作的好处在于可以增加灵活性，例如我们原本有一个二阶的贝塞尔曲线，想要改变曲线形状只能修改三个顶点的位置。但是通过升阶，把它变成三阶的四阶的乃至更高阶的，我们就可以通过更多的顶点来调整这个曲线。

伯恩斯坦多项式的升阶公式如下：
在这里插入图片描述
那么就可以写成如下形式：

展开来看：

经过一系列操作，我们新的式子就可以写成：

也就是说头尾两个控制点不变，中间生成n个新的控制点，这样原本n+1个控制点就变成了n+2个控制点，实现了升阶操作。

升阶公式最终如下：
在这里插入图片描述
当升阶的次数越多，顶点形成的线条就会越逼近曲线本身。在微积分里有个Weierstrass逼近定理：一个连续函数总是可以用一个多项式无限逼近，我们的升阶操作就等于证明了这个定理（这定理本质上就可以用伯恩斯坦多项式来证明）。

但是这里其实有个问题，我们之前说切向量的时候，说贝塞尔曲线顶点和终点的切向量是第一个和最后一个控制点的方向。但是如上图，我们曲线没有变，等于切向量没有变。但是随着控制点增加，控制点的方向一直在变，这不是违背了吗？实际上我们将一个n阶贝塞尔曲线升阶，不管升多少次，本质上（代数上）还是n阶的，升阶只是给它换了个表现形式而已。

降阶操作

降阶的目的是寻找一组新的控制点定义的曲线，使得误差最小。

降阶是升阶的逆过程，我们假设n阶贝塞尔曲线是n-1阶贝塞尔曲线的升阶结果，根据根据升阶公式我们可以得到：
在这里插入图片描述
那么我们就可以得到下面两个递推公式（下面的公式1用Pi表示降阶后的顶点，公式用Pi）：

这两个公式问题都很多，误差很大，公式1在P0处是精确的，越往后误差越大，而公式2在Pn是精确的，越往前误差越大。

因此作者提出，前半部分按第一个公式算，后半部分用第二个公式来算。如下：
在这里插入图片描述
后面有大佬提出了一种改进方法，把这两个公式线性插值一下，如下：

逐段的贝塞尔曲线

我们来看看贝塞尔曲线存在的一些问题，如下图一个十阶的贝塞尔曲线
在这里插入图片描述
这里会存在一些问题：

首先二阶三阶的贝塞尔曲线我们很容易通过控制点就想象出曲线的大致形状，但是当控制点很多的时候，我们很难想像出曲线大致的形状。
其次我们无法做到只修改曲线的某一部分，因为只要动一个控制点，整个曲线都会跟着发生变化，正所谓牵一发而动全身。
高阶的贝塞尔曲线往往会导致拐来拐去的，导致并不光滑，但是例如飞机汽车的表面都是光滑的。

因此在CAD中，通常不鼓励使用高阶的贝塞尔曲线，而是用低阶的贝塞尔曲线（通常是四阶的）相连，如下图：
在这里插入图片描述
如图，由三段四阶的贝塞尔曲线连成，分别是ABCD,HIJK和WXYZ三组控制点。

连续性

如上图：在K(或W）点有很明显的拐点，这种现象我们称之为不连续。在微积分的观念里，我们对曲线上的某一点左右进行求导，如果导数相等，那么就是连续的。对于这种连续性我们称为传统的连续性。那么对应K（或者W）点的左右导数根据贝塞尔曲线导数的性质我们很容易就得到3(K-J)和3(W-X)，很明显是不相等的，所以不连续。
但是这里我们可以举一个反例，如下图：
在这里插入图片描述
我们以AB两点作为控制点形成贝塞尔曲线，那么该曲线自然就是AB线段，同样的，我们以BC两点再形成一个贝塞尔曲线为BC线段，AB和BC形成一条直线段。此时B点左右的导数分别为(B-A)和(C-B)，这两个值明显是不相等的，按微积分里的定义就是不连续。可是它明明是一条直线，怎么会出现不连续呢？这个反例意味着传统的连续性不适用与CAD和图形学中形状的连续性，因此提出了几何连续性的概念。

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

那么我们就有如下定义，假设我们有P0,P1,P2…Pn,和Q0，Q1，Q2、、、Qm形成的两条贝塞尔曲线，如果：
在这里插入图片描述
连续性的效果如下：

贝塞尔曲面

我们接下来看看下面这种情况：
在这里插入图片描述
如图，我们有四条一模一样的贝塞尔曲线，图中的小红点是这些曲线在t从0到1变化时对应的点。
那么假设我们现在把这些小红点互相联系那，会形成啥呢？没错，会形成一个曲面，如下图：

出现了很多的摩尔纹。
但是这样形成的曲面其实有很大的问题，它仅仅只是在一个方向弯曲，在与其垂直的方向上仍旧还是笔直的，简单的修改下控制点的位置，就会出现失真的问题：
在这里插入图片描述
那么我们就要想办法来优化它，如果用着个四个小红点再做控制点，不就可以在垂直方向上又形成一条光滑的曲线了么？

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

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

贝塞尔曲面上的任意一点

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

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

要求贝塞尔曲面上一点
在这里插入图片描述
那么第i条n阶的贝塞尔曲线上的点P(u)即为：

那么我们就可以得到 m 个 P(u)，我们可以标记为P(u)0…P(u)m，这些 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阶连续（G0连续或者C0连续）**的表示如下：

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

每个控制点都一样，那么得到的曲线也都是一样的，所以曲线上对应v位置的点也都一样，所以也可以写成如下形式：
在这里插入图片描述
我们不管如何去一定其他控制点，链接处都保持一样的曲线，这就是所谓的0阶连续。
而对于1阶连续（G1或者C1），同样是对切向量进行限制。

在两个曲面的连接处上的任意一点，我们其实可以得到三个切线。而只要三个切线共面，我们就称之为1阶集合连续。如果除了小红点方向的切线，另外两个切线的值相当方向相反的话，那么就是一阶连续。

三角域贝塞尔曲面

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

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

如下图，如果我们只有如下三个顶点，那么每条边自然都是直线，形成的也就是三角形平面。
在这里插入图片描述
但是如果每条边变成多个顶点呢？如下图，每条边不就可以由这些顶点变成一条贝塞尔曲线了么。

在这里插入图片描述
三条边变曲线好搞，我们要怎么让整个三角形面都变成光滑的曲面呢？一样的思路，我们要求出三角形内任意一点的坐标，利用重心坐标，我们对三角形内任意点记作P(u，v,w)。
然后我们再定义一些规范，首先每条边的顶点数要相同，不能一条边4个顶点，另一条边10个顶点。这是在为难人。
如果每条边都有n+1个顶点，我们称之为n阶的三角域贝塞尔曲面。

为了方便理解，我们首先假设每条边被各自的顶点分成了等长的线段，那么我们就可以求出每个顶点的重心坐标了(至于哪个顶点对应u，哪个对应v，哪个对w，随你喜欢)，如下：
在这里插入图片描述
得到重心坐标后，我们把这些值同时乘以阶数n，得到的值就作为控制点的下标。例如上图中的三角形是3阶的，那么就把这些重心坐标的值都乘以3，例如重心坐标(1/3, 0, 2/3)的控制点，乘以3后为（1,0，2），那么该控制点就标记位P1.0.3，这样我们可以得到下图：
在这里插入图片描述
也就是说我们的控制点可以用Pi，j，k来表示，其中i+j+k=n。这其实就是我们三角域的伯恩斯坦多项式，因此对于三角域贝塞尔曲面上的任意一点P(u，v,w)我们可以得到这样一个公式：