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

Question

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[[;; #[[1]]]],
 AspectRatio -> Automatic, ImageSize -> Medium, 
 PlotStyle -> Directive[Black, PointSize[0.01], Opacity[#[[2]]]]
] & /@ {{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[[;;#[[1]]]],
  AspectRatio->Automatic,ImageSize->Medium,
  PlotStyle->Directive[Black,PointSize[#[[2]]],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.

Answer 1

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:

Answer 2

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

Now the full answer:

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[#[[1]]]}, 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][[1]] /. 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

Answering @Oleksandr comments, the following does not assume an exponential model:

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

Answer 3

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

Answer 4

In version 10, ListPlot chooses point size automatically. They have used this exact example to demonstrate it here.