Codes for Blending Functions of NURBS

Recently I have been interested in computing the value of the blending functions of NURBS. I provided two codes for computing the blending functions here. One uses the recursive method and the other uses the non-recursive method.

Recursive Method
pros : Easier to implement; No need for extra loops
cons : Need extra memory for function calls

Non-recursive Method
pros : Faster than recursive method; No need for ask for extra memory for function calls
cons : Need an array and need to use for-loop

Files : main.m(for testing) , BlendingU.m(recursive) ,
BlendingFun.m(non-recursive)

In the following, I will show how to obtain Np3 (p=1,2,3, ..., 8) and get the plot of the Blending function which will be shown in the end. Note one can change the 3 to 1 and 2 to see the difference. I mean to obtain Np1 and Np2.

main.m

clear
clc

%% Knot vector & Normalization
U=[0 0 0 1 2 3 4 4 5 5 5];
U=U/max(U);

%% Non-recursive
for u=0:0.0001:1
   for p=1:1:8
      N=BlendingFun(p,3,u,U);
      subplot(211),plot(u,N);hold on
   end
end
hold off

%% Recursive
for u=0:0.0001:1
   for p=1:1:8
      N=BlendingU(p,3,u,U);
      subplot(212),plot(u,N);hold on
   end
end

BlendingU.m

%% Recursive method to compute the blending function
%% Pros : Easier to implement; No need for extra loops
%% Cons : Need extra memory for function calls

function N=BlendingU(n1,pp,uu,Ubar)

   if pp==1
      if uu>Ubar(1,n1) & uu<Ubar(1,n1+1)
         N=1;
      else
         N=0;
      end
   else
      N1N=BlendingU(n1,pp-1,uu,Ubar);
      N2N=BlendingU(n1+1,pp-1,uu,Ubar);
      N=0;
      if N1N~=0 & Ubar(1,n1+pp-1)-Ubar(1,n1)~=0
         N1=(uu-Ubar(1,n1))/(Ubar(1,n1+pp-1)-Ubar(1,n1));
         N=N1*N1N;
      end
      if N2N~=0 & Ubar(1,n1+pp)-Ubar(1,n1+1)~=0
         N2=(Ubar(n1+pp)-uu)/(Ubar(1,n1+pp)-Ubar(1,n1+1));
         N=N+N2*N2N;
      end
   end

BlendingFun.m

%% None-recursive method to compute the blending function
%% Only the associated k elements in U are used in this function
%% Need a 1 by (1+k)*k*0.5 array, say ARRAY in this function
%% Pros : Faster than recursive method; No need for ask for extra memory for
%% function calls
%% Cons : Need an array and need to use for-loop

function N=BlendingFun(i,k,u,U)

   U=U(i:i+k);
   ARRAY=zeros(1,(1+k)*k*0.5);

   for ind=0:1:k-1
      if u>=U(ind+1) & u<U(ind+1+1)
         ARRAY(ind+1)=1;
      end
   end

   layer=2;
   pt=k+1;

   while layer<=k
      for x=0:1:(k+1)-layer-1
         if U(x+1+layer-1)-U(x+1)~=0 & ARRAY(pt+x-(k+1-layer)-1)~=0
            ARRAY(pt+x)=((u-U(x+1))/(U(x+1+layer-1)-U(x+1)))*ARRAY(pt+x-(k+1-layer)-1);
         end
         if U(x+1+1+layer-1)-U(x+1+1)~=0 & ARRAY(pt+x-(k+1-layer))~=0
            ARRAY(pt+x)=ARRAY(pt+x)+((U(x+1+1+layer-1)-u)/(U(x+1+1+layer-1)-U(x+1+1)))*ARRAY(pt+x-(k+1-layer));
         end
         if x==(k+1)-layer-1
            layer=layer+1;
            pt=pt+(k+2)-layer;
         end
      end
   end

   %% Return the last element of ARRAY
   N=ARRAY((1+k)*k*0.5);