Using NURBS library in Matlab

Get the NURBS library from here
Extract the files to your hardisk, for example, C:\Matlab
So, your NURBS lib folder will be on C:\Matlab\nurbs-1.3.1 folder,

Create new m-file, TestNURBSlib.m, the control points data are taken from An Introduction to NURBS: With Historical Perspective book, pp. 215

Here is the control points data :

Control Points	X	Y	Z	w
B1,1	-100	-100	0	1
B1,2	-100	-50	0	1
B1,3	-100	0	0	1
B1,4	-100	50	0	1
B1,5	-100	100	0	1

B2,1	-50	-100	0	1
B2,2	-50	-50	25	1
B2,3	-50	0	50	1
B2,4	-50	50	100	1
B2,5	-50	100	0	1

B3,1	0	-100	0	1
B3,2	0	-50	25	1
B3,3	0	0	50	1
B3,4	0	50	25	1
B3,5	0	100	0	1

B4,1	50	-100	0	1
B4,2	50	-50	25	1
B4,3	50	0	150	1
B4,5	50	100	0	1
B4,4	50	50	25	1

B5,1	100	-100	0	1
B5,2	100	-50	0	1
B5,3	100	0	0	1
B5,4	100	50	0	1
B5,5	100	100	0	1

NOTE: I was modified B2,4 [-50 50 25] into B2,4 [-50 50 100]

Arrange the control points data into a matrix in the form as shown below,

pnts(:,:,1st row Control Point ) = [ B1,1-x B1,2-x B1,3-x B1,4-x B1,5-x;
    B1,1-y B1,2-y B1,3-y B1,4-y B1,5-y;
    B1,1-z B1,2-z B1,3-z B1,4-z B1,5-z;
    B1,1-w B1,2-w B1,3-w B1,4-w B1,5-w;
pnts(:,:,2nd row Control Point ) = [B2,1-x B2,2-x B2,3-x B2,4-x B2,5-x;
    B2,1-y B2,2-y B2,3-y B2,4-y B2,5-y;
    B2,1-z B2,2-z B2,3-z B2,4-z B2,5-z;
    B2,1-w B2,2-w B2,3-w B2,4-w B2,5-w;

The complete source code :

clear all;
clc;

addpath nurbs-1.3.1/nurbs/inst -BEGIN;

pnts = zeros(4,5,5);
pnts(:,:,1) = [ -100 -100 -100 -100 -100;
    -100  -50 0 50 100;
    0 0 0 0 0;
    1 1 1 1 1];
pnts(:,:,2) = [ -50 -50 -50 -50 -50; % b2,4
    -100 -50 0 50 100;
    0 25 50 100 0;
    1 1 1 1 1];

pnts(:,:,3) = [0 0 0 0 0;
    -100 -50 0 50 100;
    0 25 50 25 0;
    1 1 1 1 1]; 
    
pnts(:,:,4) = [50 50 50 50 50;
    -100 -50 0 50 100;
    0 25 150 25 0;
    1 1 1 1 1];
    
pnts(:,:,5) = [ 100 100 100 100 100;
    -100 -50 0 50 100;
    0 0 0 0 0;
    1 1 1 1 1];

 knots{1} = [0 0 0 1/3 2/3 1 1 1];
 knots{2} = [0 0 0 1/3 2/3 1 1 1];

 srf = nrbmak(pnts,knots);
 nrbplot(srf,[20 20]);hold on;
 title('NURBS surface');
 
 % create plot for the control points
 for i = 1:5
        % U direction
        plot3([pnts(1,1,i) pnts(1,2,i) pnts(1,3,i) pnts(1,4,i) pnts(1,5,i)],[pnts(2,1,i) pnts(2,2,i) pnts(2,3,i) pnts(2,4,i) pnts(2,5,i)],[pnts(3,1,i) pnts(3,2,i) pnts(3,3,i) pnts(3,4,i) pnts(3,5,i)],'b.-');
        % V direction
        plot3([pnts(1,i,1) pnts(1,i,2) pnts(1,i,3) pnts(1,i,4) pnts(1,i,5)],[pnts(2,i,1) pnts(2,i,2) pnts(2,i,3) pnts(2,i,4) pnts(2,i,5)],[pnts(3,i,1) pnts(3,i,2) pnts(3,i,3) pnts(3,i,4) pnts(3,i,5)],'r.-');
 end
 
 axis([-100 100 -100 100 0 200])
 
 hold off;

Here is the resulting NURBS surface

NURBS Surface

NOTE : For NURBS surface that has control point weight (w != 1 and w > 0), The corresponding control point must be multiplied by the weight, For example, the following control points,

pnts = zeros(4,5,5);
pnts(:,:,1) = [ -100 -100 -100 -100 -100;
    -100  -50 0 50 100;
    0 0 0 0 0;
    1 1 1 1 1];
pnts(:,:,2) = [ -50 -50 -50 -50 -50; % b2,4
    -100 -50 0 50 100;
    0 25 50 100 0;
    1 1 1 1 1];

pnts(:,:,3) = [0 0 0 0 0;
    -100 -50 0 50 100;
    0 25 50 25 0;
    1 1 1 1 1]; 

% change the B4,3 weight into 30
pnts(:,:,4) = [50 50 30*50 50 50;
    -100 -50 30*0 50 100;
    0 25 30*150 25 0;
    1 1 30 1 1];
    
pnts(:,:,5) = [ 100 100 100 100 100;
    -100 -50 0 50 100;
    0 0 0 0 0;
    1 1 1 1 1];