[enter image description here][1]As part of a homework assignment, Im trying to create bezier curves (with 4 points) with a maximum finite curvature of K_x. Curvature is calculated with a defined formula k = |y''|/(1+y'^2)^1.5
FOr a set of 4 points where p0=(0,0) and p3=(0.3,0) and p1,p2 are the degree of freedom.
What I tried is to take the space of the square between [(0,0), (0,0.15), (0.3,0.15), (0,-0.15)] and scan it such that for every point p1 I pick (which because of symmetry I choose it only to be in the top right quarter) I scan all possible p2 and for each option I check whether the curve meets the condition of the curvature and if so I save the points that meet the condition.
import numpy as np
import scipy as sc
from matplotlib import pyplot as plt
def bezzier_curve(p0, p1, p2, p3, t):
b0 = b_in_calc(3, 0, t) * np.transpose(p0)
b1 = b_in_calc(3, 1, t) * np.transpose(p1)
b2 = b_in_calc(3, 2, t) * np.transpose(p2)
b3 = b_in_calc(3, 3, t) * np.transpose(p3)
b = b0 + b1 + b2 + b3
return b
def b_in_calc(n, i, t):
b = sc.special.binom(n, i)
t1 = np.power(1-t, n-i)
t2 = np.power(t, i)
return b*t1*t2
def numerical_deriv(y, x):
# finding numerical derivative of the finction
dx = np.gradient(x)
dy = np.gradient(y)
d = dy/dx
return d
k_max_allowed = 5
d = 0.3
p0 = np.zeros([1, 2])
p3 = np.array([[0.3, 0]])
p1 = np.zeros([1, 2])
p2 = np.zeros([1, 2])
t = np.linspace(0, 1, 10000)
P1_space = np.linspace(0, d/2, num=40)
P2_xspace = np.linspace(0, d, endpoint=False, num=40)
P2_yspace = np.linspace(-d/2, d/2, 40)
output_name = 'desired_coordinates_K_max_5.txt'
condition = 0 # if we started writing P1 or not
condition1 = 0 # if we already wrote at least 1 P2 point
with open(output_name, 'w',encoding='utf-8') as f:
for i in range(1, len(P1_space)):
p1[0, 0] = P1_space[i]
for j in range(1, len(P1_space)):
p1[0, 1] = P1_space[j]
for n in range(1, len(P2_xspace)):
p2[0,0] = P2_xspace[n]
for m in range(len(P2_yspace)):
p2[0,1] = P2_yspace[m]
# create bezier curve
bezz = bezzier_curve(p0, p1, p2, p3, t)
x_vec = bezz[0, :]
y_vec = bezz[1, :]
#fig=plt.figure(1,figsize=(6,6))
#ax=plt.subplot(111);
#plt.plot(x_vec, y_vec)
# interpolate points to spread them
x_interp = np.linspace(0, d, len(t), endpoint=True)
y_interp = np.interp(x_interp, x_vec, y_vec)
# derive and find curvature vector and max k
y_deriv = numerical_deriv(y_interp, x_interp)
y_deriv2 = numerical_deriv(y_deriv, x_interp)
k_vec = y_deriv2/(1 + y_deriv**2)**1.5
k_max = np.amax(k_vec)
# check if meets condition, if so write it to .txt file
if k_max <= k_max_allowed:
if condition == 0:
f.write('P1:' + str(P1_space[int(i)]) + ',' + str(P1_space[int(j)]) + '\n')
condition = 1
if condition1 ==0:
f.write('P2:')
condition1 = 1
f.write(str(P2_xspace[n]) + ',' + str(P2_yspace[m]) + '|')
condition = 0
condition1 = 0
f.write('\n')
f.close()
plt.text(5,5 , 'Complete', fontsize = 22)
plt.xlim(0, 15)
plt.ylim(0, 10)
plt.show()
My problem is when I read the file and try a quite a few points, I see that some of the points give curves with higher curvature then K_max I desired. I'd appreciate an advice on what conditions/constrains to add in order to fix this. maybe I should limit the space from which I can pick p2?
I would appreciate an answer with the numerical derivative because as I interpret the results I also use numerical derivative
I want to be able to make sure that the maximal curvature is less then K_max. I use a numerical derivation but the result didn't change when I did analytical derivation for the curve
Edit: The problem can be seen in the image I've added - where the control points marked in 'x' are causing a curve with large curvature near the left side. [1]: https://i.stack.imgur.com/eNeeF.png
