pszda
pszda

Reputation: 21

Why the 2D spline graph is incorrect?

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):

Spline segments

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

Answers (1)

pszda
pszda

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

Related Questions