Generating a “best fit” curve to a set of noisy data corresponding to a moving object

Question

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:

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.

Answer 1

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}}}

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]

... 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}]

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.