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

Improve this question
$\endgroup$

2 Answers 2

Reset to default
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$
3
  • $\begingroup$ Oh man, who would have thought about a bug in the underlying Mathematica functions. Thanks!! $\endgroup$ Commented Feb 3, 2020 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$ Commented Feb 3, 2020 at 20:06
  • $\begingroup$ Yes, that was an oversight from my side - thank you for pointing it out, I've corrected it now $\endgroup$ Commented Feb 3, 2020 at 21:17
4
$\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$
4
  • $\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$ Commented Feb 4, 2020 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$ Commented Feb 4, 2020 at 14:14
  • 1
    $\begingroup$ @johk95 That's what I would call perfect timing ;) $\endgroup$ Commented Feb 4, 2020 at 14:15
  • $\begingroup$ @Lukas and johk95, thanks for catching that; I've fixed it. $\endgroup$ Commented Feb 4, 2020 at 14:20

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.