de Castelijau算法的python实现

Beizer曲线的阶数

在不考虑分段的情况下，给定n个初始控制点，Beizer曲线的阶数是n-1

de Castelijau算法的关键

关键在于初始控制点和下一级控制点的递推公式

$P_i^k(t)=(1-t)P_i^{k-1}(t)+t{P}_{i+1}^{k-1}(t)$

注意$i$和$k$的取值范围

python代码

代码输出给定点的Beizer曲线的参数方程，并画出每一级控制多边形。

第一次用python写科学计算，欢迎交流。

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

###########确定Beizer曲线的次数，de Castelijau算法的参数，n+1个点i=0-n
t = eval(input("请输入平面Beizer曲线de castelijau算法的参数t:"))
###########给出初始控制点，并用数组储存
def cmul(*b):
    count = []
    for i in [*b]:
        count.append(i)
    return count
p = eval("cmul({})".format(input("请按照格式[x,y]输入初始控制顶点P并用逗号隔开：")))
p0 = np.array(p).T
p0 = p0.tolist()
#############给出顶点所确定的曲线的次数
n = len(p)-1
############都是列表类型的
x = p0[0]
y = p0[1]


###########计算阶乘
def jc(n):
    if n ==1 or n ==0:
        return 1
    else:
        return jc(n-1)*n
###########Berstein基函数列
u =sym.symbols('u') #定义符号变量x，即算法的t
def Bernstein(n,u):
    list = []
    for i in range(n+1):
        list.append(jc(n)*(u**i)*((1-u)**(n-i))/(jc(i)*jc(n-i)))
    return list
###########输出参数方程
sumx = 0
sumy = 0

for i in range(len(x)):
    sumx=sumx+(x[i]*Bernstein(n,u)[i])
    sumy=sumy+(y[i]*Bernstein(n,u)[i])

s=(sumx.evalf(subs = {u:0.4}),sumy.evalf(subs = {u:0.4}))
sum = [sumx,sumy]
print("曲线的参数方程为:p(u)={}".format(sum))
print("参数方程求得曲线上一点 p({})为{}".format(t,s))


################下一级控制顶点坐标函数
def nextpoint(p):
    d = len(p)#当前级控制顶点个数
    q = []#把下一级的储存在空列表中
    for i in range(d-1):#下一次控制顶点要减一
        q.append((1-t)*p[i] + t*p[i+1])
    return q

#输出每一级控制顶点坐标
a=[]
b=[]
a.append(x)
b.append(y)
while len(x)>1:
    x = nextpoint(x)
    y = nextpoint(y)
    a.append(x)
    b.append(y)
print("每一级控制顶点的横坐标依次为:{}".format(a))
print("每一级控制顶点的纵坐标依次为:{}".format(b))
print("de Castelijau求得曲线上一点 p({})为({},{})".format(t,a[n][0],b[n][0]))

################绘制曲线和控制多边形
def plot():
    c = []
    z = []
    for j in range(0,10000,1):
        j=j/10000
        c.append(sumx.evalf(subs = {u:j}))
        z.append(sumy.evalf(subs = {u:j}))
    plt.plot(c,z,'.',markersize=0.5)
    
    for k in range(n+1):
        x= np.array(a[k])
        y =np.array(b[k])
        plt.plot(x,y)
        plt.plot(x, y, '*')
plot()

结果