Extracting motion data from a list of coordinates

Question

I have a series of CSV files of timestamped coordinates (X, Y, and Z in mm). What would be the simplest way to extract motion data from them?

Measurables

The information I'd like to extract includes the following:

1. Number of direction changes
2. Initial acceleration of the first and last movements
3. ...and the bearing (angle) of these movements
4. Average speed whilst non-stationary

Ideally, I'd eventually like to be able to categorise patterns of motion, so bonus points for anyone who can suggest a way of doing this. It strikes me that one way I could do this would be to generate pictures/videos of the motion from the coordinates and ask humans to categorise them - suggestions as to how I'd do this are very welcome.

Noise

A complication is the fact that the readings are polluted with noise. In order to overcome this, each recording is prefaced with at least 20 seconds of stillness which can serve as a sort of "noise profile". I'm not sure how to implement this though.

Specifics

If it helps, the motion being recorded is that of a persons hand during a simple grabbing task. The data is generated using a magnetic motion tracker attached to the wrist. Also, I'm using C#, but I guess the maths is language agnostic.

Edits

• Magnetic tracker spec: http://www.ascension-tech.com/realtime/RTminiBIRD500_800.php
• Sample data file: http://tdwright.co.uk/sample.csv

Bounty

For the bounty, I'd really like to see some (pseudo-)code examples.

Thorpe

Do you want to get motion information per axis? – Duo

Ah, sorry. I deleted my previous comments just as you responded. – Duo

No worries, your rephrasing is useful. Thinking about it, I only really care about movement in the XY plane - height above/below the origin is less interesting. – Thorpe

Ok. Would you like to detect motion changes in the x axis and the y axis, separately? Or, would you like to detect derivations from a straight line of motion (i.e. taking into account both the x and y axis) – Duo

Based on my observation of the recordings, it seems like taking both into account would be better: often the direction changes were quite subtle. – Thorpe

You should cross post this to math.stackexchange.com – Educated

Answer 1

Let's see what can be done with your example data.

Disclaimer: I didn't read your hardware specs (tl;dr :))

I'll work this out in Mathematica for convenience. The relevant algorithms (not many) will be provided as links.

The first observation is that all your measurements are equally spaced in time, which is most convenient for simplifying the approach and algorithms. We will represent "time" or "ticks" (measurements) on our convenience, as their are equivalent.

Let's first plot your position by axis, to see what the problem is about:

(* This is Mathematica code, don't mind, I am posting this only for 
   future reference *)
ListPlot[Transpose@(Take[p1[[All, 2 ;; 4]]][[1 ;;]]), 
 PlotRange -> All,
 AxesLabel -> {Style["Ticks", Medium, Bold], 
               Style["Position (X,Y,Z)", Medium, Bold]}]

Now, two observations:

• Your movement starts around tick 1000
• Your movement does not start at {0,0,0}

So, we will slightly transform your data subtracting a zero position and starting at tick 950.

ListLinePlot[
 Drop[Transpose@(x - Array[Mean@(x[[1 ;; 1000]]) &, Length@x]), {}, 950], 
 PlotRange -> All,
 AxesLabel -> {Style["Ticks", Medium, Bold], 
               Style["Position (X,Y,Z)", Medium, Bold]}]

As the curves have enough noise to spoil the calculations, we will convolve it with a Gaussian Kernel to denoise it:

kern = Table[Exp[-n^2/100]/Sqrt[2. Pi], {n, -10, 10}];
t = Take[p1[[All, 1]]];
x = Take[p1[[All, 2 ;; 4]]];

