Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
261 views
in Technique[技术] by (71.8m points)

python - Curve curvature in 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 ]])

And the curve looks like this:

plt.scatter(a[:,0], a[:,1])

enter image description here

Question:

How can I calculate the tangential and the radial aceleration vectors at each point? I found some formulas that might be relevant:

enter image description here

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?

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

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.

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 ]])

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)])

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 ]])

which makes sense when glancing at the scatterplot of a.

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)

This yields the following array:

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])

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.

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

This yields the following numpy array:

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   ]])

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:

In [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.])

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

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]])

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.

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

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...