numpy中的曲线曲率答案

【问题标题】：Curve curvature in numpynumpy中的曲线曲率
【发布时间】：2015-02-01 23:55:51
【问题描述】：

我正在使用特殊相机以 1 秒的固定时间间隔测量对象的 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 菜鸟，我想继续。如果我可以计算每个点的曲率，我的问题也将得到解决。有人可以帮忙吗？

【问题讨论】：

好吧，除非你有一个很好的曲线参数方程系统（我不相信你这样做），否则你必须用(Delta x)/(Delta t)替换，例如x'（原谅粗略的数学符号，因为 SO 不支持 LaTeX）。因为你的时间间隔都是一秒，Delta t 是 1，所以你可以用Delta x 替换x'，同样用y'。当然，这些将是近似值，但除非您尝试进行一些曲线拟合，否则这是您可以使用现有数据集做的最好的事情。
糟糕，抱歉，我的回答有点啰嗦，因为您只是在曲率之后。希望它仍然有用。

标签： python numpy

【解决方案1】：

编辑：我在几个小时内断断续续地整理了这个答案，所以我错过了你最新的编辑，表明你只需要曲率。希望这个答案无论如何都会有所帮助。

除了做一些曲线拟合，我们逼近导数的方法是通过finite differences。值得庆幸的是，numpy 有一个 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. 在t 的每个值处，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 的散点图一起查看时，上面的向量是有意义的。特别是，我们在第五个点之后从向下转向向上，在第 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

【讨论】：

非常感谢@JackManey 帮助社区！
嗨@JackManey，很好的答案。我有一个疑问，负值对法线向量表示什么。我应该将它们添加到“a”中以获得法线到曲线点的实际值吗？
@user3515225：谢谢。至于你的问题......好吧，法线向量是一个向量。允许有负坐标。这不是问题。
这个计算是不是有点矫枉过正？据我了解，为了获得曲率，计算表达式 length_dT_dt / ds_dt 就足够了。
这是一篇了不起的帖子！非常感谢！唯一让我感到困惑的是“t_component”和“n_component”位。当我绘制曲率时，我得到了正确的东西（我猜）。 “t_component”、“n_component”和“acceleration”位是干什么用的？谢谢。