Manual ArcLength sampling on spline differs from MeshFunction sampling

Question

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'[1]];
tfn = NDSolveValue[{t'[s] == 1/Norm[dg[t[s]]], t[0] == 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?

Answer 1

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};

knots = ArrayPad[Subdivide[deg], deg, "Fixed"];

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[7], AbsoluteThickness[3]}]

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] == 0}, 
                   s, {t, 0, 1}, Method -> "Extrapolation"];

end = arc[1];

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[5], ColorData[97, 4]], 
                                    Point[Transpose[{xf /@ tvals, yf /@ tvals}]]}, 
                         MeshFunctions -> {"ArcLength"}, Mesh -> {20}, 
                         MeshStyle -> Directive[AbsolutePointSize[7], ColorData[97, 3]]], 
         PointLegend[{Directive[AbsolutePointSize[7], ColorData[97, 3]], 
                      Directive[AbsolutePointSize[5], ColorData[97, 4]]},
                     {"MeshFunctions \[Rule] \"ArcLength\"", 
                      "manually computed"}]]

Answer 2

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[1] - g[1 - d])/d];
tfn = NDSolveValue[{t'[s] == 1/Norm[dg[t[s]]], t[0] == 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]]}]