# BSplineFunction derivatives wrong if using weights?

Bug introduced in 7.0 and persisting through 10.4

In my work, I make heavy use of non-uniform rational B-spline (NURBS) functions, defined using the function BSplineFunction[] with the option defining weights. I never before questioned the results given by Mathematica, but it seem that I discovered something that seems like a bug. Let's use a simple example : a quarter of circle. The degree, knot vector, control point vector and weights used for this are :

d = 2;
kV = {0, 0, 0, 1, 1, 1};
P = {{0, 0}, {0, 1}, {1, 1}};
W = {1, 1/Sqrt[2], 1};


I defined the two parametric functions x and y this way :

x[t_] :=
BSplineFunction[P[[All, 1]],
SplineWeights -> W, SplineDegree -> d, SplineKnots -> kV][t];
y[t_] :=
BSplineFunction[P[[All, 2]],
SplineWeights -> W, SplineDegree -> d, SplineKnots -> kV][t];


The results obtained are perfect, {x[t], y[t]} is an exact quarter of circle. The problem is when I want to have the derivatives of x and y. Here is the graph I have when I plot x'[t] (blue) and the function I should have (computed by redefining all the NURBS functions from the beginning)

We can see that Mathematica derivative is in fact x'[t] = t b2, which is in reality the derivative of the Spline function defined with the same degree, knot vector and control points, but uniform weights.(which is wrong)

I would like to know if I made a mistake somewhere, or if it is really a bug of BSplineFunction[].

-
What do you mean by "(computed by redefining all the NURBS functions from the beginning)"? – Dr. belisarius Aug 20 '13 at 18:12
I mean I redefined everything in my way: picewise polynomes for the B-Spline basis functions, rationnal polygones to have the Nurbs basis, and then the curve as the somme of basis*control points. The curve defined this way is perfectly equal to the curve defined by Mathematica, but their derivatives are not equal. – user7987 Aug 20 '13 at 18:47
Thanks a lot. It's clear now. – Dr. belisarius Aug 20 '13 at 21:56
user7987, did you get an answer to your bug report from Wolfram? I had just the same problem in Mathematica 9. Mathematica seems to compute the derivates of the BSplineFunction wrong, if weights are applied (different from 1.0). The function itself (0. derivate = compute curve point in 3D) is working as desired and computes the right coordinates even if the weigths are different from 1.0. But the derivates are only correct, if all weights are 1.0. Edit: I just tried to compute the derivates numerically with ND. This leads to correct results, while using D leads to incorrect values. – gdir Mar 21 '14 at 9:40
I did fill a bug report, and I had an anwser saying that yes, indeed it was a bug, and that it will be corrected in future releases. But not a word since. I actually finished my work with Splines using a BSplineFunction I redefined, but this was extremely unsatisfactory as it multiplied by at least 10 all my computing times. And I was completely unable to use it interactively, as was planned. – user7987 Mar 28 '14 at 10:47

Yes, there seems to be a bug in there.

You still may use BSplineFunctionif you are OK with numerical results:

<< NumericalCalculus
d = 2;
kV = {0, 0, 0, 1, 1, 1}; P = {{0, 0}, {0, 1}, {1, 1}}; W = {1, 1/Sqrt[2], 1};
x[t_] := BSplineFunction[P[[All, 1]], SplineWeights -> W, SplineDegree -> d, SplineKnots -> kV][t] /; 0 < t < 1
x[r_] := 0 /; r <= 0
x[r_] := 1 /; r >= 1
Plot[{x[t], ND[x[u], u, t, Scale -> .0001]}, {t, 0, 1}, Evaluated -> True]


-
Thank you! I filled up a bug report on Wolfram/support. I will post their answer here. Meanwhile, NumericalCalculus will do. It is way faster that defining my own Nurbs function. – user7987 Aug 21 '13 at 13:40
Would you add the version header to this question if you have time? – Mr.Wizard Jan 24 at 3:12

### The definition of B-Spline curve

$$\vec{C}(u)=\sum _{i=0}^n N_{i,p}(u) \vec{P}_i \text{ }\qquad (0\leq u\leq 1)$$

