Plotting data points: Optimizing size and visuals

Question

The following problem is one that may have showed up to you while plotting a large data collection. Suppose you have a data set of say 1 million points. ListLinePlot will use all of those points to create your plot. Unfortunately, if the points are so close together they will overlap in the plot and they become irrelevant.

I was wondering if there is some quick way in Mathematica to get rid of the redundancies in the plots. The term redundancies depends of course of the image resolution (size).

Let us work with the following example

rx[n_] := Accumulate[Prepend[RandomVariate[ExponentialDistribution[1000], n], 0]]
ry[n_] := Accumulate[Prepend[RandomVariate[NormalDistribution[0, .001], n], 0]]
x = rx[100000];
y = ry[100000];
ListLinePlot[Table[{x[[i]], y[[i]]}, {i, Length[x]} ]]

Here is one of the many possible plots generated by the above

This plot was generated with 100,000 points and some of them may not even be visible since other points overlap part of the point. The question is, is there some built in function that gives us a way to remove the redundancies and thus reducing the size of the figure?

I should be more precise. If this plot is to be exported as a raster plot (jpeg, png) then there is no problem. The problem arises if we export it as a pdf or eps. The resulting vector format will have to plot all of those points and it will take a while to display it when in reality we only have a few points to display.

The very naive way I did before was to remove some points from my initial data, as follows:

ListLinePlot[Table[{x[[i]], y[[i]]}, {i, 1, Length[x], 100}]]

And you can imagine what will happen if I skip 1000 points at a time

ListLinePlot[Table[{x[[i]], y[[i]]}, {i, 1, Length[x], 1000}]]

Possible Solution 1

The first idea that came to my mind was to loop over my data and see if the points are close to each other, if they are sufficiently close (we can use a parameter to decide this) then we remove one of them.

This solution will remove points from our initial data set but we want to make it so that we do not remove important features of the plot as I did with the 3rd plot.

Possible Solution 2

The other solution was to export/rasterize the final figure I wanted (only the lines in the plot, this excludes the axes or any other objects). Once this is done then we outline the plot and make a vector figure of it again.

Keep in mind that we are trying to reduce the size of points that are being plotted while making the plot resemble what Mathematica plotted very naively.

Any ideas of how to implement any of the ideas above or an even better one to achieve the desired result?

Answer 1

This filters out your points by their EuclideanDistance. I think it makes a pretty good job preserving the curve's features with very few points:

rx[n_] :=  Accumulate[Prepend[RandomVariate[ExponentialDistribution[1000], n], 0]]
ry[n_] :=  Accumulate[Prepend[RandomVariate[NormalDistribution[0, .001], n], 0]]
c = Transpose[{rx@#, ry@#}] &@100000;
t = Table[(lp = {c[[1]]}; 
          If [EuclideanDistance[lp[[-1]], #] > parm, AppendTo[lp, #]] & /@ c;
          Column[{StringJoin["Points ", ToString@Length@lp], 
          ListLinePlot@lp}]), 
    {parm, 10^-2, 5 10^-2, 10^-2}];

Edit

Here you may see how it works at a "microscopic" scale. Coarser iterations chose points not selected in previous ones. That's how it preserves "the form".

Answer 2

Here is a modified version of belisarius' answer I ended up using. Instead of using the euclidean distance to find out if two points are close to each other I set a fixed rectangle in an area and collapse it into a line. For instance, consider the same example.

rx[n_] := Accumulate[Prepend[RandomVariate[ExponentialDistribution[1000], n], 0]]
ry[n_] := Accumulate[Prepend[RandomVariate[NormalDistribution[0, .001], n], 0]]
c = Transpose[{rx@#, ry@#}] &@100000;
ListLinePlot[c]

Here we can see that from 0 to 100 there are about 325 pixes which means that any point that is located in an interval of length 100/325 essentially lies in the same vertical line. There may be many points in a rectangle of this length so we need to find the minimum and maximum so that we can draw a line. Here is a function I made to filter the data given the length of PlotRange/pixels.

filter[c_, uup_] := Module[{lp, min, max, x},
  lp = {};
  min = c[[1]];
  max = c[[1]];
  x = c[[1, 1]];
  If[#[[1]] - x > uup,
     If[min[[1]] < max[[1]],
      AppendTo[lp, min];
      AppendTo[lp, max],
      AppendTo[lp, max];
      AppendTo[lp, min]
      ];
     min = #;
     max = #;
     x = #[[1]],
     If[#[[2]] > max[[2]], max = #];
     If[#[[2]] < min[[2]], min = #];
     ] & /@ c;
  Return[lp];
  ]

Using said value we obtain:

lp = filter[c, 100/325];
Print["Points: " <> ToString@Length@lp];
ListLinePlot[lp]

The amount of overlapping has been reduced and thus we do not see the bold look it had. Here is a side by side comparison with different parameters for the filter.

Column@Table[(
   lp = filter[c, 100/i];
   Grid[{
     {
      StringJoin["Points ", ToString@Length@lp], 
      StringJoin["Points ", ToString@Length@c]
      }, {
      ListLinePlot[lp, ImageSize -> 2.5 72],
      ListLinePlot[c, ImageSize -> 2.5 72]
      }
     }]),
  {i, 100, 1000, 200}]

Edit

I'm adding a modified version of the filter written by Thomas:

filter1[c_, uup_] := Module[{
   d = Split[c, Floor[First@#1, uup] === Floor[First@#2, uup] &],
   x, ord, y
   },
  y = (#[[2]] & /@ (Transpose /@ d));
  ord = Ordering /@ y;
  ord = {First@#, Last@#} & /@ ord;
  ord = Sort /@ ord;
  Join @@ MapThread[Part, {d, ord}]
  ]

Edit 2

This version is again much faster. It is basically a different way of calculating the vector d, the rest is identical to filter1 above:

filter2[c_, uup_] := Module[{tmp, len, d, ord, y},
  tmp = Split[Floor[c[[All, 1]], uup]];
  len = Accumulate@Flatten[{1, Length /@ tmp}];
  d = MapThread[Take[c, {#1, #2}] &, {Most@len, Rest@len - 1}];
  y = #[[All, 2]] & /@ d;
  ord = Ordering /@ y;
  ord = {First@#, Last@#} & /@ ord;
  ord = Sort /@ ord;
  Join @@ MapThread[Part, {d, ord}]
  ]

Answer 3

To reduce the number of points, I'd use Interpolation:

f = Interpolation[Transpose@{x,y}, InterpolationOrder->1];
Plot[f@x,{x,0,100}]

Answer 4

As I said in my comment, the Douglas-Peucker algorithm is a standard method for simplifying polygonal curves. This demonstration by Mark McClure contains a Mathematica implementation. Since he is a user on this site, I feel like he deserves the points for this answer. Nevertheless, in the interests of having an answer here, and hoping this falls under fair use, I quote his implementation:

dist[q : {x_, y_}, {p1 : {x1_, y1_}, p2 : {x2_, y2_}}] := With[
   {u = (q - p1).(p2 - p1)/(p2 - p1).(p2 - p1)},
   Which[
    u <= 0, Norm[q - p1],
    u >= 1, Norm[q - p2],
    True, Norm[q - (p1 + u (p2 - p1))]
    ]
   ];
testSeg[seg[points_List], tol_] := Module[
    {},
    dists = dist[#, {points[[1]], points[[-1]]}] & /@ 
      points[[Range[2, Length[points] - 1]]];
    max = Max[dists];
    If[max > tol,
     pos = Position[dists, max][[1, 1]] + 1;
     {seg[points[[Range[1, pos]]]], 
      seg[points[[Range[pos, Length[points]]]]]},
     seg[points, done]]] /; Length[points] > 2;
testSeg[seg[points_List], tol_] := seg[points, done];
testSeg[seg[points_List, done], tol_] := seg[points, done];
dpSimp[points_, tol_] := 
  Append[First /@ First /@ Flatten[{seg[points]} //. 
       s_seg :> testSeg[s, tol]], Last[points]];

(P.S. @Mark, if you want to add your own answer, I'll be happy to delete this one.)

As suggested by @whuber, for plots it's more meaningful to compare vertical differences than geometric distances. This just requires changing the dist function to

dist[q : {x_, y_}, {p1 : {x1_, y1_}, p2 : {x2_, y2_}}] := 
  With[{u = (x - x1)/(x2 - x1)}, Abs[y - (y1 + u (y2 - y1))]];

Then you can do, for example,

ListLinePlot[dpSimp[Transpose@{x,y}, 0.01]]

to get a plot that differs by at most 0.01 from the original data. Here are four plots with tolerances I picked to give just under 1000, 500, 200, and 100 points respectively.