function F2 = sRefine (M)
% sRefine  - refines the B-spline mesh by doubling the density of knots
%
%      F = sRefine (M)
%
%      M -   m x (n+1) matrix of B-spline knots
%      F - (2m)x(2n+1) matrix of B-spline knots defining the same surface
%
%  see also: sMBA

[m,n] = size(M);

F = sExtend(M,[1 1],[1 1]);

F2 = zeros(2*m, 2*(n-1)+1);

F1 = conv2( F, 1/64*[1 6 1; 6 36 6; 1 6 1] ); 
F2(2*(1:m)-(1), 2*(1:n)-(1)) = F1((1:m)+2, (1:n)+2);

F1 = conv2( F, 1/16*[0 1 1; 0 6 6; 0 1 1] ); 
F2(2*(1:m)-(1), 2*(1:n-1)    ) = F1((1:m)+2, (1:n-1)+3);

F1 = conv2( F, 1/16*[0 0 0; 1 6 1; 1 6 1] ); 
F2(2*(1:m)    , 2*(1:n)-(1)) = F1((1:m)+3, (1:n)+2);

F1 = conv2( F, 1/4 *[0 0 0; 0 1 1; 0 1 1] ); 
F2(2*(1:m)    , 2*(1:n-1)    ) = F1((1:m)+3, (1:n-1)+3);


return;