function zp = sValue( points, xp, yp )
% sValue calculates value of B-spline in a point
%	   zp = sValue( M, xp, yp )
%      M - matrix with spline node coefficients
%      xp, yp - vectors with coordinates of the points
%      zp - value of the spline surface in the points xp,yp
%

    [m,n] = size(points);

    F = sExtend(points,[1 3],[1 2]);


    zp = zeros(size(xp));
    for c = 1:length(xp),
	x = xp(c); y = yp(c);
	i = floor(x)-1; j = floor(y)-1;
	s = x-floor(x); t = y-floor(y);
	S = 1/6 * [(1-s)^3  (3*s^3-6*s^2+4) (-3*s^3+3*s^2+3*s+1)  s^3]' ;
        T = 1/6 * [(1-t)^3  (3*t^3-6*t^2+4) (-3*t^3+3*t^2+3*t+1)  t^3]' ;
        I = i+(0:3)+2;              J = j+(0:3)+2;
	zp(c) = S' * F(I,J) * T;
    end;  %for

return