# How can the opacity of a dense data set be automatically chosen?

### The problem

In a dense data set plotted with opaque points, data points can pile on top of each other. This can cause serious problems with the interpretation of the plots. I'll use as an example a scatterplot of normally distributed points.

theData = RandomReal[NormalDistribution[], {100000, 2}];


Unless one is careful, these plots end up showing a solid disk around the origin, which obscures the distribution since point densities above a certain level appear clipped. This also gives inappropriate emphasis to outliers, by plotting them with the same intensity as what might be a pile of a million points close together at the origin.

Here is the same plot shown at different sizes.

Row[ListPlot[
theData,
AspectRatio -> Automatic, ImageSize -> #, PlotStyle -> Directive[Black]
] & /@ {90, 180, 360}] An artificial black disk appears at the origin and its apparent size depends on the plot size, even though the data are identical.

Here is the same plot with different settings for PointSize.

Row[ListPlot[
theData,
AspectRatio -> Automatic, ImageSize -> Medium, PlotStyle -> Directive[Black, PointSize[#]]
] & /@ {0.01, 0.005, 0.0025}] The disk appears a different size depending on point size, even though the data are identical.

Here are plots of successively smaller subsets of the data.

Row[ListPlot[
theData[[;; #]],
AspectRatio -> Automatic, ImageSize -> Medium, PlotStyle -> Directive[Black]
] & /@ {100000, 50000, 25000}] The distributions are obviously all the same, but they look different in the plot because of the pileup of points.

### Solving pileup by changing opacity

For a given point size, one solution is to choose the minimum opacity where the darkest pixel in the plot is black. Here is an ugly, slow, procedural, brute-force way of choosing the correct opacity:

Block[
{darkestPixel},
darkestPixel=1;
idealOpacity=0;
While[
darkestPixel>0,
idealOpacity=idealOpacity+0.001;
darkestPixel=Min@Flatten@ImageData@ColorConvert[
Rasterize@ListPlot[
theData,
AspectRatio->Automatic,ImageSize->Medium,
PlotStyle->Directive[Black,PointSize[0.01],Opacity[idealOpacity]],Axes->False
],
"Grayscale"
]
]
]


Here are the resulting plots for different subsets of the data:

Row[ListPlot[
theData[[;; #[]]],
AspectRatio -> Automatic, ImageSize -> Medium,
PlotStyle -> Directive[Black, PointSize[0.01], Opacity[#[]]]
] & /@ {{100000, 0.018}, {50000, 0.042}, {25000, 0.1}}] ### Solving pileup by changing point size

Another way is to choose a fixed opacity and increase the size of the points until the darkest pixel is black.

Block[
{darkestPixel},
darkestPixel = 1;
idealPointSize = 0;
While[
darkestPixel > 0,
idealPointSize = idealPointSize + 0.001;
darkestPixel = Min@Flatten@ImageData@ColorConvert[
Rasterize@
ListPlot[theData,
AspectRatio->Automatic, ImageSize->Medium,
PlotStyle->Directive[Black,PointSize[idealPointSize],Opacity[0.05]],Axes->False
],
"Grayscale"
]
]
]


The results for different subsets of the data:

Row[
ListPlot[
theData[[;;#[]]],
AspectRatio->Automatic,ImageSize->Medium,
PlotStyle->Directive[Black,PointSize[#[]],Opacity[0.05]]
]&/@{{100000,0.003},{50000,0.01},{25000,0.015}}
] ### The question

What's the best way to quickly and automatically choose the correct opacity, point size or combination of the two so that there is no "clipping" in a dense plot? Is it possible to do this inside of PlotStyle instead of as a pre-processing step? For bonus points, explain the best way to do this for multiple data sets on the same plot.

• Adjusting Opacity/PointSize to fit a plot with all data is nice, but there might be limitations. I myself encountered one when dealing with a similar data visualization problem, where the data is so dense, I used a very small Opacity and end with a discontinuity of the color.. – Silvia Aug 8 '12 at 18:06

At the risk of stating the obvious (and also not directly answering the question):

A ListPlot with opacity less than 1 is very very similar to a histogram, so why not just use that?

theData = RandomReal[NormalDistribution[], {10000, 2}];
opts = {ColorFunction -> Function[c, GrayLevel[1 - c]],
PlotRange -> {{-4, 4}, {-4, 4}}, ImageSize -> Medium};
Row@{
DensityHistogram[theData, opts],
SmoothDensityHistogram[theData, opts]
} (When doing histogram of large data sets in mma 8 please be aware, that there is a memory leak in HistogramList. See this answer)

EDIT So to also answer the original question. From the histogram it should be possible to estimate the opacity and point size (although I would not recommend this as the rendering is much slower then a normal histogram). Here is my take on this

{bins, counts} = HistogramList[theData];
maxcounts = Max@counts;
binarea = Times @@ Part[Differences /@ bins, All, 1];
r = N@Sqrt[binarea/Pi];
s = 10;
ListPlot[theData,
PlotStyle ->
Directive[Black, Opacity[1/maxcounts*s], PointSize[r/s]],
ImageSize -> #] & /@ {90, 180, 360}


10 000 points: 100 000 points: • theData = RandomReal[NormalDistribution[], {10000, 2}]; i = Image@ ListPlot[theData, PlotStyle -> {Black, Opacity[.3]}, Axes -> False] Lighter@Lighter@i – Dr. belisarius Aug 7 '12 at 21:18
• I was in the middle of writing my own routine to count number of points per "sqaure" (ie to coarse-grain, or bin, the data) when I realized mma must have this built-in. then I saw this! +1 – acl Aug 7 '12 at 22:41
• @verde Did you forget a semicolon after the list plot? Without it it just exceeds the recursion limit... – Ajasja Aug 8 '12 at 7:48
• Nope, but there is a <return> after the first plot that is not showed in the comment – Dr. belisarius Aug 8 '12 at 14:14

tl;dr

Final results first:

(*Function Definition*)
ClearAll[opaFun];
Options[opaFun] = Options[ListPlot];
opaFun[points_, opts : OptionsPattern[]] := Module[{f, steps = 10 },
f[x_] := Min[Norm /@ Flatten[ImageData@
ListPlot[points, opts, Axes -> False, PlotStyle -> {Black, Opacity[x]}],1]]/Sqrt@3;
Return@NestWhileList[{#, f[#]} &@(#[] + 1/steps) &, {0, 1},
UnsameQ @@ ({##}[[All, 2]]) &, 2, steps][[-2, 1]]
];
(*Usage*)
opts = {AspectRatio -> 1/2, ImageSize -> 300, Axes -> False};
Timing@ListPlot[#, opts, PlotStyle -> {Black, Opacity[opaFun[#]]}] &@
RandomReal[NormalDistribution[], {10000, 2}]


## A little but time consuming experiment first. You'll not need to do it, as it is only for finding a model.

Let's see how the Min value of the plot intensity varies with Opacity:

theData = RandomReal[NormalDistribution[], {2000, 2}];
f[x_] := f[x] =
Min[Norm /@ Flatten[ImageData@Rasterize[
ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200,
PlotStyle -> {Black, Opacity[x]}, Axes -> False]], 1]];
Plot[f[x] , {x, 0, .4}, PlotRange -> Full] So, it is an exponential. (Note: in the last edit to this post I got rid of the exponential model by fitting a few points with an Interpolation, and it works pretty nice)

Let's fit it:

data = Table[{i, f[i]}, {i, 0, 1, .1}]

model = a Exp[b x];
fit = FindFit[data, model, {a, b}, x];
modelf = Function[{t}, Evaluate[model /. fit]]

Show[ListPlot@data, Plot[modelf[x], {x, 0, 1}]] Now you are ready to set the min value of the brightness of the plot to whatever you want:

(The Sqrt@3 is a normalization factor for the intensity of the {1,1,1} RGB pixel.)

Let's use it:

opac = x /. Solve[# == a E^(b x)/Sqrt@3, x] /. fit & /@ {1/2, 1/4, 1/20, 1/200}

ListPlot[theData, AspectRatio -> Automatic, ImageSize -> 200,
PlotStyle -> {Black, Opacity[#[]]}, Axes -> False] & /@ opac Edit

Let's pack that all together in a function and plot two very different point sets with the same maximum darkness.

ClearAll[opa];
Options[opa] = Options[ListPlot];
opa[desiredOpacity_, points_, opts : OptionsPattern[]] :=
Module[{f, a, b, model, fit, modelf, x},
f[x_] :=
f[x] = Min[
Norm /@ Flatten[
ImageData@
Rasterize[
ListPlot[points, Axes -> False,
PlotStyle -> {Black, Opacity[x]}]], 1]];

model = a Exp[b x];
fit = FindFit[Table[{i, f[i]}, {i, 0, 1, .1}], model, {a, b}, x];
modelf = Function[{t}, Evaluate[model /. fit]];
Return[x /. Quiet@Solve[# == modelf[x]/Sqrt@3, x][] /. fit &@
desiredOpacity];
]

theData  = RandomReal[NormalDistribution[], {2000, 2}];
theData1 = RandomReal[NormalDistribution[], {10000, 2}];

opad = opa[.5,  theData,  AspectRatio -> Automatic, ImageSize -> 200];
opad1 = opa[.5, theData1, AspectRatio -> Automatic, ImageSize -> 200];

Grid[{{ListPlot[theData, Axes -> False, PlotStyle -> {Black, Opacity[opad]},
AspectRatio -> Automatic, ImageSize -> 200],
ListPlot[theData1,Axes -> False, PlotStyle -> {Black, Opacity[opad1]},
AspectRatio -> Automatic, ImageSize -> 200]}},
Frame -> All] The same, but darker Remember that the default plot is: Edit

ClearAll[opa];
Options[opa] = Options[ListPlot];
opa[desiredOpacity_, points_, opts : OptionsPattern[]] :=
Module[{f, model, x},
f[x_] := f[x] =
Min[Norm /@ Flatten[ImageData@ Rasterize[
ListPlot[points, Axes -> False, PlotStyle -> {Black, Opacity[x]}]], 1]];
model = Interpolation[Table[{i, f[i]}, {i, 0, 1, .1}]];
x /. FindRoot[desiredOpacity == model[x]/Sqrt@3, {x, 0, 1}][]]

• Don't the model parameters depend on the distribution of the data? It seems to me that, as you'd have to fit the model for every distribution that you want to plot, the time-intensive step can't be factored out in general. – Oleksandr R. Aug 8 '12 at 9:05
• @OleksandrR. The time intensive step was an experiment. For each distribution you need an appropriate model, yes but you may "guess" it from an interpolation dome over data = Table[{i, f[i]}, {i, 0, 1, .1}] instead of fitting the data like I did here. Posting that generalization later – Dr. belisarius Aug 8 '12 at 12:20
• I think I've now convinced myself that it'll always be exponential, but the parameters will still need to be determined each time. You could also use your current Plot approach with PlotPoints/MaxRecursion rather than a fixed step. – Oleksandr R. Aug 8 '12 at 12:34
• @OleksandrR. Getting the parms for a specific distribution is fast. Unless you have to do that a thousand times, I guess it is not a problem – Dr. belisarius Aug 8 '12 at 14:10
• @OleksandrR. see last edit using Interpolating Function – Dr. belisarius Aug 8 '12 at 14:26

Maybe the clipping can be avoided by not using ListPlot. Here is my variation using Point that you might find helpful (or maybe not), no time for extensive testing of all your conditions:

dist = Rescale[1 - EuclideanDistance[{0, 0}, #] & /@ theData];

Graphics[{Thick,
Point[theData,
VertexColors -> (Blend[{{0, {Black, Opacity[.01]}}, {1, {Black,
Opacity[.1]}}}, #1] & /@ dist)]}, AspectRatio -> 1,
Axes -> True, PlotRange -> {{-4, 4}, {-4, 4}}] In version 10, ListPlot chooses point size automatically. They have used this exact example to demonstrate it here.

• That's pretty cool! I wonder what happens if you have two different datasets in the same plot -- are their point sizes adapted independently or are they the same? I can see arguments in favour of either. – Rahul Jul 9 '14 at 20:15