9
$\begingroup$

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) $$

Although I know that Mathematica owns a built-in function BSplineBasis, however, I would like to write my auxiliay function $N_{i,p}(u)$ to learn the NURBS theory and Mathematica programming.

NBSpline

Alogrithm:

enter image description here

 (*=======================Caculate N[i,0](u)================================*)
 NBSpline[i_Integer, 0, u_Symbol, U : {Sequence[_] ..}] /; 
  i <= Length[U] - 2 :=
  Piecewise[
   {{1, U[[i + 1]] <= u < U[[i + 2]]},
    {0, u < U[[i + 1]] || u >= U[[i + 2]]}}]

 (*=======================Caculate N[i,p](u)================================*)

 NBSpline[i_Integer, p_Integer, u_Symbol, U : {Sequence[_] ..}?OrderedQ] /;
  p > 0 && i + p <= Length[U] - 2 :=
  Module[{ini},
   ini = Table[NBSpline[j, 0, u, U], {j, i, i + p}];
   First@Simplify@
    Nest[
     Dot @@@
      (Thread@
       {Partition[#, 2, 1],
        With[{m = i + p - Length@# + 1},
         Table[
          {(u - U[[k + 1]])/(U[[k + m + 1]] - U[[k + 1]]),
           (U[[k + m + 2]] - u)/(U[[k + m + 2]] - U[[k + 2]])}, {k, i, i + Length@# - 2}]]}) &,
   ini, p]
]

Test

NBSpline[1, 3, u, {1, 2, 3, 4, 5, 7}] // TraditionalForm

enter image description here

In my function NBSpline I avoid the condition $u_i=u_{i+1}$, because it will occured the case $\frac{0}{0}$

In the book "The NURBS book", it defines this quotient $\frac{0}{0}$ to be zero.

Question

How to deal with the condition $\frac{0}{0}$ that I sometimes need to set it to 0 ? Namely, How to deal with the condition $u_i=u_{i+1}$ in B-Spline basis function?

$\endgroup$
9
  • 3
    $\begingroup$ Quiet@Block[{Indeterminate = 0}, 0/0] $\endgroup$
    – Kuba
    Oct 15, 2014 at 12:29
  • 1
    $\begingroup$ @Kuba,Thanks for your suggestion sincerely!:-) Then I execute your code, NBSpline[1, 3, u, {1, 1, 3, 4, 5, 7}] give a result without warning, however, NBSpline[1, 3, u, {1, 1, 1, 4, 5, 7}] and NBSpline[1, 3, u, {1, 1, 2, 2, 5, 7}] give the warning informtion Power::infy: Infinite expression 1/0 encountered. >> $\endgroup$
    – xyz
    Oct 15, 2014 at 12:36
  • 4
    $\begingroup$ Perhaps one should approach the problem of understanding $N_{i,p}$ by asking how to deal with the condition $u_i = u_{i+1}$. Redefining 0/0 is not the only way, and, given your comment, it's probably not a good way. $\endgroup$
    – Michael E2
    Oct 15, 2014 at 12:47
  • $\begingroup$ @MichaelE2, I got it!Thanks! $\endgroup$
    – xyz
    Oct 15, 2014 at 12:56
  • $\begingroup$ I guess you might be interested in this book: Curves and Surfaces for CAGD - A Practical Guide by Gerald Farin. $\endgroup$
    – Silvia
    Mar 27, 2016 at 9:11

2 Answers 2

4
+500
$\begingroup$

Here is one way to deal with repeated entries in U. One can define a function to compute the coefficient, using one rule when $u_i = u_j$ and the general formula otherwise. One might put extra conditions on the patterns in coeff below, but if the function is called only within NBSpline, then one might assume the conditions are met.

ClearAll[coeff];
coeff[u_, i_, j_, U_] /; U[[i]] == U[[j]] := 0;
coeff[u_, i_, j_, U_] := (u - U[[i]])/(U[[j]] - U[[i]])

Then change the definition of NBSpline for p != 0 as follows.

NBSpline[i_Integer, p_Integer, u_Symbol, 
   U : {Sequence[_] ..}?OrderedQ] /; p > 0 && i + p <= Length[U] - 2 :=
  Module[{ini}, ini = Table[NBSpline[j, 0, u, U], {j, i, i + p}];
  First@Simplify@
    Nest[Dot @@@ (Thread@{Partition[#, 2, 1], 
          With[{m = i + p - Length@# + 1}, 
           Table[{
             coeff[u, k + 1, k + m + 1, U], 
             coeff[u, k + m + 2, k + 2, U]},
            {k, i, i + Length@# - 2}]]}) &, ini, p]]

Example:

NBSpline[1, 3, u, {1, 2, 2, 4, 5, 7}]

Mathematica graphics

The output of NBSpline[1, 3, u, {1, 2, 3, 4, 5, 7}] agrees with the output in the question.

P.S. The pattern U : {Sequence[_] ..}?OrderedQ is equivalent to U_List?OrderedQ. You might want a check that restricts U to be a list of numbers, since an ordered list of symbols such as {a, b, c} passes the OrderedQ test. The pattern U_?(VectorQ[#, NumericQ] && OrderedQ[#] &) is one way.

$\endgroup$
4
  • $\begingroup$ +1,Good idea to deal with $u_i=u_{i+1}$:-) $\endgroup$
    – xyz
    Oct 16, 2014 at 1:54
  • $\begingroup$ BTW, Is it possible to make the order of Piecewise like $2\leq u<4, 4\leq u<5, 5\leq u<7$, rather than $5\leq u<7, 2\leq u<4, 4\leq u<5$? $\endgroup$
    – xyz
    Oct 17, 2014 at 8:20
  • 1
    $\begingroup$ The easiest way would be to sort it afterwards: MapAt[SortBy[Last], piecewisefn, 1], piecewisefn is the result of NBSpline. (Use SortBy[#, Last] & instead of SortBy[Last] if you're using V9 or earlier.) $\endgroup$
    – Michael E2
    Oct 17, 2014 at 10:01
  • $\begingroup$ I can't get the same result after running your code (version 12.0). $\endgroup$ Jan 29, 2020 at 7:17
4
$\begingroup$

Is this the behaviour you need?

Solution

Unprotect[Power];
Power[0, -1] = 1
Protect[Power]

Examples

0/0
0

Mathematica graphics

Explanation

Mathematica graphics

Revert to normal

Unprotect[Power];
ClearAll[Power];
Protect[Power];
$\endgroup$
1
  • $\begingroup$ The main probelm is the condition $u_i=u_{i+1}$, your solution can indeed deal with this case.Thanks sincerely!! +1:-) I have edited my question . $\endgroup$
    – xyz
    Oct 15, 2014 at 13:00

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.