贝塞尔曲线Bezier curve

现代机械设计方法，图形软件标准，样条曲线和曲面，图形变换

贝塞尔曲线

贝塞尔曲线简介

贝塞尔曲线(Bezier curve)，又称贝兹曲线或贝济埃曲线，是应用于二维图形应用程序的数学曲线。

特性：

• 使用$n+1$个（有序）控制点$\{P_0,P_1,P_2,...,P_n\}$来控制曲线的形状
• 曲线经过起点$P_1$和终点$P_n$，但不经过中间点$P_2-P_{n-1}$
• 贝塞尔曲线方程，方程的最高次数即是曲线的阶

$$Bez(t)=\sum _{i=0}^n{C}_n^i{P}_i(1-t{)}^{n-i}{t}^i=C_n^0{P}_0(1-t{)}^n{t}^0+C_n^1{P}_1(1-t{)}^{n-1}{t}^1+\dots \dots +C_n^{n-1}{P}_{n-1}(1-t{)}^1{t}^{n-1}+C_n^n{P}_n(1-t{)}^0{t}^n,t\in [0,1]$$

应用：

• photoshop中的钢笔工具就是应用的三次贝塞尔曲线

贝塞尔曲线推导

1. 二阶贝塞尔曲线的绘制

已知$3$个不共线的控制点$A,B,C$，线段$AB$和$BC$上各找到$D、E$ 共$2$个动点，线段$DE$上$1$个动点$F$，这$3$个动点满足 :

$$\frac{AD}{DB}=\frac{BE}{EC}=\frac{DF}{FE}$$

这就是抛物线的三切线定理，最终形成的二级贝塞尔曲线（抛物线）被直线$AB,BC,DE$相切，切点为$A,C,F$

$F$的集合$\{F\}$就是二阶贝塞尔曲线，方程可由一阶贝塞尔曲线：${}_{1}Bez(t)=(1-t)P_0+t{P}_1，t\in [0,1]$推导出：

• $D(t{)}^1=(1-t)A+tB$
• $E(t{)}^1=(1-t)B+tC$
• $F(t{)}^2=(1-t)D(t{)}^1+tE(T{)}^1=(1-t{)}^{2}A+2(t-1)B+t^{2}C$

上述方程使用点$P$代表其坐标的有序对$(x，y{)}^T$，似乎可以推广到三维空间

2. $n$阶贝塞尔曲线的绘制

• 已知$n+1$个不共线的控制点$\{P_0^0,P_1^0,P_2^0,...,P_n^0\}$
• 相邻点连接成$n$条线段$P_i^0{P}_{i+1}^0$，并在各个线段上找到$1$阶动点$P_i^1$ 共$n$个动点$\{P_0^1,P_1^1,P_2^1,...,P_{n-1}^1\}$
• 相邻点连接成$n-1$条线段$P_i^1{P}_{i+1}^1$，并在各个线段上找到$2$阶动点$P_i^2$ 共$n-1$个动点$\{P_0^2,P_1^2,P_2^2,...,P_{n-2}^2\}$
• ……
• 相邻点连接成的最后$1$条线段$P_0^{n-1}{P}_1^{n-1}$，并在线段上找到$n$阶动点$P_i^n$ 共$1$个动点$\{P_0^n\}$
• 这些动点$\{P_i^j\}$动点满足 :
  $$\frac{P_i^0{P}_i^1}{P_i^1{P}_{i+1}^0}=\frac{P_k^1{P}_k^2}{P_k^2{P}_{k+1}^1}=\dots \dots =\frac{P_0^{n-1}{P}_0^n}{P_0^n{P}_1^{n-1}}$$

• 动点集合$\{P_0^n\}$就是$n$阶贝塞尔曲线

可视化例子

一阶（两个控制点），即直线，曲线方程为一次多项式

三阶（四个控制点），曲线方程为三次多项式：

$Bez(t)=(1-t{)}^3{P}_0+3t(1-t{)}^2{P}_1+3t^2(1-t)P_2+t^3{P}_3,t\in [0,1]$

