如何使用 geomdl 或其他库向样条曲线添加边界约束？
How to add boundary constraints to a spline with geomdl or other library?
这是没有约束的样条：
from geomdl import fitting
from geomdl.visualization import VisMPL
path = [(2077.0, 712.0, 1136.6176470588234), (2077.0004154771536, 974.630482962754, 1313.735294117647), (2077.1630960823995, 1302.460574562254, 1490.8529411764707), (2078.1944091179635, 1674.693193015173, 1667.9705882352941), (2080.5096120056783, 2086.976611915444, 1845.0882352941176), (2085.1051468332066, 2711.054258877495, 2022.2058823529412), (1477.0846185328733, 2803.6223679691457, 2199.323529411765), (948.4693105162195, 2802.0390667447105, 2376.4411764705883), (383.8615403256207, 2804.843424134807, 2553.5588235294117), (-41.6669725172834, 2497.067373170676, 2730.676470588235), (-37.94311919744064, 1970.5155845437525, 2907.794117647059), (-35.97395938535092, 1576.713103381243, 3084.9117647058824), (-35.125016151504795, 1214.2319876178394, 3262.029411764706), (-35.000550767864524, 893.3910350913443, 3439.1470588235297), (-35.0, 631.2108462417168, 3616.264705882353), (-35.0, 365.60545190581837, 3793.3823529411766), (-35.0, 100.00005756991993, 3970.5)]
degree = 3
curve = fitting.interpolate_curve(path, degree)
curve.vis = VisMPL.VisCurve3D()
curve.render()
# the following is to show it under matplotlib and prepare solutions comparison
import numpy as np
import matplotlib.pyplot as plt
qtPoints = 3*len(path)
s = np.linspace(0, 1, qtPoints, True).tolist()
pt = curve.tangent(s) # returns points and tangents
spline = [u for u, v in pt] # get points, leave tangents
我要添加约束：
- x >= -35
- x <= 2077
- y <= 2802
geomdl 库不建议带约束的样条。我试过这个 hack，只是通过纠正点以保持在边界内：
path2 = [(x if x >= -35 else -35, y if y <= 2802 else 2802, z) for x,y,z in spline]
path2 = [(x if x <= 2077 else 2077, y, z) for x,y,z in path2]
curve2 = fitting.interpolate_curve(path2, 3)
pt2 = curve2.tangent(s) # returns points and tangents
spline2 = [u for u, v in pt2] # get points, leave tangents
plt.plot([u[0] for u in path], [u[1] for u in path], 'o',
[u[0] for u in spline], [u[1] for u in spline], 'b',
[u[0] for u in spline2], [u[1] for u in spline2], 'r')
plt.show()
curve2.vis = VisMPL.VisCurve3D()
curve2.render()
这两个在一起（向左转90°）：
结果不理想（红色）：
另一种方法是直接使用路径作为控制点。这是 NURBS 的结果：
from geomdl import NURBS
curve_n = NURBS.Curve()
curve_n.degree = min(degree, len(path)) # order = degree+1
curve_n.ctrlpts = path
last_knot = len(path) - curve_n.degree
curve_n.knotvector = np.concatenate((np.zeros(curve_n.degree), np.arange(0, last_knot + 1), np.ones(curve_n.degree)*last_knot)).astype(int)
curve_n.delta = 0.05
spline_n = curve_n.evalpts
plt.plot([u[0] for u in path], [u[1] for u in path], 'o',
[u[0] for u in spline_f], [u[1] for u in spline_f], 'b',
[u[0] for u in spline2], [u[1] for u in spline2], 'r',
[u[0] for u in spline_n], [u[1] for u in spline_n], 'g')
plt.show()
结果（绿色）离路径太远。
如果我使用 NURBS 点进行新的拟合，并玩弄 NURBS 度数，我会得到满意的结果：
from geomdl import fitting
from geomdl import NURBS
#from geomdl.visualization import VisMPL
import numpy as np
import matplotlib.pyplot as plt
path = [(2077.0, 712.0, 1136.6176470588234), (2077.0004154771536, 974.630482962754, 1313.735294117647), (2077.1630960823995, 1302.460574562254, 1490.8529411764707), (2078.1944091179635, 1674.693193015173, 1667.9705882352941), (2080.5096120056783, 2086.976611915444, 1845.0882352941176), (2085.1051468332066, 2711.054258877495, 2022.2058823529412), (1477.0846185328733, 2803.6223679691457, 2199.323529411765), (948.4693105162195, 2802.0390667447105, 2376.4411764705883), (383.8615403256207, 2804.843424134807, 2553.5588235294117), (-41.6669725172834, 2497.067373170676, 2730.676470588235), (-37.94311919744064, 1970.5155845437525, 2907.794117647059), (-35.97395938535092, 1576.713103381243, 3084.9117647058824), (-35.125016151504795, 1214.2319876178394, 3262.029411764706), (-35.000550767864524, 893.3910350913443, 3439.1470588235297), (-35.0, 631.2108462417168, 3616.264705882353), (-35.0, 365.60545190581837, 3793.3823529411766), (-35.0, 100.00005756991993, 3970.5)]
degree = 3
qtPoints = 3*len(path)
# fitting without constraints
curve_f = fitting.interpolate_curve(path, degree)
#curve.vis = VisMPL.VisCurve3D()
#curve.render()
s = np.linspace(0, 1, qtPoints, True).tolist()
pt = curve_f.tangent(s) # returns points and tangents
spline = [u for u, v in pt] # get points, leave tangents
# fitting with constraints, awkward hack
path2 = [(x if x >= -35 else -35, y if y <= 2802 else 2802, z) for x,y,z in spline]
path2 = [(x if x <= 2077 else 2077, y, z) for x,y,z in path2]
curve2 = fitting.interpolate_curve(path2, 3)
pt2 = curve2.tangent(s) # returns points and tangents
spline2 = [u for u, v in pt2] # get points, leave tangents
# control points = path
curve_n = NURBS.Curve()
curve_n.degree = 2 #min(degree, len(path)) # order = degree+1
curve_n.ctrlpts = path
last_knot = len(path) - curve_n.degree
curve_n.knotvector = np.concatenate((np.zeros(curve_n.degree), np.arange(0, last_knot + 1), np.ones(curve_n.degree)*last_knot)).astype(int)
curve_n.delta = 0.05
spline_n = curve_n.evalpts
# fitting without constraints on NURBS points
curve3 = fitting.interpolate_curve(spline_n, 3)
pt3 = curve3.tangent(s) # returns points and tangents
spline3 = [u for u, v in pt3] # get points, leave tangents
plt.plot([u[0] for u in path], [u[1] for u in path], 'o',
[u[0] for u in spline_f], [u[1] for u in spline_f], 'b',
[u[0] for u in spline2], [u[1] for u in spline2], 'r',
[u[0] for u in spline3], [u[1] for u in spline3], 'y',
[u[0] for u in spline_n], [u[1] for u in spline_n], 'g')
plt.show()
但它并不坚固，可能只是一个臭名昭著的 DIY。
[True if x >= -35 and x <= 2077 and y <= 2802 else False for x,y,z in spline3]
[True, False, False, False, False, False, False, False, False, False, False, False, False, False, True, True, True, True, True, True, True, False, False, False, False, True, True, True, True, True, True, True, True, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, True, True, True]
如何保持顺畅，在路径上，请尊重约束，可能与另一个图书馆？我找到 , but that solves derivatives constraints and I don't figure out how to adapt this solution. I raised also the question on a strictly mathematical point of view here.
好吧，这个话题很难，但我明白了，受到 二维边界（导数）约束样条的启发。提议的解决方案还使用 scipy.optimize.minimize.
这是完整的代码，经过一些解释：
import numpy as np
from scipy.interpolate import UnivariateSpline, splev, splprep, BSpline
from scipy.optimize import minimize
xmin = -35
xmax = 2077
ymax = 2802
def guess(p, k, s, w=None):
"""Do an ordinary spline fit to provide knots"""
return splprep(p, w, k=k, s=s)
def err(c, p, u, t, c_shape, k, w=None):
"""The error function to minimize"""
diff = (np.array(p) - splev(u, (t, c.reshape(c_shape), k))).flatten()
if w is None:
diff = (diff*diff).sum()
else:
diff = (diff*diff*w).sum() #not sure it is the good way to multiply w
return np.abs(diff)
def constraint(c, l, t, c_shape, k, eqorineq, eqinterv):
X = np.linspace(0, 1, l*20)
v = splev(X, (t, c.reshape(c_shape), k))
if eqorineq == 'ineq':
ineq_contrib = sum([(x < xmin)*(x-xmin)**2 + (x > xmax)*(x-xmax)**2 for x in v[0]] \
+ [(y > ymax)*(y-ymax)**2 for y in v[1]])
eq_contrib = 0
for i in range(len(X)):
eq_contrib += (X[i] >= eqinterv[0][0] and X[i] <= eqinterv[0][1]) * (v[0][i] - xmin)**2 \
+ (X[i] >= eqinterv[1][0] and X[i] <= eqinterv[1][1]) * (v[0][i] - xmax)**2 \
+ (X[i] >= eqinterv[2][0] and X[i] <= eqinterv[2][1]) * (v[1][i] - ymax)**2
return -(ineq_contrib + eq_contrib)
# return -1 * ineq_contrib
elif eqorineq == 'eq':
res = 0 # equality
for i in range(len(X)):
if X[i] >= eqinterv[0][0] and X[i] <= eqinterv[0][1] and v[0][i] != xmin \
or X[i] >= eqinterv[1][0] and X[i] <= eqinterv[1][1] and v[0][i] != xmax \
or X[i] >= eqinterv[2][0] and X[i] <= eqinterv[2][1] and v[1][i] != ymax :
res = 1
return res
def spline_neumann(p, k=3, s=0, w=None):
tck, u = guess(p, k, s, w=w)
t, c0, k = tck
c0flat = np.array(c0).flatten()
c_shape = np.array(c0).shape
x0 = 0 #x[0] # point at which zero slope is required
# compute u intervals for eq constraints
xmin_umin = xmin_umax = xmax_umin = xmax_umax = ymax_umin = ymax_umax = -1
for i in range(len(p[0])):
if xmin_umin == -1 and p[0][i] <= xmin : xmin_umin = u[i]
if xmin_umin != -1 and xmin_umax == -1 and p[0][i] > xmin : xmin_umax = u[i-1]
if xmax_umin == -1 and p[0][i] >= xmax : xmax_umin = u[i]
if xmax_umin != -1 and xmax_umax == -1 and p[0][i] < xmax : xmax_umax = u[i-1]
if ymax_umin == -1 and p[1][i] >= ymax : ymax_umin = u[i]
if ymax_umin != -1 and ymax_umax == -1 and p[1][i] < ymax : ymax_umax = u[i-1]
eqinterv = [[xmin_umin, xmin_umax], [xmax_umin, xmax_umax], [ymax_umin, ymax_umax]]
for i in range(len(eqinterv)):
if eqinterv[i][0] == -1 : eqinterv[i][0] = 0
if eqinterv[i][1] == -1 : eqinterv[i][1] = 1
print("eqinterv = ", eqinterv)
con = {'type': 'ineq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'ineq', eqinterv)
#'type': 'eq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'eq', eqinterv)
#'fun': lambda c: splev(x0, (t, c.reshape(c_shape), k), der=1),
#'jac': lambda c: splev(x0, (t, c, k), der=2) # doesn't help, dunno why
}
opt = minimize(err, c0flat, (p, u, t, c_shape, k, w), constraints=con)
#opt = minimize(err, c0, (p, u, t, c_shape, k, w), method='Nelder-Mead', constraints=con)
#opt = minimize(err, c0flat, (p, u, t, c_shape, k, w))
copt = opt.x.reshape(c_shape)
#return UnivariateSpline._from_tck((t, copt, k))
#return BSpline(t, k, copt)
return ((t, copt, k), opt.success)
import matplotlib.pyplot as plt
path = [(2077.0, 712.0, 1136.6176470588234), (2077.0004154771536, 974.630482962754, 1313.735294117647), (2077.1630960823995, 1302.460574562254, 1490.8529411764707), (2078.1944091179635, 1674.693193015173, 1667.9705882352941), (2080.5096120056783, 2086.976611915444, 1845.0882352941176), (2085.1051468332066, 2711.054258877495, 2022.2058823529412), (1477.0846185328733, 2803.6223679691457, 2199.323529411765), (948.4693105162195, 2802.0390667447105, 2376.4411764705883), (383.8615403256207, 2804.843424134807, 2553.5588235294117), (-41.6669725172834, 2497.067373170676, 2730.676470588235), (-37.94311919744064, 1970.5155845437525, 2907.794117647059), (-35.97395938535092, 1576.713103381243, 3084.9117647058824), (-35.125016151504795, 1214.2319876178394, 3262.029411764706), (-35.000550767864524, 893.3910350913443, 3439.1470588235297), (-35.0, 631.2108462417168, 3616.264705882353), (-35.0, 365.60545190581837, 3793.3823529411766), (-35.0, 100.00005756991993, 3970.5)]
pathxyz = [[x for x,y,z in path], [y for x,y,z in path], [z for x,y,z in path]]
n = len(path)
#std would be interesting to define as the standard deviation of the curve compared to a no noise one. No noise ==> s=0
k = 5
s = 0
sp0, u = guess(pathxyz, k, s)
sp, success = spline_neumann(pathxyz, k, s) #s=n*std
print("success = ", success)
# % of points not respecting the constraints
perfo_vs_ineq = (sum([(x < xmin) for x in v[0]]) + sum([(x > xmax) for x in v[0]]) + sum([(y > ymax) for y in v[1]]) )/len(v[0])/2
print("perfo% vs ineq constraints = ", perfo_vs_ineq)
X = np.linspace(0, 1, len(pathxyz)*10)
val0 = splev(X, sp0)
val = splev(X, sp)
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot([x for x,y,z in path], [y for x,y,z in path], [z for x,y,z in path], 'ro')
ax.plot(val0[0], val0[1], val0[2], 'b-')
ax.plot(val[0], val[1], val[2], 'r-')
plt.show()
plt.figure()
plt.plot(val0[0], val0[1], '-b', lw=1, label='guess')
plt.plot(val[0], val[1], '-r', lw=2, label='spline')
plt.plot(pathxyz[0], pathxyz[1], 'ok', label='data')
plt.legend(loc='best')
plt.show()
最后，我有2D和3D渲染图。 3D 视图显示样条曲线使用 z 轴进行平滑。这对我的用例来说并不令人满意，所以我必须在我的约束中考虑到它，但这超出了本 Q/A:
的范围
以及显示样条约束效果的 2D 视图：
蓝色曲线没有约束，红色曲线有约束。
现在对前进方向的解释：
- 没有约束的样条计算如下：
sp0, u = guess(pathxyz, k, s)
- 带约束的样条计算如下：
sp, success = spline_neumann(pathxyz, k, s)
- 然后打印
success 遵循 scipy.optimize.minimize 标准和基于不等式约束的自定义标准作为不满足约束的点的百分比：
print("success = ", success)
perfo_vs_ineq = (sum([(x < xmin) for x in v[0]]) + sum([(x > xmax) for x in v[0]]) + sum([(y > ymax) for y in v[1]]) )/len(v[0])/2
print("perfo% vs ineq constraints = ", perfo_vs_ineq)
- 约束下的优化执行者：
opt = minimize(err, c0flat, (p, u, t, c_shape, k, w), constraints=con)。它优化了由非约束求解 获得的c0flat初始化样条的系数
- return我们必须重塑
copt = opt.x 中的系数才能被 splev 和 copt = opt.x.reshape(c_shape) 使用
splev 用于评估两个样条曲线 - 不受约束和受约束 - 这里没有新内容：
X = np.linspace(0, 1, len(pathxyz)*10)
val0 = splev(X, sp0)
val = splev(X, sp)
scipy.optimize.minimize 参数和 return 值在 manual 中进行了解释。所以我在这里只解释具体的内容err 是要最小化的 cost。它旨在坚持控制点：
def err(c, p, u, t, c_shape, k, w=None):
"""The error function to minimize"""
diff = (np.array(p) - splev(u, (t, c.reshape(c_shape), k))).flatten()
if w is None:
diff = (diff*diff).sum()
else:
diff = (diff*diff*w).sum() #not sure it is the good way to multiply w
return np.abs(diff)
- 我没有用重量参数进行测试
w。这里需要理解的重要一点是，我们仅使用 u 提供的曲线坐标对控制点进行评估。成本是控制点与 scipy.optimize.minimize 尝试使用新计算系数 c 评估它们的方式之间的差异
- 然后我们得出
constraints=con 提供给 scipy.optimize.minimize 的约束定义为：
con = {'type': 'ineq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'ineq', eqinterv)
#'type': 'eq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'eq', eqinterv)
- 我只使用
ineq 特性，因为使用 eq 特性的测试在我的用例中给出了很差的结果，但如果它可以帮助某人，我已经让代码。因此，不等式和等式约束都是使用函数 constraint(c, len(p[0]), t, c_shape, k, 'ineq', eqinterv) 计算的。我更喜欢一个函数而不是一个函数列表来只执行一次样条评估。所以当然，c是被scipy.optimize.minimize评估的参数，t和k完成评估所需的(t,c,k)元组，len(p[0])是与要评估的点数成比例，'ineq' 告诉 constraint 处理不等式，eqinterv 是一个向量，我想在其中评估计算为成本的等式。在我的用例中，我记得我需要 x >= -35 and x <= 2077 and y <= 2802。我没有详细说明特定于我的用例的计算，我只强调间隔与与 u: 齐次的曲线坐标相关的点
xmin_umin = xmin_umax = xmax_umin = xmax_umax = ymax_umin = ymax_umax = -1
for i in range(len(p[0])):
if xmin_umin == -1 and p[0][i] <= xmin : xmin_umin = u[i]
if xmin_umin != -1 and xmin_umax == -1 and p[0][i] > xmin : xmin_umax = u[i-1]
if xmax_umin == -1 and p[0][i] >= xmax : xmax_umin = u[i]
if xmax_umin != -1 and xmax_umax == -1 and p[0][i] < xmax : xmax_umax = u[i-1]
if ymax_umin == -1 and p[1][i] >= ymax : ymax_umin = u[i]
if ymax_umin != -1 and ymax_umax == -1 and p[1][i] < ymax : ymax_umax = u[i-1]
- 然后计算等式的成本：
eq_contrib = 0
for i in range(len(X)):
eq_contrib += (X[i] >= eqinterv[0][0] and X[i] <= eqinterv[0][1]) * (v[0][i] - xmin)**2 \
+ (X[i] >= eqinterv[1][0] and X[i] <= eqinterv[1][1]) * (v[0][i] - xmax)**2 \
+ (X[i] >= eqinterv[2][0] and X[i] <= eqinterv[2][1]) * (v[1][i] - ymax)**2
- 不平等成本很简单：
ineq_contrib = sum([(x < xmin)*(x-xmin)**2 + (x > xmax)*(x-xmax)**2 for x in v[0]] \
+ [(y > ymax)*(y-ymax)**2 for y in v[1]])
就这些了，希望有用
这是没有约束的样条：
from geomdl import fitting
from geomdl.visualization import VisMPL
path = [(2077.0, 712.0, 1136.6176470588234), (2077.0004154771536, 974.630482962754, 1313.735294117647), (2077.1630960823995, 1302.460574562254, 1490.8529411764707), (2078.1944091179635, 1674.693193015173, 1667.9705882352941), (2080.5096120056783, 2086.976611915444, 1845.0882352941176), (2085.1051468332066, 2711.054258877495, 2022.2058823529412), (1477.0846185328733, 2803.6223679691457, 2199.323529411765), (948.4693105162195, 2802.0390667447105, 2376.4411764705883), (383.8615403256207, 2804.843424134807, 2553.5588235294117), (-41.6669725172834, 2497.067373170676, 2730.676470588235), (-37.94311919744064, 1970.5155845437525, 2907.794117647059), (-35.97395938535092, 1576.713103381243, 3084.9117647058824), (-35.125016151504795, 1214.2319876178394, 3262.029411764706), (-35.000550767864524, 893.3910350913443, 3439.1470588235297), (-35.0, 631.2108462417168, 3616.264705882353), (-35.0, 365.60545190581837, 3793.3823529411766), (-35.0, 100.00005756991993, 3970.5)]
degree = 3
curve = fitting.interpolate_curve(path, degree)
curve.vis = VisMPL.VisCurve3D()
curve.render()
# the following is to show it under matplotlib and prepare solutions comparison
import numpy as np
import matplotlib.pyplot as plt
qtPoints = 3*len(path)
s = np.linspace(0, 1, qtPoints, True).tolist()
pt = curve.tangent(s) # returns points and tangents
spline = [u for u, v in pt] # get points, leave tangents
我要添加约束：
- x >= -35
- x <= 2077
- y <= 2802
geomdl 库不建议带约束的样条。我试过这个 hack，只是通过纠正点以保持在边界内：
path2 = [(x if x >= -35 else -35, y if y <= 2802 else 2802, z) for x,y,z in spline]
path2 = [(x if x <= 2077 else 2077, y, z) for x,y,z in path2]
curve2 = fitting.interpolate_curve(path2, 3)
pt2 = curve2.tangent(s) # returns points and tangents
spline2 = [u for u, v in pt2] # get points, leave tangents
plt.plot([u[0] for u in path], [u[1] for u in path], 'o',
[u[0] for u in spline], [u[1] for u in spline], 'b',
[u[0] for u in spline2], [u[1] for u in spline2], 'r')
plt.show()
curve2.vis = VisMPL.VisCurve3D()
curve2.render()
这两个在一起（向左转90°）：
结果不理想（红色）：
另一种方法是直接使用路径作为控制点。这是 NURBS 的结果：
from geomdl import NURBS
curve_n = NURBS.Curve()
curve_n.degree = min(degree, len(path)) # order = degree+1
curve_n.ctrlpts = path
last_knot = len(path) - curve_n.degree
curve_n.knotvector = np.concatenate((np.zeros(curve_n.degree), np.arange(0, last_knot + 1), np.ones(curve_n.degree)*last_knot)).astype(int)
curve_n.delta = 0.05
spline_n = curve_n.evalpts
plt.plot([u[0] for u in path], [u[1] for u in path], 'o',
[u[0] for u in spline_f], [u[1] for u in spline_f], 'b',
[u[0] for u in spline2], [u[1] for u in spline2], 'r',
[u[0] for u in spline_n], [u[1] for u in spline_n], 'g')
plt.show()
结果（绿色）离路径太远。
如果我使用 NURBS 点进行新的拟合，并玩弄 NURBS 度数，我会得到满意的结果：
from geomdl import fitting
from geomdl import NURBS
#from geomdl.visualization import VisMPL
import numpy as np
import matplotlib.pyplot as plt
path = [(2077.0, 712.0, 1136.6176470588234), (2077.0004154771536, 974.630482962754, 1313.735294117647), (2077.1630960823995, 1302.460574562254, 1490.8529411764707), (2078.1944091179635, 1674.693193015173, 1667.9705882352941), (2080.5096120056783, 2086.976611915444, 1845.0882352941176), (2085.1051468332066, 2711.054258877495, 2022.2058823529412), (1477.0846185328733, 2803.6223679691457, 2199.323529411765), (948.4693105162195, 2802.0390667447105, 2376.4411764705883), (383.8615403256207, 2804.843424134807, 2553.5588235294117), (-41.6669725172834, 2497.067373170676, 2730.676470588235), (-37.94311919744064, 1970.5155845437525, 2907.794117647059), (-35.97395938535092, 1576.713103381243, 3084.9117647058824), (-35.125016151504795, 1214.2319876178394, 3262.029411764706), (-35.000550767864524, 893.3910350913443, 3439.1470588235297), (-35.0, 631.2108462417168, 3616.264705882353), (-35.0, 365.60545190581837, 3793.3823529411766), (-35.0, 100.00005756991993, 3970.5)]
degree = 3
qtPoints = 3*len(path)
# fitting without constraints
curve_f = fitting.interpolate_curve(path, degree)
#curve.vis = VisMPL.VisCurve3D()
#curve.render()
s = np.linspace(0, 1, qtPoints, True).tolist()
pt = curve_f.tangent(s) # returns points and tangents
spline = [u for u, v in pt] # get points, leave tangents
# fitting with constraints, awkward hack
path2 = [(x if x >= -35 else -35, y if y <= 2802 else 2802, z) for x,y,z in spline]
path2 = [(x if x <= 2077 else 2077, y, z) for x,y,z in path2]
curve2 = fitting.interpolate_curve(path2, 3)
pt2 = curve2.tangent(s) # returns points and tangents
spline2 = [u for u, v in pt2] # get points, leave tangents
# control points = path
curve_n = NURBS.Curve()
curve_n.degree = 2 #min(degree, len(path)) # order = degree+1
curve_n.ctrlpts = path
last_knot = len(path) - curve_n.degree
curve_n.knotvector = np.concatenate((np.zeros(curve_n.degree), np.arange(0, last_knot + 1), np.ones(curve_n.degree)*last_knot)).astype(int)
curve_n.delta = 0.05
spline_n = curve_n.evalpts
# fitting without constraints on NURBS points
curve3 = fitting.interpolate_curve(spline_n, 3)
pt3 = curve3.tangent(s) # returns points and tangents
spline3 = [u for u, v in pt3] # get points, leave tangents
plt.plot([u[0] for u in path], [u[1] for u in path], 'o',
[u[0] for u in spline_f], [u[1] for u in spline_f], 'b',
[u[0] for u in spline2], [u[1] for u in spline2], 'r',
[u[0] for u in spline3], [u[1] for u in spline3], 'y',
[u[0] for u in spline_n], [u[1] for u in spline_n], 'g')
plt.show()
但它并不坚固，可能只是一个臭名昭著的 DIY。
[True if x >= -35 and x <= 2077 and y <= 2802 else False for x,y,z in spline3]
[True, False, False, False, False, False, False, False, False, False, False, False, False, False, True, True, True, True, True, True, True, False, False, False, False, True, True, True, True, True, True, True, True, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, True, True, True]
如何保持顺畅，在路径上，请尊重约束，可能与另一个图书馆？我找到
好吧，这个话题很难，但我明白了，受到 scipy.optimize.minimize.
这是完整的代码，经过一些解释：
import numpy as np
from scipy.interpolate import UnivariateSpline, splev, splprep, BSpline
from scipy.optimize import minimize
xmin = -35
xmax = 2077
ymax = 2802
def guess(p, k, s, w=None):
"""Do an ordinary spline fit to provide knots"""
return splprep(p, w, k=k, s=s)
def err(c, p, u, t, c_shape, k, w=None):
"""The error function to minimize"""
diff = (np.array(p) - splev(u, (t, c.reshape(c_shape), k))).flatten()
if w is None:
diff = (diff*diff).sum()
else:
diff = (diff*diff*w).sum() #not sure it is the good way to multiply w
return np.abs(diff)
def constraint(c, l, t, c_shape, k, eqorineq, eqinterv):
X = np.linspace(0, 1, l*20)
v = splev(X, (t, c.reshape(c_shape), k))
if eqorineq == 'ineq':
ineq_contrib = sum([(x < xmin)*(x-xmin)**2 + (x > xmax)*(x-xmax)**2 for x in v[0]] \
+ [(y > ymax)*(y-ymax)**2 for y in v[1]])
eq_contrib = 0
for i in range(len(X)):
eq_contrib += (X[i] >= eqinterv[0][0] and X[i] <= eqinterv[0][1]) * (v[0][i] - xmin)**2 \
+ (X[i] >= eqinterv[1][0] and X[i] <= eqinterv[1][1]) * (v[0][i] - xmax)**2 \
+ (X[i] >= eqinterv[2][0] and X[i] <= eqinterv[2][1]) * (v[1][i] - ymax)**2
return -(ineq_contrib + eq_contrib)
# return -1 * ineq_contrib
elif eqorineq == 'eq':
res = 0 # equality
for i in range(len(X)):
if X[i] >= eqinterv[0][0] and X[i] <= eqinterv[0][1] and v[0][i] != xmin \
or X[i] >= eqinterv[1][0] and X[i] <= eqinterv[1][1] and v[0][i] != xmax \
or X[i] >= eqinterv[2][0] and X[i] <= eqinterv[2][1] and v[1][i] != ymax :
res = 1
return res
def spline_neumann(p, k=3, s=0, w=None):
tck, u = guess(p, k, s, w=w)
t, c0, k = tck
c0flat = np.array(c0).flatten()
c_shape = np.array(c0).shape
x0 = 0 #x[0] # point at which zero slope is required
# compute u intervals for eq constraints
xmin_umin = xmin_umax = xmax_umin = xmax_umax = ymax_umin = ymax_umax = -1
for i in range(len(p[0])):
if xmin_umin == -1 and p[0][i] <= xmin : xmin_umin = u[i]
if xmin_umin != -1 and xmin_umax == -1 and p[0][i] > xmin : xmin_umax = u[i-1]
if xmax_umin == -1 and p[0][i] >= xmax : xmax_umin = u[i]
if xmax_umin != -1 and xmax_umax == -1 and p[0][i] < xmax : xmax_umax = u[i-1]
if ymax_umin == -1 and p[1][i] >= ymax : ymax_umin = u[i]
if ymax_umin != -1 and ymax_umax == -1 and p[1][i] < ymax : ymax_umax = u[i-1]
eqinterv = [[xmin_umin, xmin_umax], [xmax_umin, xmax_umax], [ymax_umin, ymax_umax]]
for i in range(len(eqinterv)):
if eqinterv[i][0] == -1 : eqinterv[i][0] = 0
if eqinterv[i][1] == -1 : eqinterv[i][1] = 1
print("eqinterv = ", eqinterv)
con = {'type': 'ineq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'ineq', eqinterv)
#'type': 'eq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'eq', eqinterv)
#'fun': lambda c: splev(x0, (t, c.reshape(c_shape), k), der=1),
#'jac': lambda c: splev(x0, (t, c, k), der=2) # doesn't help, dunno why
}
opt = minimize(err, c0flat, (p, u, t, c_shape, k, w), constraints=con)
#opt = minimize(err, c0, (p, u, t, c_shape, k, w), method='Nelder-Mead', constraints=con)
#opt = minimize(err, c0flat, (p, u, t, c_shape, k, w))
copt = opt.x.reshape(c_shape)
#return UnivariateSpline._from_tck((t, copt, k))
#return BSpline(t, k, copt)
return ((t, copt, k), opt.success)
import matplotlib.pyplot as plt
path = [(2077.0, 712.0, 1136.6176470588234), (2077.0004154771536, 974.630482962754, 1313.735294117647), (2077.1630960823995, 1302.460574562254, 1490.8529411764707), (2078.1944091179635, 1674.693193015173, 1667.9705882352941), (2080.5096120056783, 2086.976611915444, 1845.0882352941176), (2085.1051468332066, 2711.054258877495, 2022.2058823529412), (1477.0846185328733, 2803.6223679691457, 2199.323529411765), (948.4693105162195, 2802.0390667447105, 2376.4411764705883), (383.8615403256207, 2804.843424134807, 2553.5588235294117), (-41.6669725172834, 2497.067373170676, 2730.676470588235), (-37.94311919744064, 1970.5155845437525, 2907.794117647059), (-35.97395938535092, 1576.713103381243, 3084.9117647058824), (-35.125016151504795, 1214.2319876178394, 3262.029411764706), (-35.000550767864524, 893.3910350913443, 3439.1470588235297), (-35.0, 631.2108462417168, 3616.264705882353), (-35.0, 365.60545190581837, 3793.3823529411766), (-35.0, 100.00005756991993, 3970.5)]
pathxyz = [[x for x,y,z in path], [y for x,y,z in path], [z for x,y,z in path]]
n = len(path)
#std would be interesting to define as the standard deviation of the curve compared to a no noise one. No noise ==> s=0
k = 5
s = 0
sp0, u = guess(pathxyz, k, s)
sp, success = spline_neumann(pathxyz, k, s) #s=n*std
print("success = ", success)
# % of points not respecting the constraints
perfo_vs_ineq = (sum([(x < xmin) for x in v[0]]) + sum([(x > xmax) for x in v[0]]) + sum([(y > ymax) for y in v[1]]) )/len(v[0])/2
print("perfo% vs ineq constraints = ", perfo_vs_ineq)
X = np.linspace(0, 1, len(pathxyz)*10)
val0 = splev(X, sp0)
val = splev(X, sp)
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot([x for x,y,z in path], [y for x,y,z in path], [z for x,y,z in path], 'ro')
ax.plot(val0[0], val0[1], val0[2], 'b-')
ax.plot(val[0], val[1], val[2], 'r-')
plt.show()
plt.figure()
plt.plot(val0[0], val0[1], '-b', lw=1, label='guess')
plt.plot(val[0], val[1], '-r', lw=2, label='spline')
plt.plot(pathxyz[0], pathxyz[1], 'ok', label='data')
plt.legend(loc='best')
plt.show()
最后，我有2D和3D渲染图。 3D 视图显示样条曲线使用 z 轴进行平滑。这对我的用例来说并不令人满意，所以我必须在我的约束中考虑到它，但这超出了本 Q/A:
的范围以及显示样条约束效果的 2D 视图：
蓝色曲线没有约束，红色曲线有约束。
现在对前进方向的解释：
- 没有约束的样条计算如下：
sp0, u = guess(pathxyz, k, s) - 带约束的样条计算如下：
sp, success = spline_neumann(pathxyz, k, s) - 然后打印
success遵循scipy.optimize.minimize标准和基于不等式约束的自定义标准作为不满足约束的点的百分比：
print("success = ", success)
perfo_vs_ineq = (sum([(x < xmin) for x in v[0]]) + sum([(x > xmax) for x in v[0]]) + sum([(y > ymax) for y in v[1]]) )/len(v[0])/2
print("perfo% vs ineq constraints = ", perfo_vs_ineq)
- 约束下的优化执行者：
opt = minimize(err, c0flat, (p, u, t, c_shape, k, w), constraints=con)。它优化了由非约束求解 获得的 - return我们必须重塑
copt = opt.x中的系数才能被splev和copt = opt.x.reshape(c_shape) 使用
splev用于评估两个样条曲线 - 不受约束和受约束 - 这里没有新内容：
c0flat初始化样条的系数
X = np.linspace(0, 1, len(pathxyz)*10)
val0 = splev(X, sp0)
val = splev(X, sp)
scipy.optimize.minimize参数和 return 值在 manual 中进行了解释。所以我在这里只解释具体的内容err是要最小化的 cost。它旨在坚持控制点：
def err(c, p, u, t, c_shape, k, w=None):
"""The error function to minimize"""
diff = (np.array(p) - splev(u, (t, c.reshape(c_shape), k))).flatten()
if w is None:
diff = (diff*diff).sum()
else:
diff = (diff*diff*w).sum() #not sure it is the good way to multiply w
return np.abs(diff)
- 我没有用重量参数进行测试
w。这里需要理解的重要一点是，我们仅使用u提供的曲线坐标对控制点进行评估。成本是控制点与scipy.optimize.minimize 尝试使用新计算系数 - 然后我们得出
constraints=con提供给scipy.optimize.minimize的约束定义为：
c 评估它们的方式之间的差异
con = {'type': 'ineq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'ineq', eqinterv)
#'type': 'eq', 'fun': lambda c: constraint(c, len(p[0]), t, c_shape, k, 'eq', eqinterv)
- 我只使用
ineq特性，因为使用eq特性的测试在我的用例中给出了很差的结果，但如果它可以帮助某人，我已经让代码。因此，不等式和等式约束都是使用函数constraint(c, len(p[0]), t, c_shape, k, 'ineq', eqinterv)计算的。我更喜欢一个函数而不是一个函数列表来只执行一次样条评估。所以当然，c是被scipy.optimize.minimize评估的参数，t和k完成评估所需的(t,c,k)元组，len(p[0])是与要评估的点数成比例，'ineq'告诉constraint处理不等式，eqinterv是一个向量，我想在其中评估计算为成本的等式。在我的用例中，我记得我需要x >= -35 and x <= 2077 and y <= 2802。我没有详细说明特定于我的用例的计算，我只强调间隔与与u: 齐次的曲线坐标相关的点
xmin_umin = xmin_umax = xmax_umin = xmax_umax = ymax_umin = ymax_umax = -1
for i in range(len(p[0])):
if xmin_umin == -1 and p[0][i] <= xmin : xmin_umin = u[i]
if xmin_umin != -1 and xmin_umax == -1 and p[0][i] > xmin : xmin_umax = u[i-1]
if xmax_umin == -1 and p[0][i] >= xmax : xmax_umin = u[i]
if xmax_umin != -1 and xmax_umax == -1 and p[0][i] < xmax : xmax_umax = u[i-1]
if ymax_umin == -1 and p[1][i] >= ymax : ymax_umin = u[i]
if ymax_umin != -1 and ymax_umax == -1 and p[1][i] < ymax : ymax_umax = u[i-1]
- 然后计算等式的成本：
eq_contrib = 0
for i in range(len(X)):
eq_contrib += (X[i] >= eqinterv[0][0] and X[i] <= eqinterv[0][1]) * (v[0][i] - xmin)**2 \
+ (X[i] >= eqinterv[1][0] and X[i] <= eqinterv[1][1]) * (v[0][i] - xmax)**2 \
+ (X[i] >= eqinterv[2][0] and X[i] <= eqinterv[2][1]) * (v[1][i] - ymax)**2
- 不平等成本很简单：
ineq_contrib = sum([(x < xmin)*(x-xmin)**2 + (x > xmax)*(x-xmax)**2 for x in v[0]] \
+ [(y > ymax)*(y-ymax)**2 for y in v[1]])
就这些了，希望有用