where, $\vec{P}_i$ is the control point, and the $N_ {i, p} (u)$ are the $p$-th - degree B-spline basis functions defined on the non-periodic(and non-uniform) knot vector $U$.

$$U=\{\underbrace {0,\cdots ,0}_{p+1},u_{p+1},\cdots u_{m-p-1},\underbrace {1,\cdots,1}_{p+1}\}$$

For the non-rational B-spline curve of degree $p$, its derivative is a $p-1$ degree non-rational curve. where, the new control points and knot vector are $Q_i$ and $U'$, respectively. $$\vec Q_i=p \frac{\vec P_{i+1}-\vec P_i}{u_{i+p+1}-u_{i+1}}$$

$$U'=\{\underbrace {0,\cdots ,0}_{p},u_{p+1},\cdots u_{m-p-1},\underbrace {1,\cdots,1}_{p}\}$$

I think the built-in f = BSplineFunction[2D/3D-points vector]; f' just returns a non-rational curve.

However, for the rational B-spline curve: $$\vec{C}^w(u)=\frac{\sum_{i=0}^n N_{i,p}(u)w_i\vec{P}_i}{\sum_{i=0}^n N_{i,p}(u)w_i}=\frac{\vec{A}(u)}{w(u)}$$

where, $\vec P_i=\{x_i, y_i\}$(2D curve) or $\vec P_i=\{x_i, y_i,z_i\}$(3D curve). Then

$$\left[\vec{C}^w(u)\right]'=\left[\frac{\vec A(u)}{w(u)}\right]'=\frac{\vec A'(u)w(u)-\vec A(u)w'(u)}{w^2(u)}$$

The derivative of built-in BSplineFunction[] for curve case is right when the weights is vector that contains same value. Namely, the curve is non-rational.

Here, ptsW = wi {xi, yi, zi}

BSplineDer[pts_, wgts_, {deg_, knots_}][u_?NumericQ] :=
Module[{ptsW, A, w, Au, wu, AuDer, wuDer},
ptsW = pts wgts;
A = BSplineFunction[ptsW, SplineDegree -> deg, SplineKnots -> knots];
w = BSplineFunction[wgts, SplineDegree -> deg, SplineKnots -> knots];
(*calculate the A(u) and w(u)*)
Au = A[u];
wu = w[u];
(*calculate the derivative of A(u) and w(u)*)
AuDer = A'[u];
wuDer = w'[u];
(*using the NURBS curve derivative formula*)
(wu AuDer - wuDer Au)/wu^2
]


### TEST

pts = {{1, 1}, {2, 3}, {3, -1}, {4, 1}, {5, 0}};
wgts = {1, 2, 3, 4, 5};
knots = {0, 0, 0, 0, 0.5, 1, 1, 1, 1};
f = BSplineFunction[pts, SplineDegree -> 3, SplineWeights -> {1, 2, 3, 4, 5}]

Show[
{ParametricPlot[f[t], {t, 0, 1}],
Graphics[
Table[Arrow[{f[t], f[t] + BSplineDer[pts, wgts, {3, knots}][t]/7}], {t, 0, 1, 0.1}]}],
ListPlot[f /@ Range[0, 1, 0.1],
PlotStyle -> Directive[Black, PointSize[Medium]]]},
PlotRange -> All, Axes -> False
]


Using gdir's data to test BSplineDer[]

BSplineDer[Pts, w, {5, U}][160]
f = BSplineFunction[Pts, SplineDegree -> p, SplineKnots -> U, SplineWeights -> w]
Needs["NumericalCalculus"]
ND[f[u], u, 160]


-

For anyone tired of holding their breath for the bug fix, here is how you construct a second order nurbs interpolation directly:

 b2[n_, k_, u_] :=