五阶（六个控制点），曲线方程为五次次多项式：

$$Bez(t)=(1-t{)}^5{P}_0+5t(1-t{)}^4{P}_1+10t^2(1-t{)}^3{P}_2+10t^3(1-t{)}^2{P}_3+5t^4(1-t)P_4+t^5{P}_5,t\in [0,1]$$

代码实现

N=length(control_points);
ta=zeros(N,N);%%对数组进行初始化
%%杨辉三角左右两边的值赋1

%%贝塞尔曲线方程的系数
% 杨辉三角的数的规律
% 1
% 1 1
% 1 2 1
% 1 3 3 1
% 1 4 6 4 1
for i=1:N
    ta(i,1)=1;
    ta(i,i)=1;
end
%%从第二个数开始，也就是从第三行开始，等于前列的左边加上正上方的一个
for row=2:N
    for col=2:row
        ta(row,col)=ta(row-1,col-1)+ta(row-1,col);
    end
end

%%曲线生成
for i=1:M
    t=i/M;%%确定每一个点的比例
    for k=0:N-1
        c=k;%分别确定a,b,c三个系数
        b=N-c-1;%分别确定a,b,c三个系数
        a=ta(N,k+1);%分别确定a,b,c三个系数
             
        p(i,1)=p(i,1)+a*(1-t)^b*t^c*control_points(k+1,1);%确定点的x坐标
       
        p(i,2)=p(i,2)+a*(1-t)^b*t^c*control_points(k+1,2);%确定点的y坐标
   end
  
end

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb

def inputPoints():
	controlPoints = []
	num = 1
	while True:
		print('\nenter %dst control point:'%num)
		x = input('x:')
		y = input('y:')
		#z = input('z:')
		print('Point:[%f,%f]'%(float(x),float(y)))
		i = input('Are you sure?(y or n)')
		if i=='y' or i=='Y':
			controlPoints.append([float(x),float(y)])
			inp = input('continue entering points or not?(y or n)')
			num = num + 1
			if inp == 'n':
				break
		else:
			continue

	#print(controlPoints)
	return controlPoints

def getInterpolationPoints(controlPoints, tList):
	n = len(controlPoints)-1
	interPoints = []
	for t in tList:
		Bt = np.zeros(2, np.float64)
		for i in range(len(controlPoints)):
			Bt = Bt + comb(n,i) * np.power(1-t,n-i) * np.power(t,i) * np.array(controlPoints[i])
		interPoints.append(list(Bt))

	return interPoints

if __name__ == '__main__':
#	points = inputPoints()
	points = [[1,1],[3,4],[5,5],[7,2]]
	tList = np.linspace(0,1,50)
	interPointsList = getInterpolationPoints(points, tList)
	x = np.array(interPointsList)[:,0]
	y = np.array(interPointsList)[:,1]

	plt.plot(x,y,color='b')
	plt.scatter(np.array(points)[:,0],np.array(points)[:,1],color='r')
	plt.show()

B样条曲线

贝塞尔曲线方程是对$n$个控制点$\{P_i\}$施加权重函数$W(t,n,i)$后加和得到的，即：$$B(t)=\sum _{i=0}^n{W}_{n,i}(t)P_i$$
其中，$W_{n,i}(t)=C_i^n(1-t{)}^{n-i}{t}^i$

可以发现，单独的一个控制点对整体曲线影响较大，曲线阶次过高（阶次与控制节点数量相关）

本质上就是在控制点前增加一个权重函数，然后加和；对权重函数进行修改，并使得控制点仅仅能影响局部形状，这就是B样条曲线(basic spline curve)的设计思路

B样条曲线的数学描述

对于$n+1$个控制点$P_0,P_1,\cdots ,P_n$，有一个包含$m+1$个节点的常数列表（或节点向量）$t_0,t_1,\cdots ,t_m$，其 k次$B$样条曲线表达式为（且$k=m-n-1>0$）

$$P(t)=\sum _{i=0}^n{W}_{k,i}(t)P_i，t\in [t_k,t_n)$$

