numpy 中的曲线曲率
Curve curvature in numpy
我正在使用特殊相机以 1s 的固定时间间隔测量物体的 x、y 坐标（以厘米为单位）。我有一个 numpy 数组中的数据：
a = np.array([ [ 0. , 0. ],[ 0.3 , 0. ],[ 1.25, -0.1 ],[ 2.1 , -0.9 ],[ 2.85, -2.3 ],[ 3.8 , -3.95],[ 5. , -5.75],[ 6.4 , -7.8 ],[ 8.05, -9.9 ],[ 9.9 , -11.6 ],[ 12.05, -12.85],[ 14.25, -13.7 ],[ 16.5 , -13.8 ],[ 19.25, -13.35],[ 21.3 , -12.2 ],[ 22.8 , -10.5 ],[ 23.55, -8.15],[ 22.95, -6.1 ],[ 21.35, -3.95],[ 19.1 , -1.9 ]])
曲线看起来像这样：
plt.scatter(a[:,0], a[:,1])
问题：
如何计算每个点的切向和径向加速度矢量？我发现了一些可能相关的公式：
我可以很容易地用 np.diff(a, axis=0) 计算 vx 和 vy 预测，但我是一个 numpy/python 菜鸟，我很难继续下去.如果我能计算出每个点的曲率，我的问题也就解决了。有人可以帮忙吗？
编辑：我花了几个小时断断续续地整理了这个答案，所以我错过了你最新的编辑，表明你只需要曲率。希望这个答案无论如何都会有所帮助。
除了进行一些曲线拟合外，我们近似导数的方法是通过 finite differences. Thankfully, numpy has a gradient 方法为我们进行这些差异计算，处理每个内部点的平均前一个和下一个斜率的细节，保留每个端点，等等
import numpy as np
a = np.array([ [ 0. , 0. ],[ 0.3 , 0. ],[ 1.25, -0.1 ],
[ 2.1 , -0.9 ],[ 2.85, -2.3 ],[ 3.8 , -3.95],
[ 5. , -5.75],[ 6.4 , -7.8 ],[ 8.05, -9.9 ],
[ 9.9 , -11.6 ],[ 12.05, -12.85],[ 14.25, -13.7 ],
[ 16.5 , -13.8 ],[ 19.25, -13.35],[ 21.3 , -12.2 ],
[ 22.8 , -10.5 ],[ 23.55, -8.15],[ 22.95, -6.1 ],
[ 21.35, -3.95],[ 19.1 , -1.9 ]])
现在，我们计算每个变量的导数并将它们放在一起（出于某种原因，如果我们只调用 np.gradient(a)，我们会得到一个数组列表...不确定我的头脑那里发生了什么，但我现在就解决它):
dx_dt = np.gradient(a[:, 0])
dy_dt = np.gradient(a[:, 1])
velocity = np.array([ [dx_dt[i], dy_dt[i]] for i in range(dx_dt.size)])
这为我们提供了 velocity 的以下向量：
array([[ 0.3 , 0. ],
[ 0.625, -0.05 ],
[ 0.9 , -0.45 ],
[ 0.8 , -1.1 ],
[ 0.85 , -1.525],
[ 1.075, -1.725],
[ 1.3 , -1.925],
[ 1.525, -2.075],
[ 1.75 , -1.9 ],
[ 2. , -1.475],
[ 2.175, -1.05 ],
[ 2.225, -0.475],
[ 2.5 , 0.175],
[ 2.4 , 0.8 ],
[ 1.775, 1.425],
[ 1.125, 2.025],
[ 0.075, 2.2 ],
[-1.1 , 2.1 ],
[-1.925, 2.1 ],
[-2.25 , 2.05 ]])
看一下 a.
的散点图就明白了
现在，对于速度，我们取速度向量的长度。然而，有一件事我们在这里没有真正记住：一切都是t的函数。因此，ds/dt 实际上是 t 的标量函数（与 t 的矢量函数相反），就像 dx/dt 和 dy/dt 一样。因此，我们将 ds_dt 表示为每个一秒时间间隔的 numpy 值数组，每个值对应于每秒速度的近似值：
ds_dt = np.sqrt(dx_dt * dx_dt + dy_dt * dy_dt)
这会产生以下数组：
array([ 0.3 , 0.62699681, 1.00623059, 1.36014705, 1.74588803,
2.03254766, 2.32284847, 2.57512136, 2.58311827, 2.48508048,
2.41518633, 2.27513736, 2.50611752, 2.52982213, 2.27623593,
2.31651678, 2.20127804, 2.37065392, 2.8487936 , 3.04384625])
当您查看 a 的散点图上的点之间的间隙时，这又是有一定意义的：物体加快速度，在转弯时稍微减速，然后再加速备份更多。
现在，为了找到单位切线向量，我们需要对 ds_dt 做一个小的变换，使其大小与 velocity 的大小相同（这有效地允许我们将向量值函数 velocity 除以标量函数 ds_dt 的（表示）：
tangent = np.array([1/ds_dt] * 2).transpose() * velocity
这会产生以下 numpy 数组：
array([[ 1. , 0. ],
[ 0.99681528, -0.07974522],
[ 0.89442719, -0.4472136 ],
[ 0.5881717 , -0.80873608],
[ 0.48685826, -0.87348099],
[ 0.52889289, -0.84868859],
[ 0.55965769, -0.82872388],
[ 0.5922051 , -0.80578727],
[ 0.67747575, -0.73554511],
[ 0.80480291, -0.59354215],
[ 0.90055164, -0.43474907],
[ 0.97796293, -0.2087786 ],
[ 0.99755897, 0.06982913],
[ 0.9486833 , 0.31622777],
[ 0.77979614, 0.62603352],
[ 0.48564293, 0.87415728],
[ 0.03407112, 0.99941941],
[-0.46400699, 0.88583154],
[-0.67572463, 0.73715414],
[-0.73919634, 0.67349 ]])
注意两件事：1. 在 t 的每个值处，tangent 指向与 velocity 相同的方向，以及 2. 在 [=26= 的每个值处]，tangent是单位向量。确实：
在[12]中：
In [12]: np.sqrt(tangent[:,0] * tangent[:,0] + tangent[:,1] * tangent[:,1])
Out[12]:
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1.])
现在，由于我们取切向量的导数除以它的长度来得到单位法向量，所以我们用同样的技巧（为了方便，将 tangent 的分量分开）：
tangent_x = tangent[:, 0]
tangent_y = tangent[:, 1]
deriv_tangent_x = np.gradient(tangent_x)
deriv_tangent_y = np.gradient(tangent_y)
dT_dt = np.array([ [deriv_tangent_x[i], deriv_tangent_y[i]] for i in range(deriv_tangent_x.size)])
length_dT_dt = np.sqrt(deriv_tangent_x * deriv_tangent_x + deriv_tangent_y * deriv_tangent_y)
normal = np.array([1/length_dT_dt] * 2).transpose() * dT_dt
这为我们提供了 normal 的以下向量：
array([[-0.03990439, -0.9992035 ],
[-0.22975292, -0.97324899],
[-0.48897562, -0.87229745],
[-0.69107645, -0.72278167],
[-0.8292422 , -0.55888941],
[ 0.85188045, 0.52373629],
[ 0.8278434 , 0.56095927],
[ 0.78434982, 0.62031876],
[ 0.70769355, 0.70651953],
[ 0.59568265, 0.80321988],
[ 0.41039706, 0.91190693],
[ 0.18879684, 0.98201617],
[-0.05568352, 0.99844847],
[-0.36457012, 0.93117594],
[-0.63863584, 0.76950911],
[-0.89417603, 0.44771557],
[-0.99992445, 0.0122923 ],
[-0.93801622, -0.34659137],
[-0.79170904, -0.61089835],
[-0.70603568, -0.70817626]])
请注意，法向量表示曲线转向的方向。当结合 a 的散点图查看时，上面的矢量就有意义了。特别是，我们在第 5 个点后从向下转向向上，在第 12 个点后我们开始向左（相对于 x 轴）转向。
最后，为了得到加速度的切向分量和法向分量，我们需要 s、x 和 y 关于 t 的二阶导数，然后我们可以获得曲率和其余分量（请记住它们都是 t 的标量函数）：
d2s_dt2 = np.gradient(ds_dt)
d2x_dt2 = np.gradient(dx_dt)
d2y_dt2 = np.gradient(dy_dt)
curvature = np.abs(d2x_dt2 * dy_dt - dx_dt * d2y_dt2) / (dx_dt * dx_dt + dy_dt * dy_dt)**1.5
t_component = np.array([d2s_dt2] * 2).transpose()
n_component = np.array([curvature * ds_dt * ds_dt] * 2).transpose()
acceleration = t_component * tangent + n_component * normal
我正在使用特殊相机以 1s 的固定时间间隔测量物体的 x、y 坐标（以厘米为单位）。我有一个 numpy 数组中的数据：
a = np.array([ [ 0. , 0. ],[ 0.3 , 0. ],[ 1.25, -0.1 ],[ 2.1 , -0.9 ],[ 2.85, -2.3 ],[ 3.8 , -3.95],[ 5. , -5.75],[ 6.4 , -7.8 ],[ 8.05, -9.9 ],[ 9.9 , -11.6 ],[ 12.05, -12.85],[ 14.25, -13.7 ],[ 16.5 , -13.8 ],[ 19.25, -13.35],[ 21.3 , -12.2 ],[ 22.8 , -10.5 ],[ 23.55, -8.15],[ 22.95, -6.1 ],[ 21.35, -3.95],[ 19.1 , -1.9 ]])
曲线看起来像这样：
plt.scatter(a[:,0], a[:,1])
问题：
如何计算每个点的切向和径向加速度矢量？我发现了一些可能相关的公式：
我可以很容易地用 np.diff(a, axis=0) 计算 vx 和 vy 预测，但我是一个 numpy/python 菜鸟，我很难继续下去.如果我能计算出每个点的曲率，我的问题也就解决了。有人可以帮忙吗？
编辑：我花了几个小时断断续续地整理了这个答案，所以我错过了你最新的编辑，表明你只需要曲率。希望这个答案无论如何都会有所帮助。
除了进行一些曲线拟合外，我们近似导数的方法是通过 finite differences. Thankfully, numpy has a gradient 方法为我们进行这些差异计算，处理每个内部点的平均前一个和下一个斜率的细节，保留每个端点，等等
import numpy as np
a = np.array([ [ 0. , 0. ],[ 0.3 , 0. ],[ 1.25, -0.1 ],
[ 2.1 , -0.9 ],[ 2.85, -2.3 ],[ 3.8 , -3.95],
[ 5. , -5.75],[ 6.4 , -7.8 ],[ 8.05, -9.9 ],
[ 9.9 , -11.6 ],[ 12.05, -12.85],[ 14.25, -13.7 ],
[ 16.5 , -13.8 ],[ 19.25, -13.35],[ 21.3 , -12.2 ],
[ 22.8 , -10.5 ],[ 23.55, -8.15],[ 22.95, -6.1 ],
[ 21.35, -3.95],[ 19.1 , -1.9 ]])
现在，我们计算每个变量的导数并将它们放在一起（出于某种原因，如果我们只调用 np.gradient(a)，我们会得到一个数组列表...不确定我的头脑那里发生了什么，但我现在就解决它):
dx_dt = np.gradient(a[:, 0])
dy_dt = np.gradient(a[:, 1])
velocity = np.array([ [dx_dt[i], dy_dt[i]] for i in range(dx_dt.size)])
这为我们提供了 velocity 的以下向量：
array([[ 0.3 , 0. ],
[ 0.625, -0.05 ],
[ 0.9 , -0.45 ],
[ 0.8 , -1.1 ],
[ 0.85 , -1.525],
[ 1.075, -1.725],
[ 1.3 , -1.925],
[ 1.525, -2.075],
[ 1.75 , -1.9 ],
[ 2. , -1.475],
[ 2.175, -1.05 ],
[ 2.225, -0.475],
[ 2.5 , 0.175],
[ 2.4 , 0.8 ],
[ 1.775, 1.425],
[ 1.125, 2.025],
[ 0.075, 2.2 ],
[-1.1 , 2.1 ],
[-1.925, 2.1 ],
[-2.25 , 2.05 ]])
看一下 a.
现在，对于速度，我们取速度向量的长度。然而，有一件事我们在这里没有真正记住：一切都是t的函数。因此，ds/dt 实际上是 t 的标量函数（与 t 的矢量函数相反），就像 dx/dt 和 dy/dt 一样。因此，我们将 ds_dt 表示为每个一秒时间间隔的 numpy 值数组，每个值对应于每秒速度的近似值：
ds_dt = np.sqrt(dx_dt * dx_dt + dy_dt * dy_dt)
这会产生以下数组：
array([ 0.3 , 0.62699681, 1.00623059, 1.36014705, 1.74588803,
2.03254766, 2.32284847, 2.57512136, 2.58311827, 2.48508048,
2.41518633, 2.27513736, 2.50611752, 2.52982213, 2.27623593,
2.31651678, 2.20127804, 2.37065392, 2.8487936 , 3.04384625])
当您查看 a 的散点图上的点之间的间隙时，这又是有一定意义的：物体加快速度，在转弯时稍微减速，然后再加速备份更多。
现在，为了找到单位切线向量，我们需要对 ds_dt 做一个小的变换，使其大小与 velocity 的大小相同（这有效地允许我们将向量值函数 velocity 除以标量函数 ds_dt 的（表示）：
tangent = np.array([1/ds_dt] * 2).transpose() * velocity
这会产生以下 numpy 数组：
array([[ 1. , 0. ],
[ 0.99681528, -0.07974522],
[ 0.89442719, -0.4472136 ],
[ 0.5881717 , -0.80873608],
[ 0.48685826, -0.87348099],
[ 0.52889289, -0.84868859],
[ 0.55965769, -0.82872388],
[ 0.5922051 , -0.80578727],
[ 0.67747575, -0.73554511],
[ 0.80480291, -0.59354215],
[ 0.90055164, -0.43474907],
[ 0.97796293, -0.2087786 ],
[ 0.99755897, 0.06982913],
[ 0.9486833 , 0.31622777],
[ 0.77979614, 0.62603352],
[ 0.48564293, 0.87415728],
[ 0.03407112, 0.99941941],
[-0.46400699, 0.88583154],
[-0.67572463, 0.73715414],
[-0.73919634, 0.67349 ]])
注意两件事：1. 在 t 的每个值处，tangent 指向与 velocity 相同的方向，以及 2. 在 [=26= 的每个值处]，tangent是单位向量。确实：
在[12]中：
In [12]: np.sqrt(tangent[:,0] * tangent[:,0] + tangent[:,1] * tangent[:,1])
Out[12]:
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1.])
现在，由于我们取切向量的导数除以它的长度来得到单位法向量，所以我们用同样的技巧（为了方便，将 tangent 的分量分开）：
tangent_x = tangent[:, 0]
tangent_y = tangent[:, 1]
deriv_tangent_x = np.gradient(tangent_x)
deriv_tangent_y = np.gradient(tangent_y)
dT_dt = np.array([ [deriv_tangent_x[i], deriv_tangent_y[i]] for i in range(deriv_tangent_x.size)])
length_dT_dt = np.sqrt(deriv_tangent_x * deriv_tangent_x + deriv_tangent_y * deriv_tangent_y)
normal = np.array([1/length_dT_dt] * 2).transpose() * dT_dt
这为我们提供了 normal 的以下向量：
array([[-0.03990439, -0.9992035 ],
[-0.22975292, -0.97324899],
[-0.48897562, -0.87229745],
[-0.69107645, -0.72278167],
[-0.8292422 , -0.55888941],
[ 0.85188045, 0.52373629],
[ 0.8278434 , 0.56095927],
[ 0.78434982, 0.62031876],
[ 0.70769355, 0.70651953],
[ 0.59568265, 0.80321988],
[ 0.41039706, 0.91190693],
[ 0.18879684, 0.98201617],
[-0.05568352, 0.99844847],
[-0.36457012, 0.93117594],
[-0.63863584, 0.76950911],
[-0.89417603, 0.44771557],
[-0.99992445, 0.0122923 ],
[-0.93801622, -0.34659137],
[-0.79170904, -0.61089835],
[-0.70603568, -0.70817626]])
请注意，法向量表示曲线转向的方向。当结合 a 的散点图查看时，上面的矢量就有意义了。特别是，我们在第 5 个点后从向下转向向上，在第 12 个点后我们开始向左（相对于 x 轴）转向。
最后，为了得到加速度的切向分量和法向分量，我们需要 s、x 和 y 关于 t 的二阶导数，然后我们可以获得曲率和其余分量（请记住它们都是 t 的标量函数）：
d2s_dt2 = np.gradient(ds_dt)
d2x_dt2 = np.gradient(dx_dt)
d2y_dt2 = np.gradient(dy_dt)
curvature = np.abs(d2x_dt2 * dy_dt - dx_dt * d2y_dt2) / (dx_dt * dx_dt + dy_dt * dy_dt)**1.5
t_component = np.array([d2s_dt2] * 2).transpose()
n_component = np.array([curvature * ds_dt * ds_dt] * 2).transpose()
acceleration = t_component * tangent + n_component * normal