Clustering of space-time data

Question

Below is an example of a gaze sequence I recorded during a 3 seconds display. That is, where the eye was at every millisecond. While we should have 3000 points, we are missing some due to blinking.

Fixation or visual fixation is the maintaining of the visual gaze on a single location. I need to extract those fixations. That is group of gazes that are both contiguous in time and space. Below is the location of fixations. Of course we'll have to implement thresholds.

In this file available for download, you will find 7 sequences there if you simply open and run the attached notebook.

gazeSeq[1] is the first out of 7 sequences made out of 3000 sublists representing each gaze record as {x,y,time}:

gazeSeq[1][[1]] 
Out[1]= {-0.562, 0.125, 1000.}

where time goes from 1000 to 4000 corresponding to the 3 seconds of display. As said above some data points might be missing.

While I tried to use GatherBy, I could not manage to include the condition of "time contiguity" and would get grouping of gazes that did not happen during the same interval during the trial.

Answer 1

Finding the cluster centers is the hard part. There are zillions of ways to do this, such as standardizing $(x,y,t)$ and applying some (almost any) kind of cluster analysis. But these data are special: the eye movement has a measurable speed. The gaze is resting if and only if the speed is low. The threshold for "low" is physically determined (but can also be found in a histogram of the speeds: there will be a break just above 0). That yields the very simple solution: fixations occur at the points of low speed.

It's a good idea to smooth the original data slightly before estimating the speeds:

