3次B样条曲线反求控制顶点
1. 说明
- 代码不完善，因为方程组的边界条件不太清楚，总是有问题。。。。
2. 边界条件取二次样条在端点处切矢
- 该边界条件指的是新增方程组条件，端点处的导数为二次样条在端点处的导数。
# 基函数
def N(i,k,T,t): # 曲线N_ik
if k == 0:
if tT[i+1]:
return 0
else:
return 1
else:
result = (t-T[i])/(T[i+k]-T[i])*N(i,k-1,T,t) + (T[i+k+1]-t)/(T[i+k+1]-T[i+1])*N(i+1,k-1,T,t)
return result
def Inter_func(P,judge): # 抛物线插值函数 传入参数为(坐标点,判断条件) 判断条件：0表示前三个，1表示后三个
x = P.tolist()[0]
y = P.tolist()[1]
if judge == 0:
value = 2*x[0]*((x[1]-x[2])*y[0]+(x[2]-x[0])*y[1]+(x[0]-x[1])*y[2])-np.power(x[1],2)*(y[0]-y[2])-np.power(x[2],2)*(y[1]-y[0])-np.power(x[0],2)*(y[2]-y[1])
value = value/((x[0]-x[1])*(x[0]-x[2])*(x[1]-x[2]))
else:
value = 2*x[2]*((x[1]-x[2])*y[0]+(x[2]-x[0])*y[1]+(x[0]-x[1])*y[2])-np.power(x[1],2)*(y[0]-y[2])-np.power(x[2],2)*(y[1]-y[0])-np.power(x[0],2)*(y[2]-y[1])
value = value/((x[0]-x[1])*(x[0]-x[2])*(x[1]-x[2]))
return value
# 曲线绘制
def B_curve_inverse(k,n,P_,zero_num):
P = P_
# 取前三个插值点
P_first = np.array(P)[:,:3]
# 前三个插值点传入抛物线插值函数，返回一个在第一个点处的导数
p0_diff = Inter_func(P_first,0) # P0点的导数值
p0 = (np.array([P_first[0][0],P_first[1][0]])-np.array([P_first[0][0],p0_diff])/3).tolist()
P = np.insert(P,0,p0,axis=1) #将该点的导数插入
# 取后三个插值点
P_last = np.array(P)[:,-3:]
# 后三个插值点传入抛物线插值函数，返回一个在最后一个点处的导数
pn_diff = Inter_func(P_last,1) # Pn点的导数值
pn = (np.array([P_last[0][2],P_last[1][2]])+np.array([P_last[0][2],pn_diff])/3).tolist()
P = np.insert(P,len(P[0]),pn,axis=1) #将该点的导数插入
# 获取系数矩阵
coef = []
for i in range(len(P[0])):
if i == 0:
col = [-1,1]+np.repeat(0,zero_num+1).tolist()# 将一行系数的0进行填充
coef.append(col)
elif i == len(P[0])-1:
col = np.repeat(0,zero_num+1).tolist()+[-1,1]# 将一行系数的0进行填充
coef.append(col)
else:
# print(zero_num-i+1)
col = np.repeat(0,i-1).tolist()+[1/6,4/6,1/6]+np.repeat(0,zero_num-i+1).tolist()
coef.append(col)
coef1 = np.matrix(coef).I # 系数矩阵取逆
V = coef1*np.matrix(P).T # 获得控制顶点
V = V.T.tolist() # 将控制顶点矩阵变为2*(n+2)的list
# 控制顶点的x=V[0]，控制顶点的y=V[1]
# 绘制曲线
T = np.linspace(0,1,k+len(V[0])+1) # T存储节点
t_x = np.linspace(0,1,100) # t_x存储每一个t值
X = V[0]
Y = V[1]
x = [] # 用来存储曲线的x值
y = [] # 用来存储曲线的y值
for i in range(len(V[0])-k): # for循环用作获取第几条曲线段的数值
result = pd.DataFrame(t_x,columns=['t'])
for j in range(k+1):
result1 = []
for t in t_x:
result1.append(N(i+j,k,T,t))
result['N_{0}{1}'.format(i+j,k)] = result1 # 将N_ij存入dataframe
# 把Dataframe中 T[i+k]<=t<=T[i+k+1] 的数据取出来 保存为 matrix, 然后用matrix*np.matrix(X[i:i+k+1]).T获取曲线的x值
# lambda x: x>T[i+k] and x= 或 <= 曲线有时会有问题 比如 V = [[0,0,2,3,2,1],[1,3,3,2,1,1]]
Ni_matrix = np.matrix(result[result['t'].apply(lambda x: x>=T[i+k] and x<=T[i+k+1])].iloc[:,1:]) # Ni_matrix 是一个 t_ba*j 维矩阵
x = x+((Ni_matrix*np.matrix(X[i:i+k+1]).T).T).tolist()[0]
y = y+((Ni_matrix*np.matrix(Y[i:i+k+1]).T).T).tolist()[0]
fig = plt.figure()
# plt.plot(V[0],V[1],marker='o',markerfacecolor='white')
# print(V)
plt.plot(P_[0],P_[1],'o')
plt.plot(x,y)
plt.show()
return coef
# 参数传入
# P = np.array([[0,1,10,12,15],[0,10,8,6,1]]).astype(float).tolist() # 插值点
P = np.array([[1,2,3,4,4,5,7,8],[1,2,5,2,1.5,1,0.5,-1.5]]).astype(float).tolist() # 插值点
# P = np.array([[1,2,3,4,5],[0,10,8,6,1]]).astype(float).tolist() # 插值点
n =len(P[0]) # 插值点个数
k = 3 # 样条曲线次数
zero_num = n+2-k # 矩阵需要补0的个数
# print(zero_num)
coef = B_curve_inverse(k,n,P,zero_num)
3. 边界条件取0
- 该边界条件指的是新增方程组条件，使得
P
0
=
P
1
,
P
n
−
1
=
P
n
P_0=P_1,P_{n-1}=P_n
P0=P1,Pn−1=Pn。
# 基函数
def N(i,k,T,t): # 曲线N_ik
if k == 0:
if tT[i+1]:
return 0
else:
return 1
else:
result = (t-T[i])/(T[i+k]-T[i])*N(i,k-1,T,t) + (T[i+k+1]-t)/(T[i+k+1]-T[i+1])*N(i+1,k-1,T,t)
return result
def Inter_func(P,judge): # 抛物线插值函数 传入参数为(坐标点,判断条件) 判断条件：0表示前三个，1表示后三个
x = P.tolist()[0]
y = P.tolist()[1]
if judge == 0:
value = 2*x[0]*((x[1]-x[2])*y[0]+(x[2]-x[0])*y[1]+(x[0]-x[1])*y[2])-np.power(x[1],2)*(y[0]-y[2])-np.power(x[2],2)*(y[1]-y[0])-np.power(x[0],2)*(y[2]-y[1])
value = value/((x[0]-x[1])*(x[0]-x[2])*(x[1]-x[2]))
else:
value = 2*x[2]*((x[1]-x[2])*y[0]+(x[2]-x[0])*y[1]+(x[0]-x[1])*y[2])-np.power(x[1],2)*(y[0]-y[2])-np.power(x[2],2)*(y[1]-y[0])-np.power(x[0],2)*(y[2]-y[1])
value = value/((x[0]-x[1])*(x[0]-x[2])*(x[1]-x[2]))
return value
# 曲线绘制
def B_curve_inverse(k,n,P_,zero_num):
P = P_
# 取前三个插值点
P_first = np.array(P)[:,:3]
# 前三个插值点传入抛物线插值函数，返回一个在第一个点处的导数
p0_diff = Inter_func(P_first,0) # P0点的导数值
p0 = (np.array([P_first[0][0],P_first[1][0]])-np.array([P_first[0][0],p0_diff])/3).tolist()
# P = np.insert(P,0,p0,axis=1) #将该点的导数插入
P = np.insert(P,0,[0,0],axis=1)
# 取后三个插值点
P_last = np.array(P)[:,-3:]
# 后三个插值点传入抛物线插值函数，返回一个在最后一个点处的导数
pn_diff = Inter_func(P_last,1) # Pn点的导数值
pn = (np.array([P_last[0][2],P_last[1][2]])+np.array([P_last[0][2],pn_diff])/3).tolist()
# P = np.insert(P,len(P[0]),pn,axis=1) #将该点的导数插入
P = np.insert(P,len(P[0]),[0,0],axis=1)
# 获取系数矩阵
coef = []
for i in range(len(P[0])):
if i == 0:
col = [-1,1]+np.repeat(0,zero_num+1).tolist()# 将一行系数的0进行填充
coef.append(col)
elif i == len(P[0])-1:
col = np.repeat(0,zero_num+1).tolist()+[-1,1]# 将一行系数的0进行填充
coef.append(col)
else:
# print(zero_num-i+1)
col = np.repeat(0,i-1).tolist()+[1/6,4/6,1/6]+np.repeat(0,zero_num-i+1).tolist()
coef.append(col)
coef1 = np.matrix(coef).I # 系数矩阵取逆
V = coef1*np.matrix(P).T # 获得控制顶点
V = V.T.tolist() # 将控制顶点矩阵变为2*(n+2)的list
# 控制顶点的x=V[0]，控制顶点的y=V[1]
# V = [[1,1.5,3,3.8,3.2,5,8],[0,2.6,3,2.2,-0.3,-0.8,2]]
# 绘制曲线
T = np.linspace(0,1,k+len(V[0])+1) # T存储节点
t_x = np.linspace(0,1,100) # t_x存储每一个t值
X = V[0]
Y = V[1]
x = [] # 用来存储曲线的x值
y = [] # 用来存储曲线的y值
for i in range(len(V[0])-k): # for循环用作获取第几条曲线段的数值
result = pd.DataFrame(t_x,columns=['t'])
for j in range(k+1):
result1 = []
for t in t_x:
result1.append(N(i+j,k,T,t))
result['N_{0}{1}'.format(i+j,k)] = result1 # 将N_ij存入dataframe
# 把Dataframe中 T[i+k]<=t<=T[i+k+1] 的数据取出来 保存为 matrix, 然后用matrix*np.matrix(X[i:i+k+1]).T获取曲线的x值
# lambda x: x>T[i+k] and x= 或 <= 曲线有时会有问题 比如 V = [[0,0,2,3,2,1],[1,3,3,2,1,1]]
Ni_matrix = np.matrix(result[result['t'].apply(lambda x: x>=T[i+k] and x<=T[i+k+1])].iloc[:,1:]) # Ni_matrix 是一个 t_ba*j 维矩阵
x = x+((Ni_matrix*np.matrix(X[i:i+k+1]).T).T).tolist()[0]
y = y+((Ni_matrix*np.matrix(Y[i:i+k+1]).T).T).tolist()[0]
fig = plt.figure()
# plt.plot(V[0],V[1],marker='o',markerfacecolor='white')
# # print(V)
plt.plot(P_[0],P_[1],'o')
plt.plot(x,y)
plt.show()
return coef,V,np.matrix(P).T
# 参数传入
# P = np.array([[0,1,2,3,4],[0,2,6,4.5,3]]).astype(float).tolist() # 插值点
P = np.array([[1,2,3,4,4.5,5,7,8],[1,2,5,2,1.3,1,0.5,-1.5]]).astype(float).tolist() # 插值点
# P = np.array([[1,2,3,4,5],[0,10,8,6,1]]).astype(float).tolist() # 插值点
n =len(P[0]) # 插值点个数
k = 3 # 样条曲线次数
zero_num = n+2-k # 矩阵需要补0的个数
# print(zero_num)
coef,V,p = B_curve_inverse(k,n,P,zero_num)
