我正在寻找一种对样条曲线上的点进行采样的快速方法，这样通过这些点的多边形或线串不会超过原始样条曲线的给定弦误差。我有一个我前段时间写的算法，它产生了图片中的结果（如果有兴趣，请参阅下面的代码；我不希望任何人仔细研究它）。它工作正常，但速度不是很快（在我的计算机上生成该图形大约需要 2 秒）。有没有更简单的东西，可能内置到 numpy 或 scipy 中来完成这个？
谢谢
import numpy as np
from scipy.optimize import brentq
from scipy.interpolate import splev
def get_rhos(ts, tck):
"""Get (signed) rhos (1/rad of curvature) for a given
set of t values.
"""
tanvs = np.array(splev(ts, tck, der=1)).T
accvs = np.array(splev(ts, tck, der=2)).T
if tanvs.ndim == 1:
tanvs = tanvs.reshape(1, -1)
accvs = accvs.reshape(1, -1)
crossp = np.cross(accvs, tanvs, axis=1)
tanvms = np.array([np.sqrt(np.dot(v, v)) for v in tanvs])
rhos = crossp / tanvms**3
return rhos
def calc_rad(pt0, pt1, pt2, calcdrop=False):
"""Calculate a radius from three points on the arc.
Lifted from http://www.physicsforums.com/showthread.php?t=173847
"""
pt0 = np.array(pt0)
pt1 = np.array(pt1)
pt2 = np.array(pt2)
v0 = pt1 - pt0
v1 = pt2 - pt0
v2 = pt2 - pt1
a = np.sqrt(np.dot(v0, v0))
b = np.sqrt(np.dot(v1, v1))
c = np.sqrt(np.dot(v2, v2))
R = (a*b*c) / np.sqrt( 2 * a**2 * b**2
+ 2 * b**2 * c**2
+ 2 * c**2 * a**2
- a**4 - b**4 - c**4)
if calcdrop:
# Calculate arc drop
drop = R - np.sqrt(R**2 - (b/2.)**2)
return R, drop
else:
return R
def chordal_sample(tck, chordaltol, oversample=10):
"""Given a spline definition and a chordal tolerance
(intol/outol), get the t-values for the spline such
that, when adjacent points are connected, the chordal
tolerance is not violated.
Accomplishes this by bracketing a solution, then using
the brentq solver to find the point where the chordal
error equals the chordal tolerance.
Note that a few extra points may be inserted where there
are inflections in the cubic; these are sometimes missed
by the arc-radius-calculating portion of the code.
"""
# This is the function we'll need when we have to
# go searching for the answer via brentq
def makeerrfunc(st, spt, tck, chordaltol):
def errfunc(et):
mt = (st + et) / 2.0
mpt = np.array(splev(mt, tck))
ept = np.array(splev(et, tck))
_, arcdrop = calc_rad(spt, mpt, ept, calcdrop=True)
diff = arcdrop - chordaltol
return diff
return errfunc
# Make sure we're sampling enough points
# TODO: How can we be sure?
ts = np.linspace(0, 1, oversample * len(tck[1][0]))
newts = [0]
# Loop through the time values
for nt in ts:
st = newts[-1]
rts = ts[ts > st] # Only consider remaining time values
# Step through adjacent pairs of time values and find
# ones that bracket the solution.
for et0, et1 in zip(rts[0:-1], rts[1:]):
# Get a 'middle time' that we can use to calc
# a 'middle point' for our arc calculations
mt0 = (st + et0) / 2.
mt1 = (st + et1) / 2.
# Interpolate points at the critical t values
ipts = np.array(splev([st, mt0, et0, mt1, et1], tck))
spt, mpt0, ept0, mpt1, ept1 = ipts.T
_, arcdrop0 = calc_rad(spt, mpt0, ept0, calcdrop=True)
_, arcdrop1 = calc_rad(spt, mpt1, ept1, calcdrop=True)
# Have we bracketed the solution yet? If so, use
# brentq to find a better one within the bracketed
# range, then move on to a new start t.
if arcdrop0 > chordaltol: # Check the initial pair
errfunc = makeerrfunc(st, spt, tck, chordaltol)
mdt = brentq(errfunc, st, et0)
newts.append(mdt)
break
if arcdrop0 <= chordaltol and arcdrop1 > chordaltol:
errfunc = makeerrfunc(st, spt, tck, chordaltol)
mdt = brentq(errfunc, et0, et1)
newts.append(mdt)
break
# Check for the existence of an inflection point
# in the bracketed range by checking the signs
# of the two calculated curvatures and looking for
# a reversal.
if get_rhos(et0, tck)[0] * get_rhos(et1, tck)[0] < 0:
newts.append((et0 + et1) / 2.0)
break
if et1 == 1.0: # No more points to try
newts.append(1.0)
break
return newts
if __name__ == '__main__':
import matplotlib.pyplot as plt
from scipy.interpolate import splprep
# Create a hi-res sample spline. Start with some
# low-res points and then resample at a higher
# res.
XY = np.array([[0.0, 1.0, 2.0, 3.0, 2.0, 1.0, 0.0],
[0.0, -1.0, -0.5, 0.0, 2.5, 1.2, 2.0]])
tck, u = splprep(XY, s=0)
XY = splev(np.linspace(0, 1, 400), tck)
tck, u = splprep(XY, s=0)
# Get a set of t values that will plot out
# a linestring with no more than 0.1 chordal
# error to the original.
ts = chordal_sample(tck, 0.1)
fig, ax = plt.subplots()
# Plot the hi-res spline
ax.plot(*XY)
# Plot the approximated spline
ax.plot(*np.array(splev(ts, tck)), marker='o')
ax.axis('equal')
ax.grid()
plt.show()
