Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
What do you mean by "(computed by redefining all the NURBS functions from the beginning)"? –  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. –  belisarius Aug 20 '13 at 21:56
1  
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 at 9:40
    
@gdir The OP hasn't been on the site since September, don't count on a quick answer –  belisarius Mar 21 at 13:52

3 Answers 3

up vote 4 down vote accepted

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

share|improve this answer
    
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

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 BSpline 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 BSpline curve. Polynomial BSplines are processed without any problems.

share|improve this answer
1  
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 at 19:13

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

enter image description here

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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