data = Transpose[Import["f:/temp/gazeSeq_1.dat"]];
smooth = MovingAverage[#, 5] & /@ data;
delta = Differences /@ smooth;
speeds = Prepend[Norm[Most[#]]/Last[#] & /@ Transpose[delta], 0];
Histogram[Log[# + 0.002] & /@ speeds]

The bimodality is clear. The gap is around $\exp(-6)-0.002\approx 0.0005$. Whence

w = Append[smooth, speeds];
ListPlot[#[[1 ;; 2]] & /@ Select[Transpose[w], Last[#] < 0.0005 &], 
 PlotStyle -> PointSize[0.015]]

There are the gaze fixations. having found them, the clustering is (almost) trivial to do (because each fixation now exists as a contiguous sequence of observations in the original data, from which it is readily split off: this respects the time component as well as the spatial ones). This method works beautifully for all seven of the sample datasets.

One advantage of this approach is that it can detect clusters of very short gaze fixations, even those of just two points in the dataset. These would likely go unnoticed by most general-purpose or ad hoc solutions. Of course you can screen them out later if they are of little interest.

Answer 2

I think FindClusters is ideal tool. It just needs slight tweaking. One of your data sets:

data = gazeSeq[3][[All, ;; 2]];

This works:

Show[
   ListPlot[#, PlotStyle -> PointSize[.003]],
   Graphics[{Red, Thick, Circle[#, .3]}& /@ Mean/@ #], 
   AspectRatio -> Automatic, PlotRange -> All, 
   Frame -> True, Axes -> False
   ]& @ Select[FindClusters[data, 15], Length[#] > 100 &]

Another method (accounting for time too) is pure data manipulation. Let's try another data set:

data = gazeSeq[7][[All, ;; 2]];

Look at EuclideanDistance and see that peaks split your data set in the clusters you want:

ListLinePlot[eudata = EuclideanDistance @@@ Partition[data, 2, 1], 
PlotRange -> All, AspectRatio -> 1/5]

It is obvious from this plot that for all your data clusters have more than ~100 points with EuclideanDistance being more than ~0.1. Use this to cluster your data:

clustered = #[[All, 2]] & /@ Select[SplitBy[Transpose[{eudata, 
Most[data]}], #[[1]] < .1 &], Length[#] > 100 &];

And get the result you need:

 Show[ListPlot[#, PlotStyle -> PointSize[.003]], 
 Graphics[{Red, Thick, Circle[#, .3]} & /@ Mean /@ #], 
 AspectRatio -> Automatic, PlotRange -> All, Frame -> True, 
 Axes -> False] &@clustered

Answer 3

Here is an image processing solution that gives you the following result for the 7 different datasets:

Approach:

First, we plot it without any frames or axes and convert to a binary image.

plotRange = Function[xy, #[gazeSeq[2][[All, xy]]] & /@ {Min, Max}] /@ {1, 2};
img = Image@ListPlot[gazeSeq[2][[All, ;; 2]], Axes -> False, PlotRange -> plotRange];
c = img // Binarize // ColorNegate

Then we remove those components that have less than a certain number of pixels after eroding it:

s = DeleteSmallComponents[c~Erosion~1, 20];

Now that we've narrowed down to the clusters, get the centroids of the clusters. You can then visualize the locations of the clusters with Graphics primitives.

m = ComponentMeasurements[s, {"Centroid"}];
Graphics[{Red, Thick, Circle[#, 5] & /@ (m[[All, 2, 1]])}]

Finally, we need to convert the centroids in the pixel coordinates to the plot coordinates and plot it with ListPlot, which will give you the figure shown in the beginning (in this case, dataset 2)

p = With[{dim = ImageDimensions@img}, 
        {Rescale[#[[1]], {1, dim[[1]]}, plotRange[[1]]], 
         Rescale[#[[2]], {1, dim[[2]]}, plotRange[[2]]]} & /@ m[[All, 2, 1]]
    ]

ListPlot[gazeSeq[2][[All, ;; 2]], AspectRatio -> 1,
   Epilog -> First@Graphics[{Red, Thick, Circle[#, 0.3] & /@ p}]]

Answer 4

Try this: First, set the distance threshold.

d = 0.1;

The main function uses Fold, which, along with its companion FoldList and MapThread, is one of the most useful "functional" functions in the language.

test = 
  Fold[If[EuclideanDistance[Most@#2, Mean[Most /@ Last[#1]]] < d, 
     Join[Most[#1], {Join[Last[#1], {#2}]}], Join[#1, {{#2}}]] &,
   {Take[gazeSeq[1], 2]}, Drop[gazeSeq[1], 2]];

Essentially what this is doing is taking each point sequentially, comparing it to the mean point of the "last" cluster already built. If the Euclidean distance is within some prespecified threshold, then the new point joins that cluster. Otherwise it starts its own cluster. Notice I am stripping out the time value with Most before calculating the distance - this is important!

You end up with a bunch of clusters and a bunch of relatively disconnected points. (Output is truncated.)

Length /@ test

{84, 4, 9, 4, 3, 260, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 9, 123, 18, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 189, 45, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 5, 2, 2, 67, 141, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 6, 14, 263,... 11, 5, 48, 91, 48}

So you only want the ones that are actually clusters.

clusters = Select[test, Length[#] > 10 &];

Which look like this.

ListPlot[Map[Most, clusters, {-2}], Frame -> True]

It would be pretty straightforward to build this up into a nice function with the distance threshold as a parameter.

There is probably a neat image processing way to do this too, but I'll leave that to Heike.

---

edit in response to whuber's suggestion

Here is an alternative that captures the idea that speed between points matters as well as distance from the center of the cluster. A new cluster starts if either the point is too far away from the center or if it is too far from the previous point. Notice that I've used a smaller threshold for the pairwise sequential distance test as the distance from the centre.

test2 = Fold[
   If[EuclideanDistance[Most@#2, Mean[Most /@ Last[#1]]] < d || 
      EuclideanDistance[Most@#2, Most[#1[[-1, -1]]] ] < 0.5 d, 
     Join[Most[#1], {Join[Last[#1], {#2}]}], Join[#1, {{#2}}]] &,
   {Take[gazeSeq[1], 2]}, Drop[gazeSeq[1], 2]];

This version does a bit better on avoiding overlapping fixations that are really a single fixation.

 Length /@ test2

{84, 4, 15, 261, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 156, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 234, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 5, 2, 2, 208, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 283,... 192}

clusters2 = Select[test2, Length[#] > 10 &];

ListPlot[Map[Most, clusters2, {-2}], Frame -> True]