Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Please note that this is an extension of a previous question: Interpolating a sparse list of two-dimensional coordinates

We have an ordered array of coordinates corresponding to the position of an object in a set of camera frames $(f_1,...,f_N) \in F$, which we're tracking with an imperfect object recognition algorithm. As a result, in any given frame, we sometimes fail to identify the object, and we sometimes identify multiple copies of the object. We're left with noisy data that looks something like this:

OrderedArray = {{{70.8938, 216.539},{70.89,216.54}}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {{71.0656,216.573}}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {{67.6546, 220.338}}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {{70.9211, 216.364}}, {{70.9184, 216.346}}, {{70.9079, 216.349}}, {{70.9046, 216.335}}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {},{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {{70.951, 216.705}}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {{70.9621,216.586}}, {{70.918, 216.576}}, {{70.9116, 216.559}}, {{70.9189,216.581}}, {{70.9115, 216.565}}, {{70.9294, 216.552}}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {{67.0276, 218.154}}};

Where the position of a coordinate in the array corresponds to a particular image frame. Here's a listplot output for a larger version of my data set:

enter image description here

I'd like to generate a "best guess" for the motion of my object from frame to frame provided this noisy data, and with the constraint that the object cannot move more than a Euclidean distance $D$ from frames $f_i$ to $f_{(i+1)}$.

How might I best use the intepolation tools in Mathematica 9 to "guess" the shape of the curve followed by the object? I'd ideally like to specify that the object follows a linear trajectory, or, say, a trajectory characterized by a polynomial of bounded degree.


Update: In response to Thies Heidecke's comment on filters, here, I'm tracking a moving object that changes direction fairly slowly. The curve fitting the data should be fairly smooth.

share|improve this question
    
To clarify, is the x axis a time/frame unit and the y axis is a position or are both x and y positions in the two dimensions of the image? –  bobthechemist Jul 27 '13 at 9:51
    
@bobthechemist The coordinates correspond to physical positions of the object on a two-dimensional (i.e. $R^2$) plane. The "time" axis corresponds to the position of the coordinate in the array. "{}" entries in the array, for example, correspond to time points where the object wasn't observed. –  Sparse Pine Jul 27 '13 at 9:55
1  
Two often used tools for that are Kalman Filters and Particle Filters. Since your data is fuzzy, interpolation is not really the best option, since interpolation assumes that your data samples are all correct and tries to fill in the values inbetween (in time). –  Thies Heidecke Jul 27 '13 at 10:47
    
@ThiesHeidecke Hmm... that sounds promising. –  Sparse Pine Jul 27 '13 at 10:53
    
Extending @ThiesHeidecke comment: The choice of the best algorithm depends on the model you use to describe the movement of whatever it is you're tracking. The Kalman filter is the optimal estimator (in a mean-square-error sense) for Gaussian (and therefore linear) systems, although they can be extended to (approximately) deal with nonlinear systems (eg. extended KF). I don't know much about particle filters but from what I understand they are not subject to these restrictions. –  sebhofer Jul 27 '13 at 11:05

1 Answer 1

up vote 0 down vote accepted

Well, I spent a bit of time playing around with this problem and I have an answer, but not necessarily a good one.

First, I take your OrderedArray, append a time unit and then delete the empty entries.

culled = DeleteCases[
   Table[{t, OrderedArray[[t]]}, {t, 1, Length[OrderedArray]}], 
   x_ /; x[[2]] == {}];
culled[[All, 2]] = Mean /@ culled[[All, 2]];

Mr. Wizard and Kuba have more elegant ways of handling your data as described in your other question. I just went with what worked for me. An important point to make with respect to your data is that there seems to be a fair amount of bouncing around. Here are your data plotted with arrows, with the right-hand graphic showing a closeup of the dense region:

Graphics[Arrow /@ 
    Table[{culled[[i, 2]], culled[[i + 1, 2]]}, {i, 1, 
      Length[culled] - 1}], 
   PlotRange -> #] & /@ {All, {{70.8, 71.2}, {216.3, 216.8}}}

Mathematica graphics

I think this randomness(?) is yet another reason why there is a problem fitting your data (in addition to the noise due to data acquisition). In any case, I am going to assume that the x and y positions are both functions of time and can be represented by the linear expressions xpos = mx t + bx and ypos = my t + by. To solve your problem para-metrically, I am using a trick I learned from Oleksander's answer to another question.. In brief, we make a model that has the four parameters (mx, my, bx, by) and two variables (i and t) where i indicates whether we are fitting the xpos or ypos.

transformedData = 
 Map[{{#[[1, 1]], #[[2, 1]], #[[2, 2, 1]]}, {#[[1, 2]], #[[2, 
        1]], #[[2, 2, 2]]}} &, {ConstantArray[{1, 2}, Length[culled]],
      culled} // Transpose]~Flatten~{1, 2}
model[mx_, my_, bx_, by_][i_, t_] := 
  Thread[{mx, my} t + {bx, by}, List][[i]];

We can now run a NonlinearModelFit that will throw a bunch of errors, but will give us a solution nonetheless:

nlm = NonlinearModelFit[transformedData, 
  model[mx, my, bx, by][i, t], {mx, my, bx, by}, {i, t}]
nlm["BestFitParameters"]
(* {mx -> -0.00448236, my -> 0.00042076, bx -> 71.0688, by -> 216.817} *)

After all this work, we can compare the fit (in blue) to your data (in red):

fittedpoints = ( ({mx, my} t + {bx, by}) /. 
      nlm["BestFitParameters"]) /. t -> culled[[All, 1]] // 
   Transpose // Quiet
ListPlot[{culled[[All, 2]], fittedpoints}, 
 BaseStyle -> PointSize[0.02], PlotStyle -> {Red, Blue}, 
 PlotRange -> All]

Mathematica graphics

... and see that the fit is pretty awful. The two possibilities are my solution stinks, or your representative data st... is incomplete :-). To check, I used an answer to this question about recovering points from an image and applied it to your complete data set. Now this isn't an exact representation since I have no idea what data are dropped and therefore don't have a good estimate of the time sequence, so I just assumed that the order of the points provided by the image digitization was the "correct" order.

Just for completeness, here's how I grabbed the data, which is stored in curvLoc:

image = Import["http://i.stack.imgur.com/Ub1DT.png"];
idata = Round[ImageData[image], 1];
col = DeleteDuplicates[Flatten[Round[ImageData[image], 1], 1]];
Graphics[{RGBColor[#], Disk[]}, ImageSize -> Tiny] & /@ col
binImage = Image@Replace[idata, {col[[3]] -> 1, _ :> 0}, {2}]
curve = ImageApply[{0, 0, 0} &, binImage, 
  Masking -> ColorNegate[Binarize[GaussianFilter[binImage, 5]]]]
curvLoc = (Reverse /@ 
    Position[ImageData[curve, DataReversed -> True], {1., 1., 1.}]);

I then applied the the same procedure as above:

t1 = Table[{i, curvLoc[[i]]}, {i, 1, Length[curvLoc]}];
tcurve = Map[{{#[[1, 1]], #[[2, 1]], #[[2, 2, 1]]}, {#[[1, 2]], #[[2, 
         1]], #[[2, 2, 2]]}} &, {ConstantArray[{1, 2}, Length[t1]], 
      t1} // Transpose]~Flatten~{1, 2};
nlm2 = NonlinearModelFit[tcurve, 
  model[mx, my, bx, by][i, t], {mx, my, bx, by}, {i, t}]
nlm2["BestFitParameters"]

(* {mx -> -0.356423, my -> 0.435831, bx -> 564.9, by -> -14.0573} *)

pfit  = ListPlot[Table[model[mx, my, bx, by][All, t] /. 
nlm2["BestFitParameters"], {t, 0, 1000}], Epilog -> 
{Red, PointSize[0.02], Point@curvLoc}]

Mathematica graphics

And we get a decent fit. Note, I didn't do the transformation as described at the beginning of this answer but it's easy enough to reproduce if desired.

Concluding thoughts

So this is one, somewhat painful, method of modeling the position of you object. It is cumbersome but provides you with some idea of how the particle is moving in the x and y directions. You can extract a distance traveled from the slopes, mx and my multiplied by the time between two points and putting those numbers into EuclideanDistance. Depending on what the application is, you might be better off performing a linear best fit of the x,y data presented. As mentioned in the comments, the scatter and the possible random-walk nature of your data might require a fairly advanced treatment. Again, it depends largely on what you want to do with the data, and whether the scatter in your data points is due to data acquisition or particle movement.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.