如何用geomdl或其他库为花键添加边界约束？
问题描述 投票：0回答：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 库没有提出带约束的花键。我试过这个黑客，只是通过修正点来保持在边界内。
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]
如何使它保持平稳，在路径上，并尊重约束条件，请问，可能用另一个库？我发现 这个但这解决了导数约束，我不知道如何调整这个解决方案。我也从严格的数学角度提出了这个问题。此处.
1个回答
0
投票
投票
好吧，题目很难，但我得到了它，灵感来自于 这个 二维边界（导数）约束花键。所提出的解决方案还利用了 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的渲染。三维视图显示，花键使用z轴来平滑。这对我的用例来说并不满意，所以我必须在我的约束条件中考虑到这一点，但这不在本次QA的范围内。
二维视图显示了约束条件对花键的影响 And the 2D view that shows the constraints effects on the spline:
蓝色的曲线是没有约束的，红色的是有约束的。
现在解释一下走向。
- 没有约束条件的花键是用..:
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的非约束性解法得到的 - 它返回的系数在
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参数和返回值在 手册. 所以我在这里只解释一下具体的内容。 err是 成本 以减少。它的建立是为了坚持控制点。
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. 成本是控制点之间的差异，以及如何用新计算的系数来评估它们c受scipy.optimize.minimize - 然后，我们再来看看提供给
scipy.optimize.minimize由constraints=con定义为。
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测试以来的实际情况equalities在我的使用案例中给出了很差的结果，但我已经让代码，如果它可能会帮助别人。所以，不等式和等式约束都是用函数来计算的。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]])
就是这样，希望有用。
最新问题
- 使用接受闭包的方法在 Rust 中创建对象安全特征
- PactSwift 如何在 CI 上构建
- 当我部署在服务器中时，Angular Service Worker 离线缓存不起作用
- 通过MySQLdb创建函数
- 我应该删除元标记中的尾部斜杠吗？
- 将任何深度转换为 CoreML
- CDK 拖放预览位于模式对话框后面
- 使用 Sharpie 生成 ApiDefinition.cs 时，使用 Swift 协议构建的 Xamarin iOS 绑定库失败
- 如何在React js中修复移动设备上的网站高度和宽度？
- 从交互式网格中获取所选行的 ID 并将值设置为另一个页面的项目
- Terraform - 如何获取 azurerm 密钥保管库访问策略的应用服务对象 ID？
- 意外激活 UFW 防火墙后如何恢复对亚马逊 EC2 实例的 SSH 访问？
- 用于屏蔽信用卡和帐号的正则表达式[已关闭]
- 如何计算java卡小程序CAP文件的NVM大小
- 如何禁用 TestNG 控制台日志记录 [main] INFO org.testng.internal.Utils -
- Kafka Broker 身份验证失败 - 凭据无效
- Blazor、会话存储和弃用参考库
- React + Bootstrap 5.3：如何在 Bootstrap 的 SASS 变量上应用条件 CSS 选择器覆盖？
- 使用 Business Central Dynamics AL 获取营销文本
- Nmap 主机发现
© www.soinside.com 2019 - 2024. All rights reserved.