其中，$W_{k,i}(t)$为k次$B$样条曲线的基函数（调和函数，权函数），并满足递推式：

$$k=0\mathrm{时}，B_{k,i}(t)=\{\begin{array}{l}1，t\in [t_i,t_{i+1})\\ 0，otherwise\end{array}$$

$$k>0\mathrm{时}，B_{k,i}(t)=\frac{t-t_i}{t_{i+k}-t_i}{B}_{k-1,i}(t)+\frac{t_{i+k+1}-t}{t_{i+k+1}-t_{i+1}}{B}_{k-1,i+1}(t)$$

显然，$W_{k,i}(t)$为一次函数，即高次B样条基函数为若干低次B样条基函数的线性组合。$N_{i,k}\left(t\right)$的次数k与控制节点的个数n无关，因此B样条曲线自由度更大

B样条曲线的次数指基函数多项式$B_{k,i}\left(t\right)$的最高次数

节点向量（kont vector）是一个一维单调非递减序列，元素$t_i$称为节点（kont），区间$[t_i,t_{i+1})$称为第 $i$ 个节点区间（knot range），节点在样条曲线上的映射$P(t_i)$称为曲节点（knot point）

在节点向量中，若某节点$t_i$出现$l$次，则称$t_i$是重复度为l的多重节点，否则为简单节点。与贝塞尔曲线不同，仅当B样条曲线首末节点重复度为$k+1$时，曲线本身才穿过首末控制点。

基函数$B_{k,i}\left(t\right)$在区间$[t_i,t_{i+k+1}]$上非零，因为该区间上总存在不为零的零阶基函数$B_{0,i}$该区间称为支撑区间，对应样条曲线上的区段称为支撑曲线。由于$B_{k,i}\left(t\right)$直接与控制节点${\mathit{p}}_i$相乘，所以${\mathit{p}}_i$只影响其支撑区间$[t_i,t_{i+k+1}]$上对应支撑曲线的形状。

所以B样条曲线可视为若干段贝塞尔曲线的拼接，是贝塞尔曲线的推广，相邻贝塞尔曲线间存在若干重合节点，保留了对称性、几何不变性、变差伸缩性等优良特性。

B样条曲线分类

节点向量（kont vector）是一个一维单调非递减序列

根据节点分布的性质，可以将样条曲线进行分类

• 均匀 B 样条：节点均匀分布，所有节点区间等长
• 准均匀 B 样条：在开始和结束处的节点可重复，中间节点均匀分布；节点向量中的首末节点重复度为$k+1$，其余节点沿数轴方向等距均匀分布且重复度为1。当$k=n$时，B样条基函数$B_{k,i}\left(t\right)$退化为伯恩斯坦多项式，即B样条曲线退化为贝塞尔曲线。
• 非均匀 B 样条：节点非均匀分布，可任意分布
• 分段贝塞尔曲线
• NURBS曲线：B样条无法描述圆锥曲线，为解决此问题，产生了非均匀有理B样条（non-uniform rational b-spline, NURBS）

(准/非)均匀B样条曲线

对于$n+1$个控制点$P_0,P_1,\cdots ,P_n$，有一个包含$m+1$个节点的常数列表（或节点向量）$t_0,t_1,\cdots ,t_m$

其 k次均匀$B$样条曲线节点表达式为（且$k=m-n-1>0$）
例如，令$n=4，m=9，k=4$，则$knots=[0，1/9，2/9，3/9，4/9，5/9，6/9，7/9，8/9，1]=\{i/m\}$ ，其中$i=0,1,3,\dots \dots ,m$

其 k次准均匀$B$样条曲线节点表达式为（且$k=m-n-1>0$）
例如，令$n=4，m=7，k=2$，则$knots=[0，0，1/5，2/5，3/5，4/5，1，1]$

其 k次非均匀$B$样条曲线节点表达式为（且$k=m-n-1>0$）
例如，令$n=4，m=7，k=2$，则可令$knots=[0，1/9，2.5/9，4/9，5.6/9，7/9，8/9，1]$ 或者$knots=[0，0，2/9，4/9，5/9，1，1，1]$

代码实现

import numpy as np
import matplotlib.pyplot as plt

# 计算在某一特定t下的 B_{i,k}
def getBt(controlPoints, knots, t):
	# calculate m,n,k
	m = knots.shape[0]-1
	n = controlPoints.shape[0]-1
	k = m - n - 1
	# initialize B by zeros 
	B = np.zeros((k+1, m))

	# get t region
	tStart = 0
	for x in range(m+1):
		if t==1:
			tStart = m-1
		if knots[x] > t:
			tStart = x-1
			break
	 
	# calculate B(t)
	for _k in range(k+1):
		if _k == 0:
			B[_k, tStart] = 1
		else:
			for i in range(m-_k):
				if knots[i+_k]-knots[i]== 0:
					w1 = 0
				else:
					w1 = (t-knots[i])/(knots[i+_k]-knots[i]) 
				if knots[i+_k+1]-knots[i+1] == 0:
					w2 = 0
				else:
					w2 = (knots[i+_k+1]-t)/(knots[i+_k+1]-knots[i+1])
				B[_k,i] = w1*B[_k-1, i] + w2*B[_k-1, i+1]
	return B

# 绘制 B_{i,k}(t)函数
def plotBt(Bt,num, i,k):
	print(k,i)
	Bt = np.array(Bt)
	tt = np.linspace(0,1,num)
	yy = [Bt[t,k,i] for t in range(num)]
	plt.plot(tt, yy)

# 根据最后一列（最高阶次）的 B(t)，即权重，乘以控制点坐标，从而求出曲线上点坐标
def getPt(Bt, controlPoints):
	Bt = np.array(Bt)
	ptArray = Bt.reshape(-1,1) * controlPoints
	pt = ptArray.sum(axis = 0)
	return pt

# 绘制出生成的样条曲线: useReg 表示是否使用曲线有效定义域[t_k, t_{m-k}]
def main1(useReg = False):
	controlPoints = np.array([[50,50], [100,300], [300,100], [380,200], [400,600]])
	knots = np.array([0,1/9,2/9,3/9,4/9,5/9,6/9,7/9,8/9,1])
	m = knots.shape[0]-1
	n = controlPoints.shape[0]-1
	k = m - n - 1
	print('n:',n)
	print('m:',m)
	print('k:',k)
    
	for t in np.linspace(0,1,100):
		if useReg and not(t >= knots[k] and t<= knots[n+1]):
			continue
		Bt = getBt(controlPoints, knots, t)
		Pt = getPt(Bt[k, :n+1], controlPoints)
		plt.scatter(Pt[0],Pt[1],color='b')
	plt.scatter(controlPoints[:,0], controlPoints[:,1],color = 'r')
	plt.show()

# 绘制 B_{i,k} 变化图:如果不给定{i,k}则显示所有B{i,k}(t)图像
def main2(i=-1,k=-1):
	controlPoints = np.array([[50,50], [100,300], [300,100], [380,200], [400,600]])
	knots = np.array([0,1/9,2/9,3/9,4/9,5/9,6/9,7/9,8/9,1])
	m = knots.shape[0]-1
	n = controlPoints.shape[0]-1
	k = m - n - 1
	print('n:',n)
	print('m:',m)
	print('k:',k)
	B = []
	num = 100 # 离散点数目
	for t in np.linspace(0,1,num):
		Bt = getBt(controlPoints, knots, t)
		B.append(list(Bt))

	figure1 = plt.figure('B_{i,k}')
	if i==-1:
		fig = []
		for i in range(n+1):
			for k in range(k+1):
				plotBt(B,num, i,k)
				fig.append('B_{%d,%d}'%(i,k))
	else:
		plotBt(B,num, i,k)
		fig.append('B_{%d,%d}'%(i,k))
	plt.legend(fig)
	plt.show()   
    
if __name__ == '__main__':
    main1()
    main2()

Nurbs曲线

ISO规定，PHIGS Plus的扩充部分，Bezier、有理Bezier、均匀B样条和非均匀B样条都被统一到NURBS 中。

B样条曲面、及其特例的Bezier曲面都不能精确表示除抛物面以外的二次曲面，而只能给出近似表示

在曲线曲面描述中，B样条方法更多地以非均匀类型出现，而均匀、准均匀、分段Bezier三种类型又被看成是非均匀类型的特例，所以习惯上称之为非均匀有理B样条（Non-Uniform Rational B-Splines）方法，简称为NURBS方法

NURBS方法提出的主要理由是，寻找与描述自由曲线曲面的B样条方法相统一的，而又能精确表示二次曲线曲面的数学方法

• 非均匀B样条采用分段参数整数多项式，而NURBS方法采用分子分母分别是分段参数多项式函数与分段多项式的分式表示，是有理的
• 与有理Bezier方法一样，NURBS方法引入了权因子和分母NURBS方法是在有理Bezier方法与非有理B样条方法的基础上发展起来的

NURBS的优点主要表现在以下几个方面：

• 将初等曲线曲面与自由曲线曲面的表达方式统一起来
• 增加了权因子，有利于对曲线曲面形状的控制和修改
• 非有理B样条、有理与非有理Bezier方法是NURBS的特例
• 在比例、旋转、平移、错切以及平行和透视投影变换下是不变的

NURBS曲线有三种表示方法：

• 分式表示是有理的由来，它说明：NURBS曲线是非有理与有理Bézier和非有理B样条曲线的推广：但却难以从中了解更多的性质。
• 在有理基函数表示形式中，可从有理基函数的性质清楚地了解NURBS曲线的性质。
• 齐次坐标表示形式说明：NURBS曲线是它的控制顶点的齐次坐标或带权控制点在高一维空间里所定义的非有理B样条曲线在$\omega =1$超平面上的投影。
• 不仅包含了明确的几何意义，而且也说明：非有理B样条曲线的大多数算法都可以推广应用于NURBS曲线。

核心是其精心设计的数据结构，包括控制点网格、权重数组、 knot向量等，利用NumPy进行数值计算，NURBS-Python库具有高度的兼容性和性能，NURBS-Python支持多种格式的导入和导出，如IGES和STEP

矢量函数

一条$p$次NURBS曲线可以表示为有理矢量函数形式：$C(u)=\frac{\sum _{i=0}^n{\omega }_i{N}_{i,p}(u)P_i}{\sum _{i=0}^n{\omega }_i{N}_{i,p}(u)}$

式中，$P_i$是曲线的控制点，${\omega }_i$是权因子，$N_{i,p}(u)$是定义在非周期节点矢量$U$上的$p$次B样条基函数，$N_{i,p}(u)$定义为：

$N_{i,0}(u)=\{\begin{array}{c}1;u_i\le u\le u_{i+1}\\ 0;ohterwise\end{array}$

$N_{i,p}(u)=\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(u)$

节点矢量$U$为：$U=\{a_0,a_1,\dots ,a_p,u_{p+1},\dots ,u_{m-p-1},b_{m-p},\dots ,b_m\}$，一般情况取$a=0,b=1$，其值单调不减，且${\omega }_i>0$

为求简化，可令：

$R_{i,p}(u)=\frac{{\omega }_i{N}_{i,p}(u)}{\sum _{j=0}^n{\omega }_j{N}_{j,p}(u)}$，称$\{R_{i,p}(u)\}$为有理基函数

则$C(u)=\sum _{i=0}^n{P}_i{R}_{i,p}(u),0\le u\le 1$，称$u$为分段有理函数

示例

7个控制点$\{P_0,P_1,P_2,P_3,P_4,P_5,P_6\}$，7个权重系数$\{{\omega }_0,{\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4,{\omega }_5,{\omega }_6\}=\{1,1,1,3,1,1,1\}$，3次NURBS曲线，节点矢量$U=\{0,0,0,0,\frac{1}{4},\frac{1}{2},\frac{3}{4},1,1,1,1\}$

则可计算出$N_{i,p}(u)$或者$R_{i,p}(u)$

Frenet框架

曲线的曲率和挠率，切矢量、法矢量、次法矢量