如何使用 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

我要添加约束:

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 中的系数才能被 splevcopt = 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评估的参数,tk完成评估所需的(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]])

就这些了,希望有用