如何找到参数化曲线与形状的相交时间？

Question

我有一条按时间参数化的曲线，它与一个形状（在本例中只是一个矩形）相交。在 , I used shapely 之后确定对象相交的位置，但是从那里开始，我很难为何时找到一个好的解决方案。目前，我通过找到最接近交叉点（space）的曲线点，然后使用它的时间戳来粗略地估算时间。

但我相信应该有更好的解决方案，例如通过求解多项式方程，可以使用 numpy 多项式的 root 方法。我只是不确定该怎么做，因为我猜你需要以某种方式引入公差，因为曲线很可能永远不会采用与 shapely 确定的完全相同的交点坐标。

这是我的代码：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Ellipse
from matplotlib.collections import LineCollection
from shapely.geometry import LineString, Polygon


# the parameterized curve
coeffs = np.array([
    [-2.65053088e-05, 2.76890591e-05],
    [-5.70681576e-02, -2.69415587e-01],
    [7.92564148e+02, 6.88557419e+02],
])
t_fit = np.linspace(-2400, 3600, 1000)
x_fit = np.polyval(coeffs[:, 0], t_fit)
y_fit = np.polyval(coeffs[:, 1], t_fit)
curve = LineString(np.column_stack((x_fit, y_fit)))

# the shape it intersects
area = {'x': [700, 1000], 'y': [1300, 1400]}
area_shape = Polygon([
    (area['x'][0], area['y'][0]),
    (area['x'][1], area['y'][0]),
    (area['x'][1], area['y'][1]),
    (area['x'][0], area['y'][1]),
])

# attempt at finding the time of intersection
intersection = curve.intersection(area_shape).coords[-1]
distances = np.hypot(x_fit-intersection[0], y_fit-intersection[1])
idx = np.where(distances == min(distances))
fit_intersection = x_fit[idx][0], y_fit[idx][0]
t_intersection = t_fit[idx]
print(t_intersection)

# code for visualization
fig, ax = plt.subplots(figsize=(5, 5))
ax.margins(0.4, 0.2)
ax.invert_yaxis()

area_artist = Rectangle(
    (area['x'][0], area['y'][0]),
    width=area['x'][1] - area['x'][0],
    height=area['y'][1] - area['y'][0],
    edgecolor='gray', facecolor='none'
)
ax.add_artist(area_artist)

