How to implement the sample-point process like the built-ins of Mathematica?

Question

Consider this:

pts = {{0, 0}, {1, 1}, {2, -1}, {3, 0}, {4, -2}, {5, 1}};
f = BSplineFunction[pts]

I can use ParametricPlot to visualize this B-spline curve:

Show[
  Graphics[{Red, Point[pts], Green, Line[pts]}, Axes -> True], 
  ParametricPlot[f[t], {t, 0, 1}]]

points = {{0, 0}, {1, 1}, {2, -1}, {3, 0}, {4, -2}, {-5, 1}};
g = BSplineFunction[points];
Show[Graphics[{Red, Point[pts], Green, Line[pts]}, Axes -> True], 
 ParametricPlot[g[t], {t, 0, 1}, AspectRatio -> Automatic]]

But when I sample by hand, I will do the following operation:

curvePts = f /@ Range[0, 1, .01];
ListPlot[curvePts]

However, when I double-click the first graph, I discovered that they are different:

In addition, I notice that

g = Sin[#] &;
{ListPlot[g /@ Range[0, 10, .1]], Plot[Sin[x], {x, 0, 10}]}

Question

• How do I sample points like Mathematica does, according to the steepness of the curve?

In this answer, I used Uniform Sampling Method.

Answer 1

Plot uses two different algorithms depending on whether PerformanceGoal is set to Quality or Speed. Yaroslav Bulatov wrote here, i.e. in the link provided by Szabolcs in a comment above, that:

Plot starts with 50 equally spaced points and then inserts extra points in up to MaxRecursion stages... According to Stan Wagon's Mathematica book, Plot decides whether to add an extra point halfway between two consecutive points if the angle between two new line segments would be more than 5 degrees.

It turns out that this corresponds to the algorithm used with PerformanceGoal -> "Speed". Remember to set the MaxRecursion option as well to compare with the plots below. In the third edition, the section on adaptive plotting in Stan Wagon's Mathematica in Action can be found on page 28.

One possible implementation of this algorithm is this:

addPoint[f_][{x1_, x2_}] := Module[{midPoint, v1, v2},
  midPoint = (x1 + x2)/2;
  v1 = {x1, f[x1]} - {midPoint, f[midPoint]};
  v2 = {midPoint, f[midPoint]} - {x2, f[x2]};

  If[VectorAngle[v1, v2] > 5 Degree, Unevaluated@Sequence[x1, midPoint], x1]]

addPoints[f_][pts_] := Append[Developer`PartitionMap[addPoint[f], pts, 2, 1], Last@pts]

addPoints[f_][pts_] takes a list of x values and a function f and adds more x values to the list according to the criteria mentioned by Yaroslav.

In order to test the algorithm we can do this:

plotPoints = 50;
maxRecursions = 4;
{min, max} = {0, 10 Pi};
initialPts = N@Table[x, {x, min, max, (max - min)/(plotPoints - 1)}];

(* Find the points corresponding to the Sin function *)
steps = NestList[addPoints[Sin], initialPts, maxRecursions];

(* Visualization: *)
visualizePts[f_, {min_, max_}][pts_] := Plot[
  f[x], {x, min, max},
  Mesh -> {Thread[{pts, Directive[Red, PointSize[Medium]]}]},
  ImageSize -> 300
  ]

Partition[visualizePts[Sin, {min, max}] /@ steps, 2] // Grid

Answer 2

Wanna listen to a story? :)

It was around 2002 when I finally became fed up with ParametricPlot3D[] and its inability to adaptively plot space curves. Recall that this was the old Graphics[] system where all the pictures were effectively done in PostScript. Thus, I set out to look for a way to adaptively plot curves in general.

I was at the time very familiar with the trio of MaxBend, PlotDivision, and PlotPoints, and set out to look for an algorithm that used these quantities. Eventually, I came up with a solution that used Nest[]. Nevertheless, I found my solution to be a bit clunky in places, so I hesitated on putting it up to MathSource. But I continued to use my function privately in place of ParametricPlot3D[], using it both for scientific and artistic purposes.

Around 2006, while very idly browsing through the built-in Mathematica packages, I finally chanced upon the source of Graphics`Spline`​. Apparently, an adaptive plotting algorithm was sitting right there all along! I then re-wrote my custom ParametricPlot3D[] to use modified versions of these internal functions (as they were not made general enough to handle the 3D case). I was still thinking about submitting my custom routine to MathSource, but 1.) my modification to handle the 3D case felt a bit too simple, and 2.) I was quite uncomfortable with the recursive implementation (it involved Fold[]-ing a function through a list of points, where that function being folded was then defined in terms of the original calling function). So I held off, still poking and prodding at the routine to see if there was a way to write the method without recursive calls.

Finally, I did figure out a way to implement the function using ReplaceRepeated[] (//.), and after some testing I felt confident that I could submit to MathSource. I prepared documentation and examples meticulously, bundled it all up in a nice *.zip file, and was just about ready to call it a night, with the plan to submit the thing tomorrow morning…

…and then Mathematica 6 was released the next day. A quick look at the capabilities of the brand new graphics system told me at once that my package has been rendered obsolete.

Oh well.

---

The point of that rather long-winded story above is that now, in 2015, I was able to find some of my old adaptive plotting code during a look-through of a rarely-used e-mail address. It would seem that I had sent it to a friend who was also dissatisfied with ParametricPlot3D[]. The timestamp indicates I had sent it on August 2006. Ah, memories…

Thus, consider this post to be a very much belated release of early 2000's code. ;)

I will omit the other fancy stuff and wrapper routines behind my custom ParametricPlot3D[], and just present the engine of the adaptive algorithm. Here it is:

bendAngle[{u1_, pt1_?VectorQ}, {u2_, pt2_?VectorQ}, {u3_, pt3_?VectorQ}, min_] := 
    Module[{norm = Sqrt[#.#] &, v1 = pt1 - pt2, v2 = pt2 - pt3, n1, n2}, 
           If[Max[u3 - u2, u2 - u1] <= min, 0,
              n1 = norm[v1]; n2 = norm[v2];
              2 ArcTan[norm[v1 n2 + n1 v2], norm[v1 n2 - n1 v2]]]]

insertPoints[pt1 : {u1_, _}, pt2 : {u2_, _}, pt3 : {u3_, _}, fun_] := 
      Module[{i1 = u1 + (u2 - u1)/2, i2 = u2 + (u3 - u2)/2},
             {pt1, {i1, fun[i1]}, pt2, {i2, fun[i2]}, pt3}]

plotAdaptive[fun_, {tmin_, tmax_}, maxb_?NumericQ, npts_Integer, div_Integer] := 
    Block[{t}, 
          Table[N[{t, fun[t]}], {t, tmin, tmax, (tmax - tmin)/(npts - 1)}] //. 
          p : {u___, a : {_?NumericQ, _?VectorQ}, b : {_?NumericQ, _?VectorQ},
                     c : {_?NumericQ, _?VectorQ}, v___} /; 
          bendAngle[a, b, c, (tmax - tmin)/(npts div)] > maxb Degree :> 
          Join[{u}, insertPoints[a, b, c, fun], {v}]]

If the code seems clunky in places, I apologize; I will readily admit that my coding ability at the time was not too good. :) One tic of note is that I had never bothered to use Norm[], as I had wanted to maintain compatibility with versions older than version 5. (Also. the code did begin to take shape during the time of version 4.)

Before demonstrating the code, let's set the mood, and step through a time machine…

<<Version5`Graphics`

Ah, PostScript. I can't say I miss it a lot, but the old pictures had a charming quality to them. :) Here's the usual sine example:

sinx[x_] := {x, Sin[x]};
With[{pts = plotAdaptive[sinx, {0, 2 Pi}, 10., 25, 4][[All, -1]]}, 
     Show[Graphics[{Line[pts], AbsolutePointSize[3], Point /@ pts}]]]

Here's a clothoid:

clothoid[t_] := {FresnelS[t], FresnelC[t]};
With[{pts = plotAdaptive[clothoid, {-5, 5}, 10., 25, 4][[All, -1]]}, 
     Show[Graphics[{Line[pts], AbsolutePointSize[3], Point /@ pts}, 
          AspectRatio -> Automatic]]]

And finally, one of the reasons why I had written those routines, the curve of Archytas:

archytas[t_] := {Cos[t]^4/(1 + Sin[t]^2)^2,
                 (Cos[t] Sin[2 t])/(1 + Sin[t]^2)^2,
                 (Sqrt[2] Cos[t] Sin[t])/(1 + Sin[t]^2)};
With[{pts = plotAdaptive[archytas, {-Pi, Pi}, 10., 25, 4][[All, -1]]}, 
     Graphics3D[{Line[pts], AbsolutePointSize[3], Point /@ pts}, 
                Boxed -> False] // Show]

It still works in 2015. I'm touched.

Thanks for reading my long post, if you got this far! Oops, I forgot. Let's go back to modern times:

<<Version6`Graphics`

Answer 3

Earlier in the summer I had written the following for How to obtain adaptive sampling as in Plot function?. It is something like J. M.'s technique. But instead of a new version of Mathematica coming along, my brain turned on and I discovered a workaround using FunctionInterpolation. I've been considering posting the following as an answer to that question, but perhaps it fits here better.

As I remarked in a comment to the above-mentioned question, how to adapt sampling depends on the goal. The goal is usually an approximation of some sort. It might be a curve, a function, an integral, and so forth. An approximation implies an error. The adaptation therefore should be made with a view to reducing the error estimation. There are different ways to look at error. The following subdivides each subinterval until the error estimate for each subinterval meets the goal specified by PrecisionGoal and AccuracyGoal.

The idea is fairly simple. An "active" subinterval has a special (undefined) head step. At each recursive step, the head is replaced by split, which decides whether the interval meets the error tolerance (... /; err[seg] < ag + pg Abs[fc]) or should be split in two. (One could split into ten or some other number, I suppose, but splitting into two is easy.) Intervals with head step keep getting split until none are left, using ReplaceAll like J. M.

Some remarks: (1) The use of Experimental`CreateNumericalFunction was to handle WorkingPrecision, as well as other potential alternatives I was exploring. It is not necessary here. It could be replaced by Function, but since it was already in place, I left it. (2) I used Plot to get an initial sampling. Plot has been refined over the years, and I trust it. I added the initial point of the interval (Plot samples very close to the initial point but not at the point), because the purpose was to approximate a function over the closed interval. (3) There is an option "ErrorNorm", with which an error measure can be passed to the function mesh. The error norm will be passed a list of three points representing the interval with its midpoint.

ClearAll[mesh, mesh`imesh];
Options[mesh] = {PrecisionGoal -> 6, AccuracyGoal -> 6, 
   WorkingPrecision -> MachinePrecision, "InitialPoints" -> Automatic,
   "ErrorNorm" -> Automatic};
mesh[f_, {var_, a_, b_}, opts : OptionsPattern[]] := mesh`imesh[f, {var, a, b}, opts];

Begin["mesh`"];
ClearAll[mesh`step, mesh`split, mesh`err, mesh`err0, mesh`imesh, mesh`initialpoints];

mesh`initialpoints = Flatten[{0., 
    Cases[Plot[1, {x, 0, 1}, PlotPoints -> 11, MaxRecursion -> 0], 
     Line[p_] :> p[[All, 1]], Infinity]}];

Options[imesh] = Options[mesh];

(* default: overestimate Δf at midpoint from curve to chord *)
err0[{{_, fa_}, {_, fc_}, {_, fb_}}] := Abs[fc - (fa + fb)/2]/4.;

imesh[f_, {var_, a0_, b0_}, OptionsPattern[]] := 
  Block[{step, split, err,
      nf = Experimental`CreateNumericalFunction[{var}, f, {}, 
        WorkingPrecision -> OptionValue@WorkingPrecision]},
   err = OptionValue["ErrorNorm"] /. Automatic -> err0;
   With[{pg = 10.^-OptionValue[PrecisionGoal],
         ag = 10.^-OptionValue[AccuracyGoal]},
    step[seg : {{a_, fa_}, {c_, fc_}, {b_, fb_}}] /; err[seg] < ag + pg Abs[fc] :=
      Sequence[{a, fa}, {c, fc}];
    step[{{a_, fa_}, {c_, fc_}, {b_, fb_}}] := Sequence[
      split[{{a, fa}, {#, nf[{#}]} &[(a + c)/2], {c, fc}}], 
      split[{{c, fc}, {#, nf[{#}]} &[(b + c)/2], {b, fb}}]];
    ];
   Append[#, {N[b0, OptionValue[WorkingPrecision]], nf[{b0}]}] &[
    split[{{#1, nf[{#1}]}, {(#1 + #2)/2, nf[{(#1 + #2)/2}]}, {#2, nf[{#2}]}}] & @@@
      Partition[
       SetPrecision[OptionValue["InitialPoints"] /. 
         Automatic -> Rescale[initialpoints, {0, 1}, {a0, b0}], 
        OptionValue[WorkingPrecision]],
       2, 1] //. split -> step]
   ];

End[];

Examples:

Linear interpolation. The default method tries to make the error of a linear interpolation of the points meet the precision and accuracy goals (with a rather crude but fast error estimate).

mesh[x^2, {x, -1, 4}, PrecisionGoal -> 2, AccuracyGoal -> 2]
(*
{{-1., 1.}, {-1., 1.}, {-1., 0.999999}, {-0.759549, 0.576914},...,{3.74217, 
  14.0039}, {4., 16.}}
*)

Plot[x^2, {x, -1, 4}, 
 Epilog -> {Red, Point[mesh[x^2, {x, -1, 4}, PrecisionGoal -> 2, AccuracyGoal -> 2]]}]

Subdivide according to the angle change. In this case the angle is not a relative error, so the precision goal should be set to Infinity; in such a case, AccuracyGoal determines when subdivision halts. The VectorAngle of the Differences between the three consecutive points serves as the error function. An AccuracyGoal of 1 corresponds to a bend of (at most) 0.1 radians or a little less than 6 degrees.

Plot[x^2, {x, -1, 4}, 
 Epilog -> {Red, 
   Point[mesh[x^2, {x, -1, 4}, PrecisionGoal -> Infinity, AccuracyGoal -> 1, 
     "ErrorNorm" -> Function[x, Abs[VectorAngle @@ Differences[x]]]]]}]

Note: One difference (I think) between this and J. M.'s method, is that subdivision is done per subinterval: If $[a, b]$ and $[b,c]$ are divided at the midpoints into sequences {a,d,b} and {b,e,c}, the sequence {d,b,c} will not be examined. This is similar in a way to the meshing of Plot3D, ContourPlot and DensityPlot, but the real reason is that the algorithm was developed to minimize the error in approximating a function, not a curve. Whether a method is better often depends on what application one has in mind.

Answer 4

Being highly astonished about all the very interesting answers, I only want to add an approach using NMinimize. Not very fast, but straightforward.

The idea behind is to minimize the squared cross product of neighboring normalized line-segments, which is equal to Sin[angle between lignsegments]^2.

Given an arbitrary function g[t], 0<t<10 (for example g = Sin[#] &; ), in the first step the inner points of the polygon have to be defined

n = 100; (* number of points *)
t0 = 0;
t1 = 10;
punkte = Join[{{t0, g[t0]}},Table[{t[i], g[t[i]]}, {i, 1, n - 2}] , {{t1, g[t1]}}]; (* polygon *)
ds = Rest[punkte] - Most[punkte]; (* line segments*)
ds2 = Map[Sqrt[#.#]/2 &, Rest[ds]] + Map[Sqrt[#.#]/2 &, Most[ds]]; (* weights knod i*)
Δα =MapThread[(1 -   (#1.#2)^2/( #1.#1 #2.#2 )) & ,{Rest[ds],Most[ds]}]; (* squared crossproduct *)

Minimization:

ti = punkte[[2 ;; -2, 1]];(* unknown parameters t[i]*) 
opt = Quiet@FindMinimum[{    Δα.ds2  }, 
Table[ {t[i], t0 + (t1 - t0) i/(n - 1)}, {i, 1, n - 2}] (* initial points*)];
(* FindMinimum gives local minimum *)

Show [{Plot[g[t], {t, t0, t1}],ListPlot[punkte /. opt[[2]], PlotStyle -> Red]}]

Checking the angle between adjacent line elements

Δαi =MapThread[ArcCos[#1.#2/Sqrt[#1.#1 #2.#2]] & , {Rest[ds ], Most[ds ]} /. opt[[2]]];
ListPlot[Δαi/Degree,AxesLabel -> {index, "Degree"}] 

shows angle values smaller than 3 Degrees which conforms well with the Plot option PerformanceGoal -> "Speed".

Further improvement could be achieved by changing the number of points n.