19
$\begingroup$

Bug introduced in 7.0 and persisting through 12.0 or later

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)

enter image description here

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[].

Improve this question
$\endgroup$
6
  • $\begingroup$ What do you mean by "(computed by redefining all the NURBS functions from the beginning)"? $\endgroup$
    – Dr. belisarius
    Commented Aug 20, 2013 at 18:12
  • $\begingroup$ 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. $\endgroup$
    – user7987
    Commented Aug 20, 2013 at 18:47
  • 1
    $\begingroup$ 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. $\endgroup$
    – gdir
    Commented Mar 21, 2014 at 9:40
  • $\begingroup$ 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. $\endgroup$
    – user7987
    Commented Mar 28, 2014 at 10:47
  • $\begingroup$ @J. M. can we now directly use second/third BSplineFunction derivatives in MMA 12.2`? $\endgroup$
    – ABCDEMMM
    Commented Jun 6, 2021 at 22:06

5 Answers 5

Reset to default
11
$\begingroup$

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]

Mathematica graphics

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – user7987
    Commented Aug 21, 2013 at 13:40
  • $\begingroup$ Would you add the version header to this question if you have time? $\endgroup$
    – Mr.Wizard
    Commented Jan 24, 2016 at 3:12
9
$\begingroup$

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

correct derivatives

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

Improve this answer
$\endgroup$
8
$\begingroup$

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[
   {Red, Dashed, Arrowheads[0.03], 
    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
]

enter image description here

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]

enter image description here

Improve this answer
$\endgroup$
6
$\begingroup$

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 built-in.

plot

Improve this answer
$\endgroup$
5
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • 5
    $\begingroup$ Although I reported the issue to Wolfram support in March for version 9.01, it's still wrong in version 10.0.0.0. $\endgroup$
    – gdir
    Commented Jul 11, 2014 at 19:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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