【Python】四点插值细分法

知识来源：Nira DYN, David LEVJN：A 4-point interpolatory subdivision scheme for curve design 论文

• 最左边的 p − 1 = p n , p − 2 = p n − 1 p_{-1} = p_n,p_{-2} = p_{n-1} p−1​=pn​,p−2​=pn−1​可以在代码中使用列表的倒序下标 p − 1 p_{-1} p−1​ = p[-1]， p − 2 p_{-2} p−2​=p[-2]；
• 最右边的 p n + 1 = p 0 , p n + 2 = p 1 p_{n+1} = p_0,p_{n+2}=p_1 pn+1​=p0​,pn+2​=p1​可以在代码中对超出n的记为m-n；
• 故每个点列都需要增加3个末端点；
• 注意再添加控制点的时候需要输入float浮点数数据类型，int整数数据类型可能会产生计算错误，详情可见最后一张图。

# 导入相关包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pylab import xticks,yticks
import math
import random

plt.rcParams["font.sans-serif"]=["SimHei"] #设置字体
plt.rcParams["axes.unicode_minus"]=False #该语句解决图像中的“-”负号的乱码问题

x = [-0.5,1,3,5,4,5,3,1,-0.5,0,-0.5]
y = [2,1,2.4,1,0,-1,-2.4,-1,-2,0,2]

def point_cp(points,w):  # points 包含了连续的四个点，且下标取i时，连续的四个点为i-1, i, i+1, i+2
    p = (1/2+w)*(points[1]+points[2])-w*(points[0]+points[3])
    return p

def add_points(P,w):
    n = len(P)
    for k in range(10):
        n = len(P)
        p_new = []  # 存储新增节点
        for i in range(n):
            if i==n-2:
                points = []
                for j in range(3):
                    points.append(P[i+j-1])
                points.append(P[0])
                p = point_cp(points,w)
                p_new.append(p)
            elif i==n-1:
                points = []
                for j in range(2):
                    points.append(P[i+j-1])
                points.append(P[0])
                points.append(P[1])
                p = point_cp(points,w)
                p_new.append(p)
            else:
                points = []
                for j in range(4):
                    points.append(P[i+j-1])
                p = point_cp(points,w)
                p_new.append(p)
        for m in range(n,0,-1):
            P = np.insert(P,m,p_new[m-1])
        
    return P

x = [-0.5,1,3,5,4,5,3,1,-0.5,0]
y = [2,1,2.4,1,0,-1,-2.4,-1,-2,0]
P_x = [-0.5,1,3,5,4,5,3,1,-0.5,0]
P_y = [2,1,2.4,1,0,-1,-2.4,-1,-2,0]
Px = add_points(P_x,1/16)
Py = add_points(P_y,1/16)

x = [0,3.0,6.0,6.0,4.0,4.0,3.0,2.0,2.0,0.0,0]
y = [1,1.4,1,0,0,0.3,0.5,0.3,0,0,1]
P_x = [0,3.0,6.0,6.0,4.0,4.0,3.0,2.0,2.0,0.0]
P_y = [1,1.4,1,0,0,0.3,0.5,0.3,0,0]

# x = [-0.5,1,3,5,4,5,3,1,-0.5,0]
# y = [2,1,2.4,1,0,-1,-2.4,-1,-2,0]
fig = plt.figure()
plt.plot(x,y,'-.')
plt.plot(Px,Py)

x = [0,3.0,6.0,6.0,4.0,4.0,3.0,2.0,2.0,0.0,0]
y = [1,1.4,1,0,0,0.3,0.5,0.3,0,0,1]
# P_x = [0,3.0,6.0,6.0,4.0,4.0,3.0,2.0,2.0,0.0]
P_x = [0,3,6,6,4,4,3,2,2,0]
P_y = [1,1.4,1,0,0,0.3,0.5,0.3,0,0]

Px = add_points(P_x,1/15)
Py = add_points(P_y,1/15)
fig = plt.figure()
plt.plot(x,y,'-.')
plt.plot(Px,Py)
plt.show()