Reputation: 21
I started getting into spline principles to use them in my studies. However, I got stuck in programming the code for Wikipedia [example] . The final graph doesn't correspond to the Wiki's example. I get two curves that are not connected, despite the introduction of the parameter t(x):
Would you mind helping me in finding the cause of the error? Thanks.
import numpy as np
import matplotlib.pyplot as plt
x0, y0, x1, y1, x2, y2 = -1, 0.5, 0, 0, 3, 3
def spline(x0, y0, x1, y1, x2, y2):
a11 = 2 / (x1 - x0)
a12 = 1 / (x1 - x0)
a13 = 0
a21 = a12
a22 = 2 * (1 / (x1 - x0) + 1 / (x2 - x1))
a23 = 1 / (x2 - x1)
a31 = a13
a32 = a23
a33 = 2 / (x2 - x1)
b1 = 3 * ((y1 - y0) / (x1 - x0) ** 2)
b2 = 3 * ((y1 - y0) / (x1 - x0) ** 2 + (y2 - y1) / (x2 - x1) ** 2)
b3 = 3 * ((y2 - y1) / (x2 - x1) ** 2)
A = np.array([[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]])
L = np.array([[-b1], [-b2], [-b3]])
AT = np.transpose(A)
ATA = np.matmul(AT, A)
invATA = np.linalg.inv(ATA)
ATL = np.matmul(AT, L)
X = - np.matmul(invATA, ATL)
[k0, k1, k2] = X
a1 = k0 * (x1 - x0) - (y1 - y0)
b1 = -k1 * (x1 - x0) + (y1 - y0)
a2 = k1 * (x2 - x1) - (y2 - y1)
b2 = -k2 * (x2 - x1) + (y2 - y1)
return a1[0], b1[0], a2[0], b2[0]
a1, b1, a2, b2 = spline(-1, 0.5, 0, 0, 3, 3)
def q1(t):
return (1 - t) * y0 + t * y1 + t *(1 - t) * ((1 - t) * a1 + t * b1)
def q2(t):
return (1 - t) * y1 + t * y2 + t *(1 - t) * ((1 - t) * a2 + t * b2)
x_list = []
y1_list = []
y2_list = []
for x in range(1000):
t = (x/1000 - x1) / (x2 - x1)
x_list.append(x)
y1_list.append(q1(t))
y2_list.append(q2(t))
plt.scatter(x_list, y1_list)
plt.scatter(x_list, y2_list)
Upvotes: 1
Views: 26
Reputation: 21
I found it! The key was the proper x-range corresponding to the 3 given points.
import numpy as np
import matplotlib.pyplot as plt
x0, y0, x1, y1, x2, y2 = -1, 0.5, 0, 0, 3, 3
def spline(x0, y0, x1, y1, x2, y2):
a11 = 2 / (x1 - x0)
a12 = 1 / (x1 - x0)
a13 = 0
a21 = a12
a22 = 2 * (1 / (x1 - x0) + 1 / (x2 - x1))
a23 = 1 / (x2 - x1)
a31 = a13
a32 = a23
a33 = 2 / (x2 - x1)
b1 = 3 * ((y1 - y0) / (x1 - x0) ** 2)
b2 = 3 * ((y1 - y0) / (x1 - x0) ** 2 + (y2 - y1) / (x2 - x1) ** 2)
b3 = 3 * ((y2 - y1) / (x2 - x1) ** 2)
A = np.array([[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]])
L = np.array([[-b1], [-b2], [-b3]])
AT = np.transpose(A)
ATA = np.matmul(AT, A)
invATA = np.linalg.inv(ATA)
ATL = np.matmul(AT, L)
X = - np.matmul(invATA, ATL)
[k0, k1, k2] = X
a1 = k0 * (x1 - x0) - (y1 - y0)
b1 = -k1 * (x1 - x0) + (y1 - y0)
a2 = k1 * (x2 - x1) - (y2 - y1)
b2 = -k2 * (x2 - x1) + (y2 - y1)
return a1[0], b1[0], a2[0], b2[0]
a1, b1, a2, b2 = spline(-1, 0.5, 0, 0, 3, 3)
def q1(t):
return (1 - t) * y0 + t * y1 + t *(1 - t) * ((1 - t) * a1 + t * b1)
def q2(t):
return (1 - t) * y1 + t * y2 + t *(1 - t) * ((1 - t) * a2 + t * b2)
x1_list = []
x2_list = []
y1_list = []
y2_list = []
for x in range(-1000, 0, 1):
t = (x/1000 - x0) / (x1 - x0)
x1_list.append(x)
y1_list.append(q1(t))
for x in range(3000):
t = (x/1000 - x1) / (x2 - x1)
x2_list.append(x)
y2_list.append(q2(t))
plt.scatter(x1_list, y1_list)
plt.scatter(x2_list, y2_list)
Upvotes: 1