Programming Tips

Mathematica in Education
Vol. 3 No.4
Copyright 1994
TELOS/Springer-Verlag

Getting Inside Plots

How to get a closer look at the points that make up a smooth-looking graph, with an application to paths in the plane.

by Stan Wagon
Macalester College
wagon@macalstr.edu

The Points Underlying a Plot

Mathematica's Plot and ParametricPlot commands work adaptively: they adapt to the curvature of a plot and compute more data points where the curviness demands it. These points are connected by straight lines to create an illusion of smoothness. The options PlotPoints and MaxBend control the resolution: PlotPoints tells how many subintervals to use initially (default is 25) and MaxBend tells what angle between segments is acceptable (default is 10û). At times, one would like to have access to the data points that define the plot. Here is one way to get them.

p = Plot[Sin[(3.3 - x)^2], {x, 0.8, 3.5}];
InputForm[p]

Out[306]//InputForm=
     Graphics[{{Line[{{0.8, -0.03317921654755681682}, 
      {0.85625, -0.3062690386043116985},{0.9125, -0.5505551120033017266}, 
      {0.96875, -0.7502623106869461109}, {0.996875, -0.8298286478357260429}, 
      {1.025, -0.8946111836030903104}, {1.053125, -0.9440618492831936536}, 
      {1.0671875, -0.9629494645818835751}, {1.08125, -0.977934380265551051}, 
     << plus 100 more pairs of reals >>    }]}}, 
  {PlotRange -> Automatic, AspectRatio -> GoldenRatio^(-1), 
    << plus 20 additional option settings >>  }]

The actual points are clearly embedded in this object. A little searching shows that they lie in p[[1, 1, 1, 1]].

Take[p[[1, 1, 1, 1]], 5]

Now we can do several things. For example, we can redo the plot with the points superimposed, to get a sense of the adaptive technique (Figure 1).

p = Plot[Sin[(3.3 - x)^2], {x, .8, 3.5},
		DisplayFunction->Identity];
Show[p,
     Graphics[{AbsolutePointSize[3],Map[Point, 
                                       p[[1,1,1,1]]]}],
     DisplayFunction -> $DisplayFunction];

Figure 1:  The adaptive nature of the plotting algorithm becomes evident when a dot is placed at each computed point on the graph.

As promised, the points pile up where the graph curves the most. Warning: If the plot has some options (e.g., PlotStyle->{GrayLevel[0.5]}), then the points live in p[[1, 2, 1, 1]]; thus if you are incorporating this into a routine that might see some options, use p[[1, -1, 1, 1]]. Note further that this technique gives us the points in sorted order. If you would like to make a movie that shows the plotting algorithm in action point by point, you will want the points in the order that Plot generates them. This can be done by defining the function to be plotted with a side condition that stores the points.

points = {};
f[x_] := Module[{y = Sin[(3.3 - x)^2]}, AppendTo[points, y]; y]
Plot[Sin[(3.3 - x)^2], {x, 0.8, 3.5}]

The points now lie in points, and this can be used to animate the plotting algorithm (details left as an exercise; use a Do-loop and Take[points, i]).

Coloring Parametric Plots

We mention another application of this sort of detective work. ParametricPlot[f[t], {t, a, b}] does not allow for the coloring of the segments according to the parameter t. But this can be a nice thing to do, since it may be difficult for the eye to follow an intricate black curve that crosses itself several times. If the curve is made of changing colors, it is easier to chase it around as the parameter increases. This very point was discussed in Matthew Thomas's column in the winter 1994 issue of this journal (p. 39). It was shown there how the coloring functions available to ParametricPlot3D allow such a colored plot to be obtained. But this 3-dimensional workaround is slow and very memory intensive. The approach we present here, which uses only 2-dimensional plotting, is much faster, and maintains the adaptive nature of ParametricPlot (ParametricPlot3D is not adaptive, but instead examines equally spaced parameter values).

Here we cannot use the p[[1, -1, 1, 1]] trick because the parameter values would not appear in the raw list of data points. So we use the AppendTo[points, ...] idea. The problem is that AppendTo is very slow. This is not a serious problem for Plot, since the number of points is usually less than about 100. But a complicated parametric plot might require 1000 or more points. To avoid AppendTo, we generate an initial table of 1500 Nulls (we assume that no more than 1500 points will be needed; the Nulls use less memory than, say, 1500 zeros), and then change table entries to pairs {t, {x[t], y[t]}} as the plot proceeds via the auxiliary function g. This necessitates keeping a counter i; we use ++i, which preincrements the value of i by 1. The switch from AppendTo to the explicit table method cuts the time down by a factor of about 2 or more.

The code below uses Sort to order the data points according to t, then Partition to transform the collection of points {t, {x, y}} into a set of pairs:  {{t, {x, y}}, {nextt, {nextx, nexty}}}. Then /. is used to collapse each such pair to {Hue[(t-a)/(b-a)], Line[{{x, y}, {nextx, nexty}}]}. Showing these graphics elements yields a sequence of colored segments whose colors change continuously with the parameter.

