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

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:

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.

• If I understand you wish to extract the points from the plot, or at least doing so would be a solution to your question? If that is the case this is a duplicate. If not please explain further. – Mr.Wizard Jul 15 '15 at 9:20
• Seen this? – J. M.'s discontentment Jul 15 '15 at 9:21
• Mathematica uses adaptive sampling. It uses denser sampling where the curvature is higher, to achieve good accuracy without having to sample too many points. Mesh -> All will show the points. MaxRecursion controls the number of refinement steps. MaxRecursion -> 0 will uniformly sample PlotPoints number of points. – Szabolcs Jul 15 '15 at 9:59
• Read here. @MrWizard Is there any chance that question can still be migrated here ... ? It would be useful. But it's probably not possible, right? – Szabolcs Jul 15 '15 at 10:01
• Actually, the hidden functions within Splines.m are effectively the same as how the old ParametricPlot[] worked. It's just a recursive improvement of an initial sample of points up until MaxBend or PlotDivisions was satisfied. I actually refined that algorithm a bit, but it was about a decade ago, around version 5.1… – J. M.'s discontentment Jul 15 '15 at 14:03

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[DeveloperPartitionMap[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 *)

(* Visualization: *)
visualizePts[f_, {min_, max_}][pts_] := Plot[
f[x], {x, min, max},
ImageSize -> 300
]

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


• Reverse engineering lawyer alert! All WR personnel to their combat stations! – Dr. belisarius Jul 15 '15 at 16:32
• @bel, pfft, did that years ago, but then again, I never was able to do a public release. I'll probably tell my story if and when I write up an answer to this question… – J. M.'s discontentment Jul 15 '15 at 16:39
• @J.M. I dunno why, but "if and when" sounds like "never ever" :D – Dr. belisarius Jul 15 '15 at 16:50
• Dammit @bel, you're spoiling my MO! :P – J. M.'s discontentment Jul 15 '15 at 16:52

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 GraphicsSpline​. 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…

<<Version5Graphics


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:

<<Version6Graphics

• Well, the PostScript renderer was amazingly fast, and FullGraphics and AbsoluteOptions actually worked on this type of graphics output. So I miss it for those reasons at least. It's nice to see it brought back (for a brief outing, anyway) nearly a decade later. – Oleksandr R. Aug 2 '15 at 3:49
• Even when the renderer was slow in showing 3D surfaces, I nevertheless found amusement in watching it render polygons from back to front. I guess that's another thing we lost with the modern system. – J. M.'s discontentment Aug 2 '15 at 3:55
• Thanks a lot for your solution:). I discovered that norm=Sqrt[#.#]& could be replaced with Norm. – xyz Aug 12 '15 at 2:20
• @Shutao, that's true. But remember, this was early 2000's code, and old Mathematica… was a bit different. – J. M.'s discontentment Aug 12 '15 at 9:02
• @Shutao, maxb, npts, and div correspond to the settings MaxBend, PlotPoints, and PlotDivision in old Mathematica. For a finer plot, one would try one or all of increasing the settings of npts and div, and decreasing maxb (note that it is a positive angle). – J. M.'s discontentment Aug 13 '16 at 2:09

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 ExperimentalCreateNumericalFunction 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, meshimesh];
Options[mesh] = {PrecisionGoal -> 6, AccuracyGoal -> 6,
WorkingPrecision -> MachinePrecision, "InitialPoints" -> Automatic,
"ErrorNorm" -> Automatic};
mesh[f_, {var_, a_, b_}, opts : OptionsPattern[]] := meshimesh[f, {var, a, b}, opts];

Begin["mesh"];
ClearAll[meshstep, meshsplit, mesherr, mesherr0, meshimesh, meshinitialpoints];

meshinitialpoints = 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 = ExperimentalCreateNumericalFunction[{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.

• @AlexeyPopkov Maybe I'm not using the terminology correctly. Plot seems to be a hybrid. There are intervals that stop getting recursively subdivided. There are others that get subdivided, but at each step of what is controlled by MaxRecursion, each active interval of the plot is subdivided more at the max/min of sine (which I would call "local") than near it. This seems clear to me in your V7.0.1 plot, which is the plot I get in V10.2. (What seems odd about Plot is that it refines the curve well past the MaxBend or Refinement -> {ControlValue -> ...} bounds in the active areas.) – Michael E2 Aug 21 '15 at 2:57
• In other words, it's not really clear to me what goal is controlling the subdivision, and it exhibits characteristics of both. (I'm beginning to think the distinction I drew is not as clear or as helpful as I had hoped.) – Michael E2 Aug 21 '15 at 2:59

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 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 plotoption PerformanceGoal -> "Speed".

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

• As I understand, you need to avaluate g while you optimize the location of t for the sample. If you already evaluate g` many times, then you could put them all on the plot and make the plot way smoother and the whole reason for having small non-uniform set of points disappears. The whole point of optimal meshing is to precicly aproximate the expensive function without evaluating it everywhere. – Johu Aug 27 '18 at 14:08
• I'm far away from an optimal meshing algorithm. My intention is to provide a utility similar to the routine used by Plot&Co – Ulrich Neumann Aug 27 '18 at 14:35