%
% --- surfsp.m ---
%
    load  nurbs.dat;  % contains coefficients of tensor-product rendering of 
			    % surface-splines

    [coeffs, dim]  = size(nurbs);
    cpp = 36;  	% coefficients per patch
    no_patches = coeffs/cpp; 
    h = 0.5; 	% step of evaluation
    dom = 3.0;
    uKnotVector = [ -dom -dom -dom -dom 0.0 0.0 dom dom dom dom];
    vKnotVector = uKnotVector;
    u = [-3:h:3];
    v = u;
    [h1, h2] = meshgrid(u);

    % trimming -- NaN is not rendered
    w = ((abs(h1)+abs(h2)) <= dom+h);
    trim_mask = w./w; 
    mat = zeros(6,6);
    for pat = 1:no_patches,
	base = (pat-1)*cpp+1;
	Coordinate3  = nurbs(base:pat*cpp,:);
	for i=1:3,
	    mat(:) = Coordinate3(:,i)';
	    tsp = tspmak(uKnotVector,vKnotVector,mat);
	    val = fnval(spmak(uKnotVector,fnval(spmak(vKnotVector,mat),v)'),u)';
	    eval(['x_' num2str(i) '= val.*trim_mask;']);
	end;
	title('Surface Splines in Matlab');
	surf(x_1, x_2, x_3);
	hold on;
    end;
    colormap(cool);