ColorParametricPlot[f_, {t_, a_, b_},  opts___] := Module[
    {data = Table[Null, {1500}], i = 0, g, fval},
  g[s_] := (fval = N[f /. t->s]; data[[++i]] = {s, fval}; fval);
  ParametricPlot[g[t], {t, a, b}, opts, Compiled->False,
         DisplayFunction->Identity];
  data = Partition[Sort[Take[data, i]], 2, 1];  
  Show[Graphics[data /. {{s_, pt_List}, next_} :> 
                 {Hue[(s-a)/(b-a)], Line[{pt, next[[2]]}]}], 
   PlotRange->All, AspectRatio->Automatic, Background->GrayLevel[0]]]

The black background is added to better show the colors. You will need a color monitor to see the glittering effects these plots can have. But using shades of gray is an acceptable alternative. The Gray->True option to ColorParametricPlot (the fuller electonic version, not the version printed here) replaces Hue[] with GrayLevel[0.8 * Sin[Pi (s-a)/(b-a)]] and eliminates the black background, thus producing a good-looking curve in continuously changing shades of gray. The example in Figure 2 compares the raw black-and-white plot, on which it is difficult to trace the actual path, with a gray-scale plot, which shows that the potential small loops near the center are not really loops at all, but more like U-turns. Figure 2 was generated by

ParametricPlot[{Cos[4 t] Cos[5 t], Sin[t] Cos[3 t]}, 
	{t, 0, 2 Pi},
	PlotPoints->200, Axes->None]

and

ColorParametricPlot[{Cos[4 t] Cos[5 t], Sin[t] Cos[3 t]},
	{t, 0, 2 Pi},
    PlotPoints->200, Gray->True, 
    PlotStyle->AbsoluteThickness[2]];

Figure 2:  A comparison of a vanilla parametric plot with one whose shading changes continuously as the parameter varies. The shaded graph clarifies the curve's behavior in the vicinity of the small loops near the center.

Seeing Speed

