# Clustering of space-time data

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.

-
Data link is not working! –  s.s.o Sep 18 '12 at 9:08
@s.s.o I found copies of the files on a backup drive and have shared them on Dropbox. –  whuber Jul 29 '14 at 15:46

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.

-
thank You very much ! It seems you know Eye-Movements pretty well ! –  500 Mar 8 '12 at 0:05
By sheer coincidence I happened to have used this velocity filter to find the gaze groups as well and ComponentMeasurements to cluster them. FindClusters didn't do the trick. I found the filtered points belonging to 1 gaze group to be not always contiguous in time, so I wonder what trivial clustering method you had in mind. –  Sjoerd C. de Vries Mar 8 '12 at 0:06
Almost any clustering method that uses speed and time along with x and y ought to work (e.g., k-means after standardizing all four variables). The hierarchical ones do a nice job (even applied directly to the original data) but seem to require a lot of fussy post-processing. I admit to being somewhat frustrated with Mathematica's clustering interface; almost any statistical package (even the command-driven R) is more congenial. For instance, R will create the dendrogram in one command and let you cut it at a threshold with another command and you're done. –  whuber Mar 8 '12 at 0:14
Big +1. Once again, I am impressed by the combination of depth and simplicity of your approach. –  Leonid Shifrin Mar 8 '12 at 0:33
@whuber, may I ask you a favor : How could I get the duration for each fixation ? –  500 Mar 8 '12 at 3:19

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


-
An out-of-the-box solution like this is attractive. But how do you obtain the parameters 100 (lengths of clusters) and 15 (number of clusters)? It would seem that using these as omnibus parameter values will (a) overlook any short-term fixations and (b) not produce correct answers in datasets with more than 15 clusters. Your solution also does not appear to take time into account. –  whuber Mar 8 '12 at 0:23
@whuber This is why I used term "tweaking", - but good points. I give some tips on the parameter values in the updated 2nd half of the post. I love simple methods with minimum code and sometimes studding and knowing your data properties can significantly reduce amount of code. Same goes for the 2nd method I used taking time in account. It's probably similar to what others are doing here. –  Vitaliy Kaurov Mar 8 '12 at 7:55

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]


-
Why doesn't this solution create an artifact at the beginning of the sequence (assuming the gaze is not yet fixed when the sequence begins)? –  whuber Mar 7 '12 at 23:50
Because I delete short "clusters". Also I think in the first sequence, the gaze was fixed. I actually should have used a slightly larger $d$ threshold. –  Verbeia Mar 8 '12 at 0:17
Thanks for pointing that out. But now I'm worried: what if short fixations actually exist in the dataset? I think there may be an easy fix, though: rather than deleting short clusters, delete those that are large (e.g., have a large total variance in x and y) and short. I think such a criterion would have universal application; that is, a scientist could decide once and for all how large and how short are generally appropriate for any eye movement study. I do still worry that your algorithm could split valid but small clusters into two and thereby fail to identify them at all. –  whuber Mar 8 '12 at 0:26
Well, yes, my algorithm (if you can call the 15 minutes I spend faffing around getting the Join function right) is sensitive to the $d$ threshold. –  Verbeia Mar 8 '12 at 1:21
Mean/@clusters is the perfect information! –  500 Mar 8 '12 at 4:24

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


-
I like this idea (it's related to "RegionGroup" in the GIS world and is quite effective), but it has an interesting and important failure mode: when the eye gazes twice at the same point at different times, this method finds only one cluster. This is exactly what the OP is concerned about. You can fix this up by re-interpreting your code as finding cluster centers in space; follow-up processing using almost any clustering algorithm based on the (x,y,time) data will then correctly find all the clusters situated in both space and time. –  whuber Mar 7 '12 at 23:46
@whuber Thanks, I was about to remark something along those lines. For example, ListPointPlot3D[FindClusters[gazeSeq[2]]] shows how well a straightforward call to FindClusters separates the clusters in space and time and from then on, it's only a matter of manipulation for the OP to find the fixation point. –  The Toad Mar 7 '12 at 23:53