# Manual ArcLength sampling on spline differs from MeshFunction sampling

Using Michael E2's answer, I performed a manual arclength parameterisation by solving the differential equation (or: finding the inverse of a function whose derivative is known).

I compared the results with the built-in arclength mesh-sampling. Unexpectedly, they are not the same (see green and red points on plot).

Using the built-in function is no solution to this problem, because I need this method for some other sampling based on specific parameterisations (other than arclength).

Here's the code:

ptsp = {{0, 0}, {0, 2}, {3, 2}, {1, -2}, {4, -2}, {4, 0}};
g = BSplineFunction[ptsp, SplineWeights -> {1, 1, 15, 15, 1, 1}, SplineDegree -> 3];

ClearAll[s, t];
dg[t_?NumericQ] := If[t - 1. <= 0, g'[t], g'];
tfn = NDSolveValue[{t'[s] == 1/Norm[dg[t[s]]], t == 0,
WhenEvent[t[s] == 1, "StopIntegration"]},
t, {s, 0, 1 + NIntegrate[Norm[g'[t]], {t, 0, 1}]}];

ParametricPlot[g[t], {t, 0, 1},
MeshFunctions -> {"ArcLength"}, Mesh -> {20-1},
MeshStyle -> {PointSize[0.015], Green},
PlotStyle -> {Black}
Epilog -> {
PointSize[0.013], Red,
Point[g /@
tfn[Rescale[Range[0, 1, 1/20], {0, 1}, First@tfn["Domain"]]]],
PointSize[0.01], Gray, Point[g /@ Range[0, 1, 1/20]]
}
]


I also tried other methods, like this answer, which gave the same wrong result.

Does anyone have an idea why this happens for this specific B-Spline? The reason for this is that the derivative of BSplineFunction is incorrect when using weights (see this question), as can be seen from the following comparison:

Plot[
{
Norm[g[s] - g[s - 0.001]]/0.001,
Norm[g'[s]]
},
{s, 0, 9},
PlotRange -> All,
PlotLegends -> {"\!$$\*FractionBox[\(\(|$$$$g \((s)$$ - g $$(s - \ 0.001)$$\)$$|$$\), $$0.001$$]\)", "|g'(s)|"}
] Using the same crude manual derivative for the original code produces the expected result:

ptsp = {{0, 0}, {0, 2}, {3, 2}, {1, -2}, {4, -2}, {4, 0}}; g = BSplineFunction[ptsp, SplineWeights -> {1, 1, 15, 15, 1, 1}, SplineDegree -> 3];

ClearAll[s, t];
d = 0.0001;
dg[t_?NumericQ] :=
If[t < 1. - d, (g[t + d] - g[t])/d, (g - g[1 - d])/d];
tfn = NDSolveValue[{t'[s] == 1/Norm[dg[t[s]]], t == 0,
WhenEvent[t[s] == 1, "StopIntegration"]},
t, {s, 0, 1 + NIntegrate[Norm[dg[t]], {t, 0, 1}]}];

ParametricPlot[g[t], {t, 0, 1}, MeshFunctions -> {"ArcLength"},
Mesh -> {20 - 1}, MeshStyle -> {PointSize[0.015], Green},
PlotStyle -> {Black},
Epilog -> {PointSize[0.01], Red,
Point[g /@
Clip@tfn[
Rescale[Range[0, 1, 1/20], {0, 1}, First@tfn["Domain"]]]],
PointSize[0.01], Gray, Point[g /@ Range[0, 1, 1/20]]}] • Oh man, who would have thought about a bug in the underlying Mathematica functions. Thanks!! Feb 3 '20 at 19:16
• I noticed in the specified range for the numerical integration you still use g'. That should be dg as well, I think. Feb 3 '20 at 20:06
• Yes, that was an oversight from my side - thank you for pointing it out, I've corrected it now Feb 3 '20 at 21:17

Here's another alternative, based on manually building the weighted B-spline from BSplineBasis[] (similar to what was done here):

deg = 3;
pts = {{0, 0}, {0, 2}, {3, 2}, {1, -2}, {4, -2}, {4, 0}};
wts = {1, 1, 15, 15, 1, 1};

xf[t_] = (pts[[All, 1]].(wts Table[BSplineBasis[{deg, knots}, j - 1, t],
{j, Length[pts]}]))/
(wts.Table[BSplineBasis[{deg, knots}, j - 1, t], {j, Length[pts]}]);
yf[t_] = (pts[[All, 2]].(wts Table[BSplineBasis[{deg, knots}, j - 1, t],
{j, Length[pts]}]))/
(wts.Table[BSplineBasis[{deg, knots}, j - 1, t], {j, Length[pts]}]);


Check:

gg = BSplineFunction[pts, SplineWeights -> wts, SplineDegree -> deg];

ParametricPlot[{gg[t], {xf[t], yf[t]}}, {t, 0, 1},
PlotStyle -> {AbsoluteThickness, AbsoluteThickness}] Then, using the method from this answer, generate the parameter values corresponding to the equispaced points:

arc = NDSolveValue[{s'[t] == Sqrt[xf'[t]^2 + yf'[t]^2], s == 0},
s, {t, 0, 1}, Method -> "Extrapolation"];

end = arc;

With[{n = 21},
tvals = (\[FormalT] /. FindRoot[arc[\[FormalT]] == end #,
{\[FormalT], #, 0, 1}]) & /@ Subdivide[n]];


Check that the corresponding points are indeed equispaced:

Max[Abs[Differences[arc[tvals], 2]]] // Chop
0


Generate and visualize the points:

Legended[ParametricPlot[{xf[t], yf[t]}, {t, 0, 1},
Epilog -> {Directive[AbsolutePointSize, ColorData[97, 4]],
Point[Transpose[{xf /@ tvals, yf /@ tvals}]]},
MeshFunctions -> {"ArcLength"}, Mesh -> {20},
MeshStyle -> Directive[AbsolutePointSize, ColorData[97, 3]]],
PointLegend[{Directive[AbsolutePointSize, ColorData[97, 3]],
Directive[AbsolutePointSize, ColorData[97, 4]]},
{"MeshFunctions \[Rule] \"ArcLength\"",
"manually computed"}]] • I think you have an off-by-one-error in there - if you use Mesh -> {20}, the points overlap perfectly (With[{n = 20}, tvals = …] also works of course). I would be pretty surprised if MeshFunctions->"ArcLength" didn't work only for some specific functions (assuming that it performs the relevant computations on the numeric result it got from plotting the function, and doesn't use some analytic method) Feb 4 '20 at 14:14
• Well but it looks like this difference is just due to an unequal number of subdivisions. In fact if you use Mesh -> {20} everything's good. (LL you won for 10 seconds!) Feb 4 '20 at 14:14
• @johk95 That's what I would call perfect timing ;) Feb 4 '20 at 14:15
• @Lukas and johk95, thanks for catching that; I've fixed it. Feb 4 '20 at 14:20