# Measure a DensityHistogram[] pair similarity

I study human vision and more specifically eye-movements.

"If we display 2 symmetrical patterns (20 min one after the other), will our gaze distribution be symmetric is my research question."

The 2 figures at the bottom row below is what is displayed to subjects for 3 seconds. One pattern, then, later in the experiment, its symmetrical transform.

Above are their respective Gaze Density histograms. That is the distribution of where their eyes were while observing the stimuli. The Blue square is the Center Of Gravity of the stimuli.

How can I measure their similarity ? If I have some ideas, I think Mathematica offers great means of image analysis that could be used here.

You could find here the data : allSymFix : 93 sublist for the 93 stimuli pairs I present, along with a manipulate to see all the histograms

allSymFix[[1,1]] are all the gaze observed on stimuli 1 original version.
allSymFix[[1,2]] on its symmetrical transform

How can I measure the similarity within each allSymFix[[original stimuli]].

I will then compare it with the similarity computed on random pairs assembled

-
What is your definition of "similar"? –  rm -rf Feb 17 '12 at 19:38
I can't open the allSymFix.dat.gz file as gunzip refuses to work on it. Besides I would suggest to make an example data file that is not 25MB in its compressed form. What statistics do you need to adhere to? As vision takes place more on a disk than a point I would start of by weighting each pixel and its neighbouring pixel with some disk function. Look at GaussianMatrix and run this over your matrix of observation values and determine for each point its weight. I might even miss understand your question. –  Matariki Feb 17 '12 at 20:02
@500 One idea would be to convert your histograms to strings of 1's and 0's using Binarize[] and then looking for overlap. Also, as Matariki suggests, provide us with a very small portion of your data, so we can work on it faster. @Matariki Mathematica uncompresses the file by itself. –  CHM Feb 17 '12 at 20:02
@Brett, sorry it was a typo. –  500 Feb 17 '12 at 20:04
@CHM Thanks, I didn't know. –  Matariki Feb 17 '12 at 20:18

## 3 Answers

We'll use SmoothKernelDistribution. Correlated pair with left data set reflected around y-axis:

lefTimagE = SmoothKernelDistribution[{-1, 1} # & /@ allSymFix[[3, 1]]];
righTimagE = SmoothKernelDistribution[allSymFix[[3, 2]]];


Visualize in 3D:

  Row@Plot3D[Evaluate[#], {x, -13, 13}, {y, -13, 13}, PlotRange -> All,
MeshFunctions -> {#3 &}, Mesh -> 15, PlotPoints -> 50] & /@
{PDF[lefTimagE, {x, y}], PDF[lefTimagE, {x, y}] PDF[righTimagE, {x, y}],
PDF[righTimagE, {x, y}]}


The middle is overlap - notice small values. Integrate to find total characteristic

NIntegrate[Evaluate[PDF[lefTimagE, {x, y}] PDF[righTimagE, {x, y}]],
{x, -13, 13}, {y, -13, 13}, Method -> "AdaptiveMonteCarlo"]


Answer: 0.00549086

Random pair:

lefTimagE = SmoothKernelDistribution[{-1, 1} # & /@ allSymFix[[3, 1]]];
righTimagE = SmoothKernelDistribution[allSymFix[[15, 2]]];


Visualize in 2D this time for verity:

Row@ContourPlot[Evaluate[#], {x, -13, 13}, {y, -13, 13},PlotRange -> All,
Mesh -> 15, PlotPoints -> 50] & /@ {PDF[lefTimagE, {x, y}], PDF[lefTimagE,
{x, y}] PDF[righTimagE, {x, y}], PDF[righTimagE, {x, y}]}


The middle is overlap. Integrate to find total characteristic

NIntegrate[Evaluate[PDF[lefTimagE, {x, y}] PDF[righTimagE, {x, y}]],
{x, -13, 13}, {y, -13, 13}, Method -> "AdaptiveMonteCarlo"]


Answer: 0.0038788

I liked Andy's analysis of the whole set for his metric. I ran it for my integral metric too:

Correlated pairs:

coRdaT = Table[NIntegrate[Evaluate[PDF[SmoothKernelDistribution[{-1, 1}
# & /@ allSymFix[[k, 1]]], {x, y}] PDF[SmoothKernelDistribution[
allSymFix[[k, 2]]], {x, y}]], {x, -13,13}, {y, -13, 13},
Method -> "AdaptiveMonteCarlo"] , {k, 1, 93}];


Random pairs:

uNcoRdaT = Table[NIntegrate[Evaluate[PDF[SmoothKernelDistribution[{-1, 1}
# & /@ allSymFix[[k, 1]]], {x, y}] PDF[SmoothKernelDistribution[
allSymFix[[RandomInteger[{1, 93}], 2]]], {x, y}]], {x, -13, 13},
{y, -13, 13}, Method -> "AdaptiveMonteCarlo"] , {k, 1, 93}];


Analysis:

SmoothHistogram[{coRdaT, uNcoRdaT}, Filling -> Axis,
PlotStyle -> {{Thick, Blue}, {Thick, Red}}]


Conclusion: on average integral of overlap for correlated pairs almost order of magnitude greater than for random pairs.

======= ARCHIVE: less reliable, needs-polishing approach =======

Here is a very simple take on this. If my understanding is correct, @500 wishes to see spatial correlation between left and right 2D patterns. I'll use SmoothDensityHistogram because it IMO gives better data representation in this case, but you can use your original approach too. The idea is to use ImageMultiply to "amplify" overlapping regions. Midle image is the overlap for a specific set of your data. Note it was ImageAdjust-ed for better visual perception. As numeric measure you have red number (computed before ImageAdjust for uniform scale) The red number is total "intensity" of overlap which could be some sort of correlation measure. BTW we also need to reflect one of the images around vertical axis, otherwise overlap will be meaningless. Here is correlated pair - data set 3, left and right images. As you can see the red number is high and the overlap does look like originals.

ili = SmoothDensityHistogram[#, Background -> Black,
ColorFunction -> GrayLevel, ImageSize -> 300,
PlotRange -> {{-13, 13}, {-13, 13}}, ImagePadding -> 0,
ImageMargins -> 0, PlotRangePadding -> 0, Mesh -> 0] & /@
allSymFix[[3]];
il = {ImageReflect[ili[[1]], Left -> Right], ili[[2]]};
Framed@Labeled[GraphicsRow[Riffle[il, (cori =
ColorConvert[ImageMultiply @@ il, "Grayscale"]) //
ImageAdjust], Spacings -> 1], ImageData[cori] // Total // Total,
Top, LabelStyle -> Directive[Red, Bold, 20]]


And here is random pairing of set 3 left image and set 15 right image. As you can see the red number is much less and the overlap does not look like originals.

ili = SmoothDensityHistogram[#, Background -> Black,
ColorFunction -> GrayLevel, ImageSize -> 300,
PlotRange -> {{-13, 13}, {-13, 13}}, ImagePadding -> 0,
ImageMargins -> 0, PlotRangePadding -> 0,
Mesh -> 0] & /@ {allSymFix[[3, 1]], allSymFix[[15, 2]]};
il = {ImageReflect[ili[[1]], Left -> Right], ili[[2]]};
Framed@Labeled[GraphicsRow[Riffle[il, (cori =
ColorConvert[ImageMultiply @@ il, "Grayscale"]) //
ImageAdjust], Spacings -> 1], ImageData[cori] // Total // Total,
Top, LabelStyle -> Directive[Red, Bold, 20]]


A word of caution: Andy's good comment made me realize there are a few things to worry about here. Most impotently, in most cases our graphics resales the data before it passes them to ColorFunction. This means that for these different data sets their maximums will look same bright on the plots:

Max /@ {allSymFix[[3, 1]], allSymFix[[15, 2]]}
*Answer:* {10.466, 11.172}


This affects correct overlap estimates.

-
this is a fascinating solution +1. My only gripe with any method that uses images is that it will be sensitive to the choice of bandwidth/binning. This one is also quite sensitive to the number of PlotPoints used. I would suggest boosting them. –  Andy Ross Feb 17 '12 at 21:56
@Vitality Kaurov, thank you very much ! –  500 Feb 18 '12 at 0:47
This is really cool, thank you. What a visual representation. –  500 Feb 18 '12 at 1:08
@AndyRoss Thanks! Good thinking btw - I updated the post ;-) –  Vitaliy Kaurov Feb 18 '12 at 8:35
@500 You are welcome, notice I updated the post. –  Vitaliy Kaurov Feb 18 '12 at 8:37

For the sake of having alternative solutions. Here I'm using BinCounts to create a grid of counts that represent the gaze densities.

An important thing is to make sure we are comparing "apples to apples" so I use the same binning for all subjects. Note that the data range from about -20 to 20 in both directions.

bspecX = {-20, 20, 1};
bspecY = {-20, 20, 1};


Now I bin the data using these specs reflecting the symmetric data about the y-axis.

hist1 = Table[BinCounts[allSymFix[[i, 1]], bspecX, bspecY], {i, 93}];
hist2 = Table[BinCounts[Transpose[{-1, 1}*Transpose[allSymFix[[i, 2]]]],
bspecX,bspecY], {i, 93}];


Here is a look at one pair of "histograms" that I am comparing..

Similar to @Vitaliy's method I look at absolute deviations between the two histograms.

scores = MapThread[Total[Abs[(#1 - #2)], Infinity] &, {hist1, hist2}];


Now I also examine these "scores" for randomly chosen pairings...

scoresr =
MapThread[
Total[Abs[(#1 - #2)], Infinity] &, {RandomSample[hist1],
RandomSample[hist2]}];


Lets compare the distributions of the scores for both the pairings and the randomly selected pairings.

SmoothHistogram[{scores, scoresr}, PlotStyle -> {Blue, Red}]


Notice that the statistic does separate the pairings (blue) from the randomly selected pairings (red) but that it isn't a perfect separation. This is probably at least partly due to the "random" pairings being somewhat artificially similar since they are taken from the same (rather small) collection of images.

Though it is statistics on statistics we can perform a formal test for location on the scores which tells us that we have enough evidence to reject the null that these two statistics have equal locations.

LocationTest[{scores, scoresr}, 0, "TestDataTable"]


-
Cool ;) I like that you analysed the whole set +1. If I find some time I'll do the same for my integral measures. –  Vitaliy Kaurov Feb 19 '12 at 9:42
@Andy, thank You very much ! –  500 Feb 19 '12 at 16:27

I know you asked specifically for the comparison of the density histograms, but while thinking about your problem some issues I want to share came to my mind.

Let's first discuss your data. I mentioned this already to Andy in the chat, that one should be aware that the data might not be in the same coordinate system. Do you know for sure, that your reference frame of your sensors is the same for both measurements? If yes, and the origin of your coordinate-system is really in the center, then you can mirror your data to compare the two measurements. Andy and Vitaliy assumed this implicitly but when you look at a plot of your first data-set, it seems a bit as they were shifted.

If you are sure about your setting, this could mean smth. like that our eye is attracted to the same objects (because the overall appearance is similar) but in detail, we prefer to visit your disks from the left, maybe since we are biased because of our reading direction.

I hope it's becoming clear now, what this means for the approach of simply mirroring the data! One completely symmetric stimuli would be a disk at the origin and let's assume, that we have a preferred side which we look first, then you would maybe measure something like this

What would you expect as result for the transformed stimuli? Right, exactly the same. But now you mirror the other stimuli and compare the two outcomes:

While the for such a simple setting the two stimuli should be the same, you will measure considerable differences between them. I was thinking a while about this problem and whether it can be prevented but my knowledge in this area is limited and what follows is only a guess.

If you would present me only a few single pixel, I would have no other chance as just to inspect them without meaning and since they are so small, I would not have the chance to prefer any side. The situation gets complicated, when you present me more pixels in a cluster, because then I see more then just the pixels. It maybe remembers me of something, some object or the silhouette of a girl and then I clearly have my preferences where to look .

Therefore, the question is: Is a disk free of a biased view? And here I would say no, because they are too big, to really look directly in the center. It's really hard to observe oneself but I would say I like the sharp edges more then the flat center. You are in the position to test this. Present your test-persons a simple disk and see what happens. If the observation is not centered, you cannot simply mirror your data and believe you will get reliable results.

Back to your question. My first plot does not show all points of your data. I assumed that when an eye is attracted to something, it sort of rests there for a while. Therefore, the movement is smaller than when the eye makes a movement. This is the reason why I filtered your data by the velocity (first differences) to obtain the spots only:

data = allSymFix[[1]];
{spots1, spots2} = Function[{data}, First[Transpose[
Select[
(Transpose[{Most[#], Norm /@ Differences[#]}]) &[data],
(Last[#] < 0.01) &]
]
]] /@ data;
spots2 = Times[{-1, 1}, #] & /@ spots2;
ListPlot[{spots1, spots2}, AspectRatio -> Automatic,
PlotRange -> {{-15, 15}, {-15, 15}},
PlotStyle -> {{Opacity[0.1], Blue}, {Opacity[0.1], Green}},
Axes -> False, Frame -> True]


Let me stick to the assumption that those spots are more important for the comparison you want. If you like you can extend the approach to all points later. The simplest thing which came to my mind was to use every point in spots1 and find the nearest point in spots2 and calculate their distances. The Total or Mean of this list of distances is then I quite good measure how similar the observations are. Or, if you like you can estimate the parameters of a gamma-distribution on those distances and obtain more information; do what you like:

distances = ParallelMap[Norm[# - First@Nearest[spots1, #]] &, spots2];
Mean[distances]
(*
0.86046
*)


For a symmetric pattern and under the assumption that disks are viewed at the center, the above measure will give a value close to zero for the simple disk above.

While I looked into your data, I was asking myself whether there is more interesting stuff to find out. One interesting point is to investigate in the movement of the eye. I don't know how exact your sensors are, but what I liked to know was: Have eyes a preferred direction in their movement. Since data is time data we can easily calculate the angle of each single eye-movement.

As I first saw the result, I could not believe it. Btw, I still doublt it and I think I do something wrong.

Histogram[
If[Norm[{#1, #2}] > 10.^-3, ArcTan[#1, #2]] & @@@
Differences[data[[1]]], 100, "PDF",
Epilog -> {Red,
Line[{{#, -.3}, {#, .3}}] & /@ Table[i, {i, -Pi, Pi, Pi/4}]}]


It seems our eyes know angles in fractions of Pi better then the most of us do. I marked the Pi/4 positions in red. I would have guessed, that we may have preferences in horizontal and vertical movement, but this? Lucky enough for the transformed stimuli, the histogram looks exactly the same.

Another thing came to my mind: If we have preferences in direction and movement, maybe we like to scan patterns from left to right or top to bottom. Since the order in your data represents the time, we could look at the movements in x- and y-direction separately. Since only x-values are affected by your transformation, it's maybe ok to assume that the y-movement of the eyes shows some similarities. What I do is, I calculated all partial sums of the y values. So if the eye stays at one point for observation, it's a line with the slope of the y-value of this point in the following plot

ListLinePlot[
FoldList[Plus, 0, Last[Transpose[#]]]/Length[#] & /@ data,
Axes -> False, Frame -> True]


One has clearly to think about the value of such a plot, but I like inspecting different sides because humans are good in recognizing patterns. Note, that equal slopes of lines in the above plot show regions, where the same object was observed. You can compare this to the plot of pure y-values and you'll see, that there it's harder to extract features even if you smooth the data. Interesting here, are for instance the region at x=7500 where it seems as there were objects observed in the same order.

-
thank you very much. Your ideas solutions are very interesting. Yes to your question on top, mirroring is correct. What you don`t have are the disks (stimuli) data. If you send me an email laeh 500 at g mail dot com I would be happy to send you the files. Again, thank You –  500 Feb 19 '12 at 16:26
@500 I would really like to, but mathematica.SE is such a time-killer and I have to do a lot stuff for myself. Do you have a complete symmetric stimuli-graphics, which you basically would not need to mirror, since it is the same? –  halirutan Feb 23 '12 at 0:00
just read that, very interesting idea, thank you ! –  500 Mar 3 '12 at 3:44