Piecewise@{
{(u - k[[n]])^2/((k[[n]] - k[[n + 1]]) (k[[n]] - k[[n + 2]])),
k[[n]] <= u < k[[n + 1]]},
{Total[
((u - k[[n + #]]) (u - k[[n + 2 + #]]))/
((k[[n + #]] - k[[n + 2 + #]]) (k[[n + 2]] - k[[n + 1]])) &
/@ {0, 1}],
k[[n + 1]] <= u < k[[n + 2]]},
{(u - k[[n + 3]])^2/((k[[n + 1]] - k[[n + 3]]) (k[[n + 2]] - k[[n + 3]])),
k[[n + 2]] <= u < k[[n + 3]]},
{0, True}};


The example:

 knot = {0, 0, 0, 1, 1, 1};
wt = {1, 1/Sqrt[2], 1};
p = {{0, 0}, {0, 1}, {1, 1}}
cu[u_] = Simplify[Divide @@ Total /@ Transpose[
MapIndexed[b2[First@#2 , knot, u] wt[[First@#2]] {#, 1} &, p]]];

Plot[ {cu[u][[1]], D[cu[u], u][[1]], D[x[u], u]}  , {u, 0, 1 },
PlotRange -> All , Evaluated -> True ]


where x[u] is the erroneous builtin..

-

One does not need a solution as drastic as george's for this case; after all, BSplineBasis[] is a built-in function. You can thus easily fall back on the definition of a NURBS curve:

x[t_] = (P[[All, 1]].(W Table[BSplineBasis[{d, kV}, j - 1, t], {j, Length[P]}]))/
(W.Table[BSplineBasis[{d, kV}, j - 1, t], {j, Length[P]}]);
y[t_] = (P[[All, 2]].(W Table[BSplineBasis[{d, kV}, j - 1, t], {j, Length[P]}]))/
(W.Table[BSplineBasis[{d, kV}, j - 1, t], {j, Length[P]}]);

Plot[{x'[t], y'[t]}, {t, 0, 1}, Frame -> True,
PlotStyle -> {RGBColor[7/19, 37/73, 22/31], RGBColor[59/67, 11/18, 1/7]}]


This also has the advantage of giving exact results if the input and the starting data (knots, weights, and control points) are all exact.

-

I have just sent a bug report to Wolfram myself. I'm running Mathematica 9.0.1 on Windows 8.1 64 Bit. As user7987 already found out, Mathematica seems to compute the wrong derivatives (functions D or ' or Derivative), if the BSplineFunction is rational (at least one weight different from 1.0). The ND function on the other hand gets the correct result.

Example:

U = { 151.214583, 151.214583, 151.214583, 151.214583, 151.214583,
151.214583, 465.1795421, 465.1795421, 465.1795421, 644.9800647,
644.9800647, 644.9800647, 710.5235674, 710.5235674, 710.5235674,
710.5235674, 710.5235674, 710.5235674}

Pts = { {297.7893569, -551.5786833, 816.897658},
{336.0970253, -551.5950999, 841.4645906},
{374.7963791, -551.6070814, 866.0075217},
{413.7840755, -551.6140964, 890.4320707},
{475.3962591, -551.6167492, 928.5260231},
{537.2023898, -551.6053382, 965.899484},
{559.7101647, -551.5992529, 979.3944287},
{590.4011027, -551.5882278, 997.6328514},
{621.0039434, -551.5733712, 1015.588207},
{629.1716392, -551.5691245, 1020.363635},
{637.3300529, -551.5645985, 1025.116651},
{645.4776368, -551.5597902, 1029.846112} }

w = { 1.0, 2.0, 1.7, 0.5, 1.2, 2.5, 3.0, 0.7, 1.0, 0.9, 1.3, 2.0}

p = 5


U is the knot vector, Pts the list of control points, w the list of weights and p the degree of the B-spline function.

f = BSplineFunction[Pts, SplineDegree -> p, SplineKnots -> U,
SplineWeights -> w]


If I try to compute the curve point at u = 160, I get the correct point coordinates:

f[160]
{307.10300280848, -551.58262342813, 822.86749296985}


If I try to get the first derivative at that point, Mathematica computes a wrong result:

f'[160]
{0.61075325333549, -0.00025349490633364, 0.39118702500718}


The numerical derivative function ND on the other hand gets the correct result:

Needs["NumericalCalculus"]
ND[f[u], u, 160]
{0.92580936962332, -0.00038604467360122, 0.59308687450074}


The same happens with higher derivatives at any point of the rational B-spline curve. Polynomial B-splines are processed without any problems.

-
Although I reported the issue to Wolfram support in March for version 9.01, it's still wrong in version 10.0.0.0. – gdir Jul 11 '14 at 19:13