贝塞尔曲线的python实现

实现代码在最后。

一、参数方程

• 直角坐标系下参数方程
• 极坐标系下参数方程

直角坐标系下的参数方程：
$x=f(t) \quad y=g(t)$x=f(t)y=g(t)

极坐标系的参数则是

$ρ=f(t) \quad θ=g(t)。$ρ=f(t)θ=g(t)。

二、二项式定理

贝塞尔其实就是比例点在空间上走过的轨迹。在介绍贝塞尔曲线之前，先复习一下高中数学——二项式定理：
$(a+b)^n=c^0_na^nb+c^1_na^{n-1}b^2+\cdots+c^i_na^{n-i}b^i+\cdots+c^n_nb^{n}=\sum_{i=0}^nc^i_na^{n-i}b^i$(a+b)n=cn0​anb+cn1​an−1b2+⋯+cni​an−ibi+⋯+cnn​bn=i=0∑n​cni​an−ibi
二项式展开式一共有$n+1$n+1项，第$i$i项为$c^i_na^{n-i}b^i$cni​an−ibi

矩阵形式：

二项展开式性质

• 对称性 $C_n^k=C_n^{n-k}$Cnk​=Cnn−k​
• 先增后减，峰值在中间，若幂为偶，在$\frac {n}{2}+1$2n​+1项最大；若幂为奇，中间两项都为最大
• 系数和为$2^n$2n
• 奇数和偶数项的系数和相等

三、贝塞尔曲线坐标分量是关于t的参数方程

3.1 一阶贝塞尔

两个控制点，1阶贝塞尔。

设$n$n维向量为$\vec {P0}$P0 $\vec {P1}$P1，那么一阶贝塞尔公式可以表示成：
$Pi(t)=(1-t)P0+tP1 \quad t\in[0,1]$Pi(t)=(1−t)P0+tP1t∈[0,1]

3.2 二阶贝塞尔

三个控制点，2阶贝塞尔。
$Pi(t)=(1-t)^2P0+2t(1-t)P1+t^2P2 \quad t\in[0,1]$Pi(t)=(1−t)2P0+2t(1−t)P1+t2P2t∈[0,1]

3.3 n阶贝塞尔

n+1个控制点，n阶
$P(t)=\sum_{i=0}^nP_iB_{i,n}(t) \quad t\in[0,1] \quad i=0,1,2, \cdots n$P(t)=i=0∑n​Pi​Bi,n​(t)t∈[0,1]i=0,1,2,⋯n

$B_{i,n}(t)$Bi,n​(t)中，$t$t的不同取值影响了各个控制点的权重。改变任意一点都将会影响最终的插补效果。

3.4 贝塞尔曲线的性质

• 各项系数和为1
• 对称性

前两个都是二项式的性质

• 递归性
  $B_{i,n}(t)=(1-t)B_{i,n-1}(t)+tB_{i-1,n-1}(t)$Bi,n​(t)=(1−t)Bi,n−1​(t)+tBi−1,n−1​(t)

• 凸包性

所有的控制点都在最小凸多边形中，不是按照顺序围城的最小多边形。

• 导数性质

可用于两个贝塞尔曲线的连接。

四、python的实现

import matplotlib.pyplot as plt
import numpy as np
import math

def line_interpolation(startp,endp,num):
    dx=(endp[0]-startp[0])/num
    dy=(endp[1]-startp[1])/num
    dz=(endp[2]-startp[2])/num

    x=np.zeros((num+1,1))
    y=np.zeros((num+1,1))
    z=np.zeros((num+1,1))
    
    for i in range(0,num+1):
        x[i]=startp[0]+dx*i
        y[i]=startp[1]+dy*i
        z[i]=startp[2]+dz*i

    return np.hstack((x,y,z))


# line up all points
def line_up(points,num):
    result=line_interpolation(points[0],points[1],num)
    for i in range(1,points.shape[0]-1):
        pi=line_interpolation(points[i],points[i+1],num)
        result=np.hstack((result,pi))
        
    return result.T


def getbin(i,n,t):
     bt=(math.factorial(n)/(math.factorial(i)*math.factorial(n-i)))*(t**i)*(1-t)**(n-i)
     return bt

points=np.array([
    [1,3,0],
    [1.5,1,0],
    [4,2,0],
    [4,3,4],
    [2,3,11],
    [5,5,9]
    ])

pis=line_up(points,1000);

fig=plt.figure()
ax = fig.gca(projection='3d')

for i in range(0,pis.shape[0]//3):
    ax.plot3D(pis[0+3*i],pis[1+3*i],pis[2+3*i],color='k')

for i in range(0,points.shape[0]):
    ax.scatter(points[i][0],points[i][1],points[i][2],marker='o',color='r')
    ax.text(points[i][0],points[i][1],points[i][2],i,size=12)


# Bezier interpolation
bts=np.zeros(points.shape[0])
matbts=np.zeros((points.shape[0],points.shape[1]))

t_step=0.01
pilst =[]
matpi=np.zeros((int(1/t_step+1),3));

for t in np.arange(0,1+0.01,0.01):
    for i in range(0,points.shape[0]):
        matbts[i]=getbin(i,points.shape[0]-1,t)*points[i]
    pi=sum(matbts).tolist()
    pilst.append(pi)
matpi=np.array(pilst)


ax.plot3D(matpi[:,0],matpi[:,1],matpi[:,2],color='r')

plt.show()

[1] 曲线篇: 贝塞尔曲线