points = np.array([x_fit, y_fit]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
z = np.linspace(0, 1, points.shape[0])
norm = plt.Normalize(z.min(), z.max())
lc = LineCollection(
    segments, cmap='autumn', norm=norm, alpha=1,
    linewidths=2, picker=8, capstyle='round',
    joinstyle='round'
)
lc.set_array(z)
ax.add_collection(lc)

ax.autoscale_view()
ax.relim()

trans = (ax.transData + ax.transAxes.inverted()).transform
intersection_point = Ellipse(
    xy=trans(fit_intersection), width=0.02, height=0.02, fc='none',
    ec='black', transform=ax.transAxes, zorder=3,
)
ax.add_artist(intersection_point)

plt.show()

为了视觉效果，这里是问题在情节中的样子：

Answer 1

最好的方法是使用插值函数来计算(x(t), y(t))。并使用一个函数来计算 d(t)：到交叉点的距离。然后我们在 d(t) 上使用 scipy.optimize.minimize 来找到 d(t) 最小的 t 值。插值将确保良好的准确性。

因此，我对您的代码进行了一些修改。

插值函数定义和距离计算
测试是否确实有交集，否则没有意义
通过最小化计算相交时间

代码（更新）：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Ellipse
from matplotlib.collections import LineCollection
from shapely.geometry import LineString, Polygon
from scipy.optimize import minimize

# Interpolate (x,y) at time t:
def interp_xy(t,tp, fpx,fpy):
        # tp: time grid points, fpx, fpy: the corresponding x,y values
        x=np.interp(t, tp, fpx)
        y=np.interp(t, tp, fpy)
        return x,y

# Compute distance to intersection:
def dist_to_intersect(t,tp, fpx, fpy, intersection):
        x,y = interp_xy(t,tp,fpx,fpy)
        d=np.hypot(x-intersection[0], y-intersection[1])
        return d



# the parameterized curve
t_fit = np.linspace(-2400, 3600, 1000)
#t_fit = np.linspace(-4200, 0, 1000)
coeffs = np.array([[-2.65053088e-05, 2.76890591e-05],[-5.70681576e-02, -2.69415587e-01],[7.92564148e+02, 6.88557419e+02],])

#t_fit = np.linspace(-2400, 3600, 1000)
#coeffs = np.array([[4.90972365e-05, -2.03897149e-04],[2.19222264e-01, -1.63335372e+00],[9.33624672e+02,  1.07067102e+03], ])

#t_fit = np.linspace(-2400, 3600, 1000)
#coeffs = np.array([[-2.63100091e-05, -7.16542227e-05],[-5.60829940e-04, -3.19183803e-01],[7.01544289e+02,  1.24732452e+03], ])

#t_fit = np.linspace(-2400, 3600, 1000)
#coeffs = np.array([[-2.63574223e-05, -9.15525038e-05],[-8.91039302e-02, -4.13843734e-01],[6.35650643e+02,  9.40010900e+02], ])

x_fit = np.polyval(coeffs[:, 0], t_fit)
y_fit = np.polyval(coeffs[:, 1], t_fit)
curve = LineString(np.column_stack((x_fit, y_fit)))


# the shape it intersects
area = {'x': [700, 1000], 'y': [1300, 1400]}
area_shape = Polygon([
    (area['x'][0], area['y'][0]),
    (area['x'][1], area['y'][0]),
    (area['x'][1], area['y'][1]),
    (area['x'][0], area['y'][1]),
])

# attempt at finding the time of intersection
curve_intersection = curve.intersection(area_shape)

# We check if intersection is empty or not:
if not curve_intersection.is_empty:   
    # We can get the coords because intersection is not empty
    intersection=curve_intersection.coords[-1]
    
    distances = np.hypot(x_fit-intersection[0], y_fit-intersection[1])

    print("Looking for minimal distance to intersection: ")
    print('-------------------------------------------------------------------------')
    # Call to minimize:
    # We pass:
    # - the function to be minimized (dist_to_intersect)
    # - a starting value to t 
    # - arguments, method and tolerance tol. The minimization will succeed when 
    #   dist_to_intersect <  tol=1e-6
    # - option: here -->  verbose
    dmin=np.min((x_fit-intersection[0])**2+(y_fit-intersection[1])**2)
    index=np.where((x_fit-intersection[0])**2+(y_fit-intersection[1])**2==dmin)
    t0=t_fit[index]
    res = minimize(dist_to_intersect, t0,  args=(t_fit, x_fit, y_fit, intersection), method='Nelder-Mead',tol = 1e-6, options={ 'disp': True})
    print('-------------------------------------------------------------------------')
    print("Result of the optimization:")
    print(res)
    print('-------------------------------------------------------------------------')
    print("Intersection at time t = ",res.x[0])    
    fit_intersection = interp_xy(res.x[0],t_fit, x_fit,y_fit)
    print("Intersection point : ",fit_intersection)
else:
    print("No intersection.")
    
    
# code for visualization
fig, ax = plt.subplots(figsize=(5, 5))
ax.margins(0.4, 0.2)
ax.invert_yaxis()

area_artist = Rectangle(
    (area['x'][0], area['y'][0]),
    width=area['x'][1] - area['x'][0],
    height=area['y'][1] - area['y'][0],
    edgecolor='gray', facecolor='none'
)
ax.add_artist(area_artist)

#plt.plot(x_fit,y_fit)
points = np.array([x_fit, y_fit]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
z = np.linspace(0, 1, points.shape[0])
norm = plt.Normalize(z.min(), z.max())
lc = LineCollection(
    segments, cmap='autumn', norm=norm, alpha=1,
    linewidths=2, picker=8, capstyle='round',
    joinstyle='round'
)
lc.set_array(z)
ax.add_collection(lc)
# Again, we check that intersection exists because we don't want to draw
# an non-existing point (it would generate an error)

if not curve_intersection.is_empty:
    plt.plot(fit_intersection[0],fit_intersection[1],'o')


plt.show()

输出：

Looking for minimal distance to intersection: 
-------------------------------------------------------------------------
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 31
         Function evaluations: 62
-------------------------------------------------------------------------
Result of the optimization:
 final_simplex: (array([[-1898.91943932],
       [-1898.91944021]]), array([8.44804735e-09, 3.28684898e-07]))
           fun: 8.448047349426054e-09
       message: 'Optimization terminated successfully.'
          nfev: 62
           nit: 31
        status: 0
       success: True
             x: array([-1898.91943932])
-------------------------------------------------------------------------
Intersection at time t =  -1898.919439315796
Intersection point :  (805.3563860471179, 1299.9999999916085)

而您的代码给出的结果不那么精确：

t=-1901.5015015 

intersection point: (805.2438793482748,1300.9671136070717)

图：

如何找到参数化曲线与形状的相交时间？

How to find the intersection time of a parameterized curve with a shape?

python

geometry

numpy

parameterization

shapely