numpy曲线曲率 [英] Curve curvature in numpy
问题描述
我正在使用专用摄像机以固定的1s时间间隔测量对象的x，y坐标(以厘米为单位).我有一个numpy数组中的数据:
I am measuring x,y coordinates (in cm) of an object with a special camera in fixed time intervals of 1s. I have the data in a numpy array:
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])
问题:
如何计算每个点的切向和径向加速度矢量?我发现了一些可能相关的公式:
How can I calculate the tangential and the radial aceleration vectors at each point? I found some formulas that might be relevant:
我可以使用np.diff(a, axis=0)轻松计算vx和vy投影，但是我是一个numpy/python noob，这是我继续研究的方法.如果我可以计算每个点的曲率，那么我的问题也将得到解决.有人可以帮忙吗?
I am able to easily calculate the vx and the vy projections with np.diff(a, axis=0) but I am a numpy/python noob and it is way over my head to continue. If I could calculate the curvature at each point, also my problem would be solved. Can somebody help?
推荐答案
编辑:我在几个小时内反复整理此答案，因此我错过了您的最新编辑，这表明您只所需的曲率.希望这个答案对您有所帮助.
EDIT: I put together this answer off and on over a couple of hours, so I missed your latest edits indicating that you only needed curvature. Hopefully, this answer will be helpful regardless.
除了进行一些曲线拟合外，我们的导数近似方法是通过有限差分.幸运的是，numpy有一个 gradient 方法，可以执行这些差异计算对我们来说，要注意对每个内部点的上一个和下一个斜率求平均，而使每个端点不受影响等细节.
Other than doing some curve-fitting, our method of approximating derivatives is via finite differences. Thankfully, numpy has a gradient method that does these difference calculations for us, taking care of the details of averaging previous and next slopes for each interior point and leaving each endpoint alone, etc.
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)，我们将得到一个数组列表...不确定我的头绪是什么) ，但我现在就解决它):
Now, we compute the derivatives of each variable and put them together (for some reason, if we just call np.gradient(a), we get a list of arrays...not sure off the top of my head what's going on there, but I'll just work around it for now):
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的以下向量:
This gives us the following vector for 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的散点图时很有意义.
which makes sense when glancing at the scatterplot of a.
现在，对于速度，我们采用速度矢量的长度.但是，这里我们没有真正记住的一件事:一切都是t 的功能.因此，ds/dt实际上是t的标量函数(与t的向量函数相反)，就像dx/dt和dy/dt一样.因此，我们将ds_dt表示为一秒钟时间间隔中每个值的numpy数组，每个值对应于每秒速度的近似值:
Now, for speed, we take the length of the velocity vector. However, there's one thing that we haven't really kept in mind here: everything is a function of t. Thus, ds/dt is really a scalar function of t (as opposed to a vector function of t), just like dx/dt and dy/dt. Thus, we will represent ds_dt as a numpy array of values at each of the one second time intervals, each value corresponding to an approximation of the speed at each second:
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散点图上的点之间的间隙时，这再次有意义:对象加快了速度，在拐角时放慢了一点，然后再加快了速度
which, again, makes some sense as you look at the gaps between the dots on the scatterplot of a: the object picks up speed, slowing down a bit as it takes the corner, and then speeds back up even more.
现在，为了找到单位正切向量，我们需要对ds_dt进行较小的转换，以使其大小与velocity相同(这实际上使我们可以对向量值函数进行除法) velocity由标量函数ds_dt的(表示):
Now, in order to find the unit tangent vector, we need to make a small transformation to ds_dt so that its size is the same as that of velocity (this effectively allows us to divide the vector-valued function velocity by the (representation of) the scalar function 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是单位矢量.确实:
Note two things: 1. At each value of t, tangent is pointing in the same direction as velocity, and 2. At each value of t, tangent is a unit vector. Indeed:
在[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的分量隔离开):
Now, since we take the derivative of the tangent vector and divide by its length to get the unit normal vector, we do the same trick (isolating the components of tangent for convenience):
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的以下向量:
This gives us the following vector for 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的散点图一起查看时，上面的矢量才有意义.特别是，我们从第五个点开始向下调高，然后在第十二个点之后开始向左(相对于x轴)转向.
Note that the normal vector represents the direction in which the curve is turning. The vector above then makes sense when viewed in conjunction with the scatterplot for a. In particular, we go from turning down to turning up after the fifth point, and we start turning to the left (with respect to the x axis) after the 12th point.
最后，要获取加速度的切线分量和法线分量，需要相对于t的s，x和y的二阶导数，然后可以获取曲率和其余部分我们的组件(请记住，它们都是t的标量函数):
Finally, to get the tangential and normal components of acceleration, we need the second derivatives of s, x, and y with respect to t, and then we can get the curvature and the rest of our components (keeping in mind that they are all scalar functions of 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
这篇关于numpy曲线曲率的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