A useful modification is to add an option that colors the curve according to the speed of the particle. Since the speed is easily computed as norm[f'[t]] where norm[v_] := Sqrt[v.v], this is straightforward, and the electronic version of ColorParametricPlot implements such a Speed->True option. Such a coloring gives a sense of where the motion is fast and where it is slow. We use lighter shades (or red) to indicate greater speed, and darker shades (or blue) for less. The code below produces the examples in Figure 3.

ColorParametricPlot[Sin[4 t] * {Sin[t], Cos[t]}, {t, 0, 2 Pi},
  Speed->True, PlotStyle->AbsoluteThickness[4]];

ColorParametricPlot[{t - 1.3 Sin[t], 1 - 1.3 Cos[t]}, {t, -Pi/2, 5 Pi},
  Speed->True, PlotStyle->AbsoluteThickness[4],
  PlotRange->{-.5, 2.5}];

Figure 3:  When the standard parametric plots of an 8-leafed rose or a cycloidal type of curve are colored according to speed, a curious phenomenon comes to the fore: the particle moves slowest when the curves are tight and fastest when the curvature is small.

It seems as if the standard parametrizations of familiar curves are fast and slow according to the curvature of the curve. Since there is a simple formula for radius of curvature Ñ (x¢(t)2 + y¢(t)2)3/2 / |x¢(t)y¢¢(t) Ð x¢¢(t)y¢(t)| Ñ we can easily investigate the speed/curvature relationship by using one-dimensional plots as follows (see Figure 4).

f[t_] := Sin[4t] {Sin[t], Cos[t]}
norm[v_] := Sqrt[v . v]
speed[t_] := Evaluate[norm[f'[t]]]

curvature[t_] := 
  Evaluate[speed[t]^3 / Abs[(f'[t]*{1,-1}) . Reverse[f''[t]]] ]

Plot[{speed[t], curvature[t]}, {t, 0, 2 Pi}, 
  PlotStyle->{AbsoluteThickness[2], AbsoluteThickness[1]}
];

Figure 4:  A comparison of the speed (thicker curve) and curvature of a standard rose. They rise and fall together, which explains why the rose's speed seems natural (i.e., it slows down at the tight curves).

We leave further investigations along these lines to the reader. A good exercise for students is to come up with a parametrized curve whose motion is contrary to what appears natural: it should be fastest when the curves are tight and slowest when they are straight; an example is shown in Figure 5.

Figure 5:  A contrived parametrized plot of a cycloid-like curve in which the speed is slowest (shading is darkest) on the straighter segments and fastest around the tight curves.

Color examples

ColorParametricPlot[{Cos[4 t] Cos[5 t], Sin[t] Cos[3 t]},
 {t, 0, 2 Pi}, PlotPoints->100, MaxBend->5, 
 Background->GrayLevel[0], Frame->True,
 FrameTicks->None, 
 FrameStyle->{AbsoluteThickness[4], Hue[58/65]}];

ColorParametricPlot[{Cos[3 t], Sin[t] Cos[5 t]},
 {t, 0, 2 Pi}, PlotPoints->25, MaxBend->5, 
 Background->GrayLevel[0], Frame->True,
 FrameTicks->None, 
 FrameStyle->{AbsoluteThickness[4], Hue[58/65]}];

ColorParametricPlot[Sin[4t]{Sin[t], Cos[2t]}, {t,0,2Pi},
PlotPoints->35, Gray->!True, Background->GrayLevel[0.5],
PlotStyle->AbsoluteThickness[5]];

Here is a neat one that shows how the speed works: fastest on the straight segments!

ColorParametricPlot[{Cos[3 t], Sin[t] Cos[5 t]},
 {t, 0, 2 Pi}, PlotPoints->100, MaxBend->5, 
 Background->GrayLevel[0], Frame->True, Speed->True,
 FrameTicks->None, PlotStyle->AbsoluteThickness[4],
 FrameStyle->{AbsoluteThickness[4], Hue[58/65]}];

Full version of ColorParametricPlot with lots of options.

Clear[ColorParametricPlot];

ColorParametricPlot::usage = "ColorParametricPlot[f, {t, a, b}] generates a parametric plot whose color changes continuously as a
function of t, thus making it easier for the eye to trace out the path. The usual options to ParametricPlot may be attached. The Speed option allows the coloring to be a function of the speed of motion instead of the parameter.";

Gray::usage = "Gray is an option that, when set to True, causes the colors to be replaced by shades of gray. Default is False.";

Speed::usage = "Speed is an option that, when True, causes the coloring to run from blue (slowest) to red (fastest) corresponding to the speed of the particle. Default is False. If Gray and Speed are both True, then the shading goes from black (slowest) to light gray (fastest).";

Off[General::spell1];

Options[ColorParametricPlot] = {Gray->False, Speed->False, 
     PlotStyle->{AbsoluteThickness[1]}};

ColorParametricPlot[f_, {t_, a_, b_},  opts___] := Module[{i = 0, g, fval,
  gr, sp, data = Table[Null, {1500}], norm, speedMin, speedMax, styles, options,
  options1, c1, c2, func1, func2, aa, bb, sFunc, d},
 
 {gr, styles, sp} = {Gray, PlotStyle, Speed} /. {opts} /.
          Options[ColorParametricPlot];
 
 options = Join[Select[{opts}, MemberQ[First /@ Options[ParametricPlot], First[#]]&], 
     {Compiled->False, DisplayFunction->Identity}];
 
 options1 = Join[Select[{opts}, MemberQ[First /@ Options[Graphics], First[#]]&], 
     {PlotRange->All, AspectRatio->Automatic}];
 
 norm[v_] := N[Sqrt[v.v]];
 
 If[sp,  speed[tt_] := Evaluate[norm[D[f, t] /. t->tt]];
         {speedMin, speedMax} = (ll=Last[Transpose[Plot[speed[t], {t, a, b}, 
      PlotRange->All, DisplayFunction->Identity][[1,1,1,1]]]]; 
      {Min[ll], Max[ll]})];

g[s_] := (fval = N[f /. t->s]; data[[++i]] = {s, fval}; fval);
  ParametricPlot[g[t], {t, a, b}, Evaluate[options]];
  data = Partition[Sort[Take[data, i]], 2, 1];  
 
{c1, c2, func1, func2, aa, bb, sFunc, d} = Which[
    gr && !sp, 
    		{0.8, N[Pi], GrayLevel[Min[Max[#,0], 1]]&, Sin, a,  	b,  Identity, 0}, 
   !gr && !sp, 
   			{1,   1,     Hue,  									     					Identity, a,   b,  Identity, 0},
   !gr &&  sp,
   		 {0.5, 1,     Hue,  			     	Identity, speedMin, speedMax, speed, 0.5},
    gr && sp,
      {0.8, N[Pi/2], GrayLevel[Min[Max[#,0], 1]]&, Sin, speedMin, speedMax, speed, 0}];

 Show[Graphics[{styles, data /. {{s_, pt_List}, next_} :> 
     {func1[d + c1 * func2[N[c2 * (sFunc[s]-aa)/(bb-aa)]]], 
            Line[{pt, next[[2]]}]}}],  Evaluate[options1]  ]]

Code for contrived function in Figure 5.

d=7;
g[t_] := (((2 t-2Pi)^d)/(2 Pi)^(d-1) + 2 Pi)/2

f[t_] := {t-1.5 Sin[t], 1 - 1.5 Cos[t]};

ColorParametricPlot[Evaluate[f[g[t]]],
   {t, 0, 2 Pi},
   Speed->True, PlotRange->All, Gray->True,
   PlotStyle->AbsoluteThickness[3], Axes->None];

The slow way, for timing comparison alluded to in text;

Clear[f];
points = {};
f[t_] := Module[{fval = Sin[4 t]{Cos[t], Sin[t]}}, 
           AppendTo[points, {t, fval}]; fval]

Print[First[Timing[ParametricPlot[f[t], {t, 0, 2 Pi}, 
DisplayFunction->Identity, Compiled->False,
AspectRatio->Automatic]]]];

 ends = Partition[Sort[points], 2, 1];
 coloredPairs = ends /. {{s_, {x_, y_}},
         {s1_, {x1_, y1_}}} :> {s, {{x, y}, {x1, y1}}};