iconEuler Examples

B-Splines

by R. Grothmann

In this notebook, we demonstrate techniques to compute such complicated
things as B-splines.

First, we start with the Haar function, which is the B-spline
of order 0 on a set of knots, 1 on [t[k],t[k+1]].
>function N0 (x,k,t) := x>=t[k] && x<t[k+1]
We can recursively apply the Cox formula for B-Splines.

Computing B-Splines

However, we must be careful not to divide through 0. Also, we must
make sure, the function works for vector x. Since t is also a vector,
we protect it from being mapped.
>function map N (x,k,n;t) ...
if n<=0 then return N0(x,k,t); endif;
if x<=t[k] then return 0;
else if x<t[k+n+1] then
  sum=0;
  if t[k]<t[k+n] then 
    sum=sum+N(x,k,n-1,t)*(x-t[k])/(t[k+n]-t[k]); 
  endif;
  if t[k+1]<t[k+n+1] then 
    sum=sum+N(x,k+1,n-1,t)*(t[k+n+1]-x)/(t[k+n+1]-t[k+1]); 
  endif;
  return sum;
endif;
return 0;
endfunction
>plot2d("N";1,2,[1,2,3,4],a=1,b=4); insimg;

Computing B-Splines

It does also work for multiple knots.
>plot2d("N";1,2,[1,1,1,4],a=1,b=4); insimg;

Computing B-Splines

Let us write a function, plotting all B-Splines for a set of
knots.
>function plotall (t,n) ...
plot2d(none,a=min(t)+epsilon,b=max(t)-epsilon,c=0,d=1); .. set plot window
loop 1 to cols(t)-n
plot2d("N";#,n,t,add=1); .. plot #-th B-spline
end
return plot();
endfunction
>plotall([1,1,1,2,3,4,4,4],2); insimg;

Computing B-Splines

The B-splines add to 1. We compute all B-Splines (row 1:m) at
x-values (columns) and add the columns to get the sum.
>x=0:0.01:5; t=[1,1,1,2,3.1,4,4,4]; k=1:5; y=sum(Nmap(x',k;2,t))';
>plot2d(x,y,thickness=2); insimg;

Computing B-Splines

B-Spline Curves

Now we write a function, which computes a B-spline curve with
support polygon through K (2xm or 3xm). We need to use Nmap to
get all B-Splines at all points. The function will work for vectors
x.
>function BSplineCurve (x,t,K) ...
n=cols(t)-cols(K)-1;
y=zeros(rows(K),cols(x));
loop 1 to rows(K)
  y[#]=sum(K[#]*Nmap(x',1:cols(K[#]);n,t))';
end
return y;
endfunction
Another function to plot the control polygon.
>function plotPolygon (K,add=1) ...
plot2d(K[1],K[2],thickness=2,add=add);
return plot2d(K[1],K[2],points=1,style="[]",add=1);
endfunction
Let us test everything.
>t=[1,1,1,2,3,3,3]; K=[0,1,1,-1;1,1,-1,-1]; ...
>x=linspace(1.001,2.999,300); ...
>y=BSplineCurve(x,t,K); ...
>plot2d(y[1],y[2],a=-1.1,b=1.1,c=-1.1,d=1.1); ...
>plotPolygon(K); insimg;

Computing B-Splines

Nurbs

Now, we can combine everything to compute Nurbs. We need a vector of
non-zero weights. We compute the three dimensional Spline and project
to z=1.
>function NurbsCurve (x,t,K,w) ...
Kw=w*K_w;
y=BSplineCurve(x,t,Kw);
y=y/y[3];
return y[1:2];
endfunction
>w=[1,3,1/2,1]; y=NurbsCurve(x,t,K,w);
The Nurbs are closer to the corners with heigher weights.
>plot2d(y[1],y[2],add=1,color=4); insimg;

Computing B-Splines

We want to show the three dimensional weighted curve and its
projection.
>y=BSplineCurve(x,t,K*w_w,[1,3,1/2,1]); ...
>h=y/y[3]; ...
>plot3d(y[1]_h[1],y[2]_h[2],y[3]_h[3],>lines,angle=-15°,>anaglyph); insimg;

Computing B-Splines

Examples Homepage