1
$\begingroup$

Description

This question comes from two questions. Namely Q1 and Q2

The defintion of B-Spline basis function as shown below:

Let $\vec{U}=\{u_0,u_1,\ldots,u_m\}$ a nondecreasing sequence of real numbers,i.e, $u_i\leq u_{i+1}\quad i=0,1,2\ldots m-1$

$$N_{i,0}(u)= \begin{cases} 1 & u_i\leq u<u_{i+1}\\ 0 & otherwise \end{cases} $$ $$N_{i,p}(u)=\frac{u-u_i}{u_{i+p}-u_i}N_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u) $$

  NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),u_] /;
   i <= Length[knots] - 2 :=
   Piecewise[
    {{1, knots[[i + 1]] <= u < knots[[i + 2]]}},0]

  coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
  coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])

  NBSpline[i_Integer, p_Integer, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &), 
    u_] /;p > 0 && i + p <= Length[knots] - 2 :=
   Module[{init, res},
    init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];
    res = First@
    Nest[
     Dot @@@
      (Thread@
        {Partition[#, 2, 1],
         With[{m = p + 2 - Length@#},
          Table[
           {coeff[u, k + 1, k + m + 1, knots],
            coeff[u, k + m + 2, k + 2, knots]}, {k, i, i + Length@# - 2}]]}) &,
    init, p]
  ]

Compare to built-in function BSplineBasis

 sortResult[x_ /; x == 0] := 0;
 sortResult[res_] := MapAt[SortBy[#, Last] &, res, 1]

 pts = {{0, 0}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
 knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};

Comparsion 1

 NBSpline[#, 0, knots, u] & /@ Range[0, 8] // Simplify

enter image description here

 PiecewiseExpand@BSplineBasis[{0, knots}, #, u] & /@ Range[0, 8]

enter image description here

Comparsion 2

sortResult /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7] // Simplify)

enter image description here

sortResult /@ (PiecewiseExpand@BSplineBasis[{1, knots}, #, u] & /@ 
Range[0, 7] // Simplify)

enter image description here


By the comparsion, I found that the result of BSplineBasis is always a closed interval, whereas the result my function NBSpline is a half closed- half open interval.

My trail:

postProcess[res_ /; res == 0] := 0;

postProcess[res_] :=
 Module[{interval, pos, expr},
  interval =
   First@Simplify@res;
  pos =
   Position[
   SortBy[interval, Last], Less];
  expr =
   ReplacePart[
    SortBy[interval, Last], pos[[-1]] -> LessEqual];
  Piecewise[expr, 0]
]

Test

 postProcess /@ (NBSpline[#, 1, knots, u] & /@ Range[0, 7])

enter image description here

Question:

Is there any good/simple method(strategy) to deal with this problem? I think my method is complex.

enter image description here

$\endgroup$
4
  • $\begingroup$ I think you have correctly implemented the formula. try specially treating the last interval If[ i<Length[knots]-2,Piecewise[..<],Piecewise[..<=]]. $\endgroup$
    – george2079
    Commented Jan 17, 2015 at 14:19
  • $\begingroup$ @george2079,This afternoon, I rewritten it as you said, however, it failed.From the result of built-in function BSplineBasis,I felt this built-in function makes the result become the closed interval. $\endgroup$
    – xyz
    Commented Jan 17, 2015 at 14:22
  • $\begingroup$ Answering "Is there any method/strategy to deal with this problem?" : NBSpline[#, 0, knots, u] & /@ Range[0, 8] /. Less :> LessEqual // Simplify $\endgroup$ Commented Jan 17, 2015 at 18:52
  • $\begingroup$ @belisarius If utilizing /. Less :> LessEqual, which will make all the sub-intervals of inequalitys become the style of $[a,b]$, $\endgroup$
    – xyz
    Commented Jan 18, 2015 at 1:51

1 Answer 1

1
$\begingroup$

a fix, specially treating the end interval:

 NBSpline[i_Integer, 0, knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &),
      u_] /; i <= Length[knots] - 2 := (
   If[ knots[[i + 2]] != knots[[-1]],
       Piecewise[{{1, knots[[i + 1]] <= u < knots[[i + 2]]}}, 0],
       Piecewise[{{1, knots[[i + 1]] <= u <= knots[[i + 2]]}}, 0] ]);
 coeff[u_, i_, j_, knots_] /; knots[[i]] == knots[[j]] := 0;
 coeff[u_, i_, j_, knots_] := (u - knots[[i]])/(knots[[j]] - knots[[i]])
 NBSpline[i_Integer, p_Integer, 
      knots_?(VectorQ[#, NumericQ] && OrderedQ[#] &), u_] /; 
           p > 0 && i + p <= Length[knots] - 2 := 
 Module[{init, res}, 
      init = Table[NBSpline[j, 0, knots, u], {j, i, i + p}];
      res = First@
      Nest[Dot @@@ (Thread@{Partition[#, 2, 1], 
          With[{m = p + 2 - Length@#}, 
              Table[{coeff[u, k + 1, k + m + 1, knots], 
              coeff[u, k + m + 2, k + 2, knots]}, {k, i, 
                  i + Length@# - 2}]]}) &, init, p]]
  pts = {{1/10, 2/10}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}};
  knots = {0, 0, 0, 1/5, 2/5, 3/5, 4/5, 1, 1, 1};
  BSplinePlot1[pts : {{_, _} ..}, knots_, opts : OptionsPattern[Plot]] :=
  Module[{p = Length@First@Split[knots] - 1, a, b},
        {a, b} = {First[knots], Last[knots]};
    f[u_] = 
        Evaluate@
           Simplify@
               Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &, 
                      pts] /. x_Piecewise :> Piecewise[x[[1]], Null];
    ParametricPlot[f[u], {u, a, b}, opts, PlotPoints -> 1000]]
    BSplinePlot1[pts, knots, ImageSize -> 600]

this bit: /. x_Piecewise :> Piecewise[x[[1]], Null] fixes the catchall behavior of Piecewise so you will get an error if you try to evaluate out of bounds. (not necessary for this example )

here is f[u] , note that we have closed just the end interval..

enter image description here

note as a practical matter @belisarius approach works as well, so long as you apply at the last stage:

 f[u_] = Evaluate@
          Simplify@
             Total@MapIndexed[NBSpline[First@#2 - 1, p, knots, u] #1 &, 
                pts] /. Less :> LessEqual /. x_Piecewise :> Piecewise[x[[1]], Null]

enter image description here

since the function is continuous it doesn't really matter if you close all the internal intervals.

$\endgroup$
1
  • $\begingroup$ Thanks for your help, In fact, I discovered a simple(but not good) method to avoid the point {0,0}. Namely, ParametricPlot[f[u], {u, a, b-10^-3}] $\endgroup$
    – xyz
    Commented Jan 20, 2015 at 6:50

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.