1

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

xyz
  • 605
  • 4
  • 38
  • 117
  • I think you have correctly implemented the formula. try specially treating the last interval If[ i<Length[knots]-2,Piecewise[..<],Piecewise[..<=]]. – george2079 Jan 17 '15 at 14:19
  • @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. – xyz Jan 17 '15 at 14:22
  • Answering "Is there any method/strategy to deal with this problem?" : NBSpline[#, 0, knots, u] & /@ Range[0, 8] /. Less :> LessEqual // Simplify – Dr. belisarius Jan 17 '15 at 18:52
  • @belisarius If utilizing /. Less :> LessEqual, which will make all the sub-intervals of inequalitys become the style of $[a,b]$, – xyz Jan 18 '15 at 01:51

1 Answers1

1

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.

george2079
  • 38,913
  • 1
  • 43
  • 110
  • 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}] – xyz Jan 20 '15 at 06:50