x1 = ListConvolve[kern, #] & /@ 
   Drop[Transpose@(x - Array[Mean@(x[[1 ;; 1000]]) &, Length@x]), {}, 
    950];

So you can see below the original and smoothed trajectories:

Now we are ready to take Derivatives for the Velocity and Acceleration. We will use fourth order approximants for the first and second derivative. We also will smooth them using a Gaussian kernel, as before:

Vel = ListConvolve[kern, #] & /@ 
   Transpose@
    Table[Table[(-x1[[axis, i + 2]] + x1[[axis, i - 2]] - 
         8 x1[[axis, i - 1]] + 
         8 x1[[axis, i + 1]])/(12 (t[[i + 1]] - t[[i]])), {axis, 1, 3}], 
    {i, 3, Length[x1[[1]]] - 2}];

Acc = ListConvolve[kern, #] & /@ 
   Transpose@
    Table[Table[(-x1[[axis, i + 2]] - x1[[axis, i - 2]] + 
         16 x1[[axis, i - 1]] + 16 x1[[axis, i + 1]] - 
         30 x1[[axis, i]])/(12 (t[[i + 1]] - t[[i]])^2), {axis, 1,  3}], 
   {i, 3, Length[x1[[1]]] - 2}];

And the we plot them:

Show[ListLinePlot[Vel,PlotRange->All,
     AxesLabel->{Style["Ticks",Medium,Bold],
                 Style["Velocity (X,Y,Z)",Medium,Bold]}],
    ListPlot[Vel,PlotRange->All]]

Show[ListLinePlot[Acc,PlotRange->All,
     AxesLabel->{Style["Ticks",Medium,Bold],
                 Style["Acceleation (X,Y,Z)",Medium,Bold]}],
     ListPlot[Acc,PlotRange->All]]

Now, we also have the speed and acceleration modulus:

ListLinePlot[Norm /@ (Transpose@Vel), 
 AxesLabel -> {Style["Ticks", Medium, Bold], 
               Style["Speed Module", Medium, Bold]}, 
 Filling -> Axis]
ListLinePlot[Norm /@ (Transpose@Acc), 
 AxesLabel -> {Style["Ticks", Medium, Bold], 
               Style["Acceleration Module", Medium, Bold]},
 Filling -> Axis]       

And the Heading, as the direction of the Velocity:

Show[Graphics3D[
 {Line@(Normalize/@(Transpose@Vel)),
 Opacity[.7],Sphere[{0,0,0},.7]},
 Epilog->Inset[Framed[Style["Heading",20],
        Background->LightYellow],{Right,Bottom},{Right,Bottom}]]]

I think this is enough to get you started. let me know if you need help in calculating a particular parameter.

HTH!

Edit

Just as an example, suppose you want to calculate the mean speed when the hand is not at rest. so, we select all points whose speed is more than a cutoff, for example 5, and calculate the mean:

Mean@Select[Norm /@ (Transpose@Vel), # > 5 &]
-> 148.085

The units for that magnitude depend on your time units, but I don't see them specified anywhere.

Please note that the cutoff speed is not "intuitive". You can search an appropriate value by plotting the mean speed vs the cutoff speed:

ListLinePlot[
 Table[Mean@Select[Norm /@ (Transpose@Vel), # > h &], {h, 1, 30}], 
 AxesLabel -> {Style["Cutoff Speed", Medium, Bold], 
               Style["Mean Speed", Medium, Bold]}]

So you see that 5 is an appropriate value.

Dramatize

Great answer, but you've managed to induce some serious Mathematica envy! Going to try and find a good Gaussian kernel algo for C#... – Thorpe

@Tom The Gaussian filter is a nice one, but perhaps other filters will do. Don't get stuck at it. – Dramatize

I think belisarius really deserves the bounty! – Salot

@gpu_drug I agree completely. Thanks @belisarius! – Thorpe

Answer 2

e solution could be as simple as a state machine, where each state represents a direction. Sequences of movements are represented by sequences of directions. This approach would only work if the orientation of the sensor doesn't change with respect to the movements, otherwise you'll need a method of translating the movements into the correct orientation, before calculating sequences of directions.

On the other end, you could use various AI techniques, although exactly what you'd use is beyond me.

To get the speed between any two coordinates:

               _________________________________
Avg Speed =   /(x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2
            --------------------------------------
                 (t2-t1)

To get the average speed for the whole motion, say you have 100 timestamped coordinates, use the above equation to calculate 99 speed values. Then sum all the speeds, and divide by the number of speeds (99)

To get the acceleration, the location at three moments is required, or the velocity at two moments.

Accel X = (x3 - 2*x + x1) / (t3 - t2)
Accel Y = (y3 - 2*y + y1) / (t3 - t2)
Accel Z = (z3 - 2*z + z1) / (t3 - t2)

Homiletics

Speed is always scalar. Velocity is a vector. – Duo

The spec is a long, dull document - what info do you want? I'll upload a sample data file if it'll help though. – Thorpe

Your average speed (which can be further generalised to include z coordinates) doesn't work if the start and end points are the same. It could be applied at each step (pair of coordinates), but this would still not account for moments of non-movement. – Thorpe

Answer 3

Note: This all assumes per axis calculations: I have no experience with two-axis particle motion.

You will have a much easier time with this if you first convert your position measurements to velocity measurements.

First step: Remove the noise. As you said, each recording is prefaced with 20 seconds of stillness. So, to find the actual measurements, search for 20 second intervals where the position doesn't change. Then, take the measurement directly after.

Second step: Calculate velocities using: (x2 - x1)/(t2 - t1); the slope formula. The interval should match the interval of the recordings.

Calculations:

Direction change:

A direction change occurs where the acceleration is zero. Use numeric integration to find these times. Integrate from 0 until a time when the result of the integration is zero. Record this time. Then, integrate from the previous time until you get zero again. Repeat until you hit the end of the data.

Initial accelerations:

These are found using the slope formula again, substituting v for x.

Average speed:

The average speed formula is the slope formula. x1 and t1 should correspond to the first reading, and x2 and t2 should correspond to the final reading.

Duo

Like I said, this isn't applicable to 2d motion. Hopefully someone who's better at math than I am can apply this to 2d motion :) – Duo