pythoninterpolationsplineapproximation
progressive iteration approximation(PIA) method for bspline not working?
I am new to this, and tried to implement this algorithm by using uniform B-Spline. And I don't know where I did it wrong, the result just doesn't come out the way it supposed to be.
I don't know if the basis is wrong or the procedure of the PIA is wrong. Is there someone that could help me out? Thank you so much!!
I use Python to implement all this.
In my understanding, the PIA is taking the given point set, P, as control points at first(iteration 0), and use these control points to calculate a b-spline, Q. Then find the difference, d, between the P and Q. Let Q+d in each iteration, until d is small enough as the threshold you set at the beginning.
I use deboor-Cox algorithm for generating the basis matrix.
def b_spline_basis(i, p, u, nodeVector):
# p means the degree of the spline
if p == 0:
if (nodeVector[i] <= u) & (u <= nodeVector[i + 1]): # if u is between two knots, the basis would be 1
result = 1.0
else:
result = 0.0
else:
# calculate non-zero intervals
length1 = nodeVector[i + p] - nodeVector[i]
length2 = nodeVector[i + p + 1] - nodeVector[i + 1]
# calculate coefficients for basis functions
if length1 == 0: # specifically 0/0
alpha = 0
else:
alpha = (u - nodeVector[i]) / length1
if length2 == 0:
beta = 0
else:
beta = (nodeVector[i + p + 1] - u) / length2
# calculate basis functions recursively
result = alpha * b_spline_basis(i, p - 1, u, nodeVector) + beta * b_spline_basis(i + 1, p - 1, u, nodeVector)
return result
And I tried the lemniscate to test whether my implementation of PIA is okay or not.
import numpy as np
import math
from bspline import b_spline
import matplotlib.pyplot as plt
import matplotlib
from bspline_basis import b_spline_basis
matplotlib.use('TkAgg')
# lemniscate with 200 points
alpha = 1
theta = np.linspace(0, 2 * np.pi, num=200)
x_real = alpha * np.sqrt(2) * np.cos(theta) / (np.sin(theta) ** 2 + 1)
y_real = alpha * np.sqrt(2) * np.cos(theta) * np.sin(theta) / (np.sin(theta) ** 2 + 1)
# draw the real points on lemniscate
plt.scatter(x_real, y_real)
# degree of bspline is 3, number of control points is 8
degree = 3
n = 8
points = []
delta = np.linspace(0, 2 * np.pi, num=8)
# x and y are the x-axis and y-axis for the control points
x = alpha * np.sqrt(2) * np.cos(delta) / (np.sin(delta) ** 2 + 1)
y = alpha * np.sqrt(2) * np.cos(delta) * np.sin(delta) / (np.sin(delta) ** 2 + 1)
plt.scatter(x, y, color='maroon')
# calculate bspline basis matrix
def bspline_basis(n, degree, knotVector):
basis = np.zeros([n, n])
for i in range(n):
j = 0
for u in delta:
basis[i][j] = b_spline_basis(i, degree, u, knotVector)
# print('knot=', knot)
# print('basis_i=', basis, 'j=',j)
j = j + 1
return basis
a = min(delta)
b = max(delta)
knotVector = [a, a, a, a, *delta[2:-2], b, b, b, b]
# basis matrix is stored in bs
bs = bspline_basis(n, degree, knotVector)
# I think if the basis is right, this plot would be a b-spline curve, but it doesn't turn out that way. I'm also confused by this.
plt.plot(np.dot(bs, np.transpose(x)), np.dot(bs, np.transpose(y)), color='red')
# the difference between real control points with calculated value
dx = x - np.transpose(np.dot(bs, np.transpose(x)))
dy = y - np.transpose(np.dot(bs, np.transpose(y)))
# norm is going to store the norm of (dx, dy)
norm = np.zeros(n)
for i in range(len(dx)):
norm[i] = math.sqrt(dx[i] ** 2 + dy[i] ** 2)
# make the biggest norm to be the error
err = max(norm)
iteration = 0
print('iteration #', iteration, ', err = ', err)
# set the threshold for the algorithm to stop
tol = 0.2
# in while loop, calculate the difference in each loop, until error is smaller than the threshold
while err > tol:
iteration = iteration + 1
x = x + dx
y = y + dy
dx = x - np.transpose(np.dot(bs, np.transpose(x)))
dy = y - np.transpose(np.dot(bs, np.transpose(y)))
for i in range(len(dx)):
norm[i] = math.sqrt(dx[i] ** 2 + dy[i] ** 2)
err = max(norm)
print('iteration #', iteration, ', err = ', err)
x_inter = np.transpose(np.dot(bs, np.transpose(x)))
y_inter = np.transpose(np.dot(bs, np.transpose(y)))
plt.show()
But the result is not even close. The err printed in each iteration gets bigger and bigger.
iteration # 0 , err = 0.8978393078534154
iteration # 1 , err = 0.5572305648715149
iteration # 2 , err = 0.8814649114823587
iteration # 3 , err = 1.406648477874589
iteration # 4 , err = 2.2515402019886657
iteration # 5 , err = 3.610001808299592
iteration # 6 , err = 5.794725750733798
iteration # 7 , err = 9.309544995196921
iteration # 8 , err = 14.966156756400013
iteration # 9 , err = 24.072299683891867
iteration # 10 , err = 38.73507669530552
iteration # 11 , err = 62.34988787737978
iteration # 12 , err = 100.3885976037046
iteration # 13 , err = 161.67015869470632
iteration # 14 , err = 260.40916333350236
iteration # 15 , err = 419.5188341631952
iteration # 16 , err = 675.9369969104991
iteration # 17 , err = 1089.2146572938898
iteration # 18 , err = 1755.3667774904786
iteration # 19 , err = 2829.2109590140344
iteration # 20 , err = 4560.398039137244
iteration # 21 , err = 7351.530766709586
iteration # 22 , err = 11851.91790312345
iteration # 23 , err = 19108.824114848438
iteration # 24 , err = 30811.492573031916
iteration # 25 , err = 49684.87189301904
iteration # 26 , err = 80124.93280280002
iteration # 27 , err = 129223.88403951934
iteration # 28 , err = 208424.68577890267
iteration # 29 , err = 336191.3189164541
iteration # 30 , err = 542318.7082430203
iteration # 31 , err = 874889.5879288138
iteration # 32 , err = 1411504.6936387809
iteration # 33 , err = 2277412.443263706
iteration # 34 , err = 3674778.915040246
...
The printed lines are too long, I won't show them all. But you get the point.
Beside, the plot is also wierd. And I just don't know where when wrong, and I upset my for days.
Is there someone can help with this? Thank you so so much!! I'm really confused right now, hoping there is someone can help me out. TAT
Solution
There are a few things we need to take care of.
First, I will put b_spline_basis into a separate file. It was almost correct but there are two changes. The intervals where the degree zero basis functions evaluate to 1 had to be adapted so that the basis functions sum up to one on the entire interval [a, b] (your version evaluated to more in the knots). This problem happens quite often, cf. e.g. here. Also, the 0/0 case needed 1 instead of 0 for alpha and beta:
def b_spline_basis(i, p, u, knotVector):
# p means the degree of the spline
if p == 0:
# The support is closed from left but open from the right ...
if (i != len(knotVector) - 2):
if ((knotVector[i] <= u) & (u < knotVector[i + 1])):
result = 1.0
else:
result = 0.0
# ... unless it is the last one, which is closed from both sides.
else:
if ((knotVector[i] <= u) & (u <= knotVector[i + 1])):
result = 1.0
else:
result = 0.0
else:
# calculate non-zero intervals
length1 = knotVector[i + p] - knotVector[i]
length2 = knotVector[i + p + 1] - knotVector[i + 1]
# calculate coefficients for basis functions
if length1 == 0: # specifically 0/0
alpha = 1 # You had 0 here.
else:
alpha = (u - knotVector[i]) / length1
if length2 == 0:
beta = 1 # You had 0 here as well.
else:
beta = (knotVector[i + p + 1] - u) / length2
# calculate basis functions recursively
result = alpha * b_spline_basis(i, p - 1, u, knotVector) + beta * b_spline_basis(i + 1, p - 1, u, knotVector)
return result
Second, I put also bspline_basis [sic] into a separate file. It is almost identical to your version but the matrix is not necessarily square in general. I would also strongly advise to rename the function; the resulting matrix is a transpose of what is usually called collocation matrix.
import numpy as np
from b_spline_basis import b_spline_basis
# calculate bspline basis matrix
def bspline_basis(n, degree, knotVector, delta):
basis = np.zeros([n, delta.size])
for i in range(n):
j = 0
for u in delta:
basis[i][j] = b_spline_basis(i, degree, u, knotVector)
# print('knot=', knot)
# print('basis_i=', basis, 'j=',j)
j = j + 1
return basis
Finally, I throw in a function for plotting a B-spline (as a curve) given its control points etc.
import numpy as np
from b_spline_basis import b_spline_basis
def plot_bspline(plt, num_samples, degree, knotVector, x_cps, y_cps, color):
beg = knotVector[0]
end = knotVector[-1]
num_cps = len(x_cps)
x_curve = np.zeros(num_samples)
y_curve = np.zeros(num_samples)
for i in range(num_samples):
for j in range(num_cps):
t_loc = i / (num_samples-1)
t = beg * (1 - t_loc) + end * t_loc
x_curve[i] += x_cps[j] * b_spline_basis(j, degree, t, knotVector)
y_curve[i] += y_cps[j] * b_spline_basis(j, degree, t, knotVector)
plt.plot(x_curve, y_curve, color=color)
Now we get back to your code, the few corrections are commented. In general, there seemed to be three sources of confusion:
- The results of
bspline_basishad to be transposed, because they are a transposed collocation matrix. - Evaluating using
bspline_basisdoes not give you a B-spline as a curve but only its values indelta. - It is important to distinguish between
x_target, y_target(values of the lemniscate that you want to approximate) andx_cps, y_cps(B-spline control points in the current iteration). You called both of themx, y.
import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('TkAgg')
from b_spline_basis import b_spline_basis
from bspline_basis import bspline_basis
from plot_bspline import plot_bspline
# lemniscate with 200 points
alpha = 1
theta = np.linspace(0, 2 * np.pi, num=200)
x_real = alpha * np.sqrt(2) * np.cos(theta) / (np.sin(theta) ** 2 + 1)
y_real = alpha * np.sqrt(2) * np.cos(theta) * np.sin(theta) / (np.sin(theta) ** 2 + 1)
# draw the real points on lemniscate
plt.scatter(x_real, y_real)
# degree of bspline is 3, number of control points is 8
degree = 3
n = 8
points = []
delta = np.linspace(0, 2 * np.pi, num=8)
# x_target and y_target are the values we want to approximate.
# They will be used as starting values for the control points as well.
x_target = alpha * np.sqrt(2) * np.cos(delta) / (np.sin(delta) ** 2 + 1)
y_target = alpha * np.sqrt(2) * np.cos(delta) * np.sin(delta) / (np.sin(delta) ** 2 + 1)
plt.scatter(x_target, y_target, color='maroon')
a = min(delta)
b = max(delta)
knotVector = [a, a, a, a, *delta[2:-2], b, b, b, b]
# basis matrix is stored in bs
bs = bspline_basis(n, degree, knotVector, delta)
# I think if the basis is right, this plot would be a b-spline curve, but it doesn't turn out that way. I'm also confused by this.
# The transpositions were wrong.
# Also, using bs does not give you a B-spline as a curve but only its values evaluated at delta, i.e., at 8 points.
plt.plot(np.dot(np.transpose(bs), x_target), np.dot(np.transpose(bs), y_target), color='red')
# If you also plot the B-spline as curve that uses x_target and y_target as control points, you will see that the red curve connects 8 of its values.
plot_bspline(plt, 100, 3, knotVector, x_target, y_target, 'green')
# Now to PIA.
# The control points in the first iteration will be the initial values.
x_cps = x_target
y_cps = y_target
# Then we have a difference between the target values and the corresponding values of our B-spline.
dx = x_target - np.transpose(np.dot(np.transpose(bs), x_cps))
dy = y_target - np.transpose(np.dot(np.transpose(bs), y_cps))
# norm is going to store the norm of (dx, dy)
norm = np.zeros(n)
for i in range(len(dx)):
norm[i] = math.sqrt(dx[i] ** 2 + dy[i] ** 2)
# make the biggest norm to be the error
err = max(norm)
iteration = 0
print('iteration #', iteration, ', err = ', err)
# set the threshold for the algorithm to stop
tol = 1e-5
# in while loop, calculate the difference in each loop, until error is smaller than the threshold
while err > tol and iteration < 100:
iteration = iteration + 1
# We change the control points ...
x_cps = x_cps + dx
y_cps = y_cps + dy
# ... and compute the difference from the target (which is constant)!
dx = x_target - np.transpose(np.dot(np.transpose(bs), x_cps))
dy = y_target - np.transpose(np.dot(np.transpose(bs), y_cps))
for i in range(len(dx)):
norm[i] = math.sqrt(dx[i] ** 2 + dy[i] ** 2)
err = max(norm)
print('iteration #', iteration, ', err = ', err)
x_inter = np.transpose(np.dot(np.transpose(bs), x_cps))
y_inter = np.transpose(np.dot(np.transpose(bs), y_cps))
# If I plot the way you did, I will again not get a B-spline as a curve but the values of the B-spline evaluated at delta. Notice that it passes the maroon points.
plt.plot(x_inter, y_inter, color='yellow')
# Let's now plot the entire B-spline as a curve. Notice that it passes through the maroon points.
plot_bspline(plt, 100, 3, knotVector, x_cps, y_cps, 'magenta')
plt.show()
If we now have a look at the plot, there is quite a lot happening:
- The blue points are the samples of the lemniscate, i.e., the input.
- The maroon points are the eight points on the lemniscate that we will be approximating (there seem to be only seven, since the first and last ones coincide).
- The green curve is the initial guess, i.e., a B-spline that uses the maroon points as its control points.
- The red polygon uses
bspline_basisto connect the values of the green B-spline in the parameter values indelta. This is a corrected version of your red curve. - The magenta curve is the final guess, i.e., a B-spline that approximates the maroon points up to
tol. - The yellow curve uses
bspline_basisto connect the values of the magenta B-spline alongdelta. This is a corrected version of yourx_inter, y_inter.
I wish you good luck with further B-spline experiments. If you are also into neural networks, you might enjoy a recent paper by my friends that investigated the connection between LSPIA and gradient descent:
Dany Rios and Bert Jüttler: LSPIA,(stochastic) gradient descent, and parameter correction. Journal of Computational and Applied Mathematics 406 (2022): 113921 (preprint)