2
$\begingroup$

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?

resulting plot

$\endgroup$
2
$\begingroup$

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)|"}
 ]

enter image description here

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

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Oh man, who would have thought about a bug in the underlying Mathematica functions. Thanks!! $\endgroup$ – johk95 Feb 3 at 19:16
  • $\begingroup$ I noticed in the specified range for the numerical integration you still use g'. That should be dg as well, I think. $\endgroup$ – johk95 Feb 3 at 20:06
  • $\begingroup$ Yes, that was an oversight from my side - thank you for pointing it out, I've corrected it now $\endgroup$ – Lukas Lang Feb 3 at 21:17
3
$\begingroup$

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

check weighted spline

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"}]]

equispaced points on a weighted B-spline

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – Lukas Lang Feb 4 at 14:14
  • $\begingroup$ 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!) $\endgroup$ – johk95 Feb 4 at 14:14
  • 1
    $\begingroup$ @johk95 That's what I would call perfect timing ;) $\endgroup$ – Lukas Lang Feb 4 at 14:15
  • $\begingroup$ @Lukas and johk95, thanks for catching that; I've fixed it. $\endgroup$ – J. M.'s technical difficulties Feb 4 at 14:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.