--- a/extra/nurbs/INDEX	Thu Jun 18 15:09:49 2015 +0000
+++ b/extra/nurbs/INDEX	Fri Jun 19 16:24:09 2015 +0000
@@ -39,6 +39,7 @@
  bspdegelev
  basisfun
  basisfunder
+ basiskntins
  findspan
  numbasisfun
  tbasisfun

--- a/extra/nurbs/NEWS	Thu Jun 18 15:09:49 2015 +0000
+++ b/extra/nurbs/NEWS	Fri Jun 19 16:24:09 2015 +0000
@@ -10,6 +10,8 @@
 
  * inst/nrbbasisfun: fixed bug to return indices in the 1D case
 
+ * inst/basiskntins: added new function to compute subdivision coefficients
+
 Summary of important user-visible changes for nurbs 1.3.9:
 -------------------------------------------------------------------
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/nurbs/inst/basiskntins.m	Fri Jun 19 16:24:09 2015 +0000
@@ -0,0 +1,172 @@
+function S = basiskntins (deg,t,u)
+
+% Compute the coefficient matrix for non-uniform B-splines subdivision.
+%
+% This represents the B-spline basis given by a coarse knot vector
+%  in terms of the B-spline basis of a finer knot vector.
+%
+% The function is implemented for the univariate case. It is based on 
+%  the paper:
+%
+% G. Casciola, L. Romani, A general matrix representation for non-uniform
+%  B-spline subdivision with boundary control, ALMA-DL, University of Bologna (2007)
+%
+% Calling Sequence:
+% 
+%    S = basisfun_kntins (deg, t, u);
+%
+%    INPUT:
+%   
+%      deg - degree of the first knot vector
+%      t   - coarse knot vector
+%      u   - fine knot vector
+%   
+%    OUTPUT:
+%   
+%      B - Value of the basis functions at the points
+%          size(B)=[numel(u),(p+1)] for curves
+%          or [numel(u)*numel(v), (p+1)*(q+1)] for surfaces
+%
+%      N - Indices of the basis functions that are nonvanishing at each
+%          point. size(N) == size(B)
+%   
+%    Copyright (C) 2015 Rafael Vazquez
+%
+%    This program is free software: you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation, either version 3 of the License, or
+%    (at your option) any later version.
+
+%    This program is distributed in the hope that it will be useful,
+%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%    GNU General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+nt = length(t);
+nu = length(u);
+S = sparse (nu-deg-1,nt-deg-1);
+[t_mult,t_single,nt_s] = knot_mult(deg,t);
+[u_mult,u_single,nu_s] = knot_mult(deg,u);
+st = deg+1;
+su = deg+1;
+row = 1;
+col = 1;
+Sl = bs2bs(deg,t,u,st,su);
+S(row:deg+row,col:deg+col) = Sl;
+t_single(nt+1) = t(nt-deg);
+i = 1;
+
+for j=1:nu_s
+  if (u_single(j) == t_single(i))
+    st = st+t_mult(i);
+    col = col+t_mult(i);
+    i = i+1;
+  end
+  su = su+u_mult(j);
+  row = row+u_mult(j);
+  Sl = bs2bs(deg,t,u,st,su);
+  S(row:deg+row,col:deg+col) = Sl;
+end
+
+end
+
+
+function [t_mult,t_single,nt_s] = knot_mult(d,t)
+epsilon = 1e-14 * (t(end) - t(1));
+
+nt = length(t);
+nt_s = 0;
+m = 1;
+for i = d+2:nt-d-1
+  if ((t(i+1) - t(i)) > epsilon)
+    nt_s = nt_s+1;
+    t_mult(nt_s) = m;
+    t_single(nt_s) = t(i);
+    m=1;
+  else
+    m = m+1;
+  end
+end
+t_single(nt_s+1)=t(nt-d);
+t_mult(nt_s+1)=0;
+
+end
+
+function S = bs2bs(d,t,u,k,l)
+
+S = zeros(d+1);
+S(1,:) = bs2bs_first_row(d,t,u,k,l);
+
+for ir=1:d
+  S(ir+1,:) = bs2bs_i_row(d,t,u,k,l,ir,S(ir,:));
+end
+end
+
+function S = bs2bs_first_row(d,t,u,k,l)
+
+S = eye(1,d+1);
+for h=1:d
+  beta_2=0.0;
+  uu=u(l+1-h);
+  for j=h:-1:1
+    d1=uu-t(k+j-h);
+    d2=t(k+j)-uu;
+    beta_1=S(j)/(d2+d1);
+    S(j+1)=d1*beta_1+beta_2;
+    beta_2=d2*beta_1;
+  end
+  S(1)=beta_2;
+end
+end
+
+function Si = bs2bs_i_row(d,t,u,k,l,ir,S)
+
+Si(1) = S(1)*(t(k+1)-u(l+ir))/(t(k+1)-u(l+ir-d));
+
+for j=1:d
+  den=t(k+j+1)-u(l+ir-d);
+  fact=(t(k+j+1)-t(k+j-d))/(t(k+j)-t(k+j-d-1));
+  Si(j+1)=(fact*(S(j)*(u(l+ir)-t(k+j-d-1))-Si(j) * ...
+    (u(l+ir-d)-t(k+j-d-1)))+S(j+1)*(t(k+j+1)-u(l+ir)))/den;
+end
+end
+
+%!test
+%! knt1 = [0 0 0 1/2 1 1 1];
+%! knt2 = [0 0 0 1/4 1/2 3/4 1 1 1];
+%! C = basiskntins (2, knt1, knt2);
+%! assert (full(C), [1 0 0 0; 1/2 1/2 0 0; 0 3/4 1/4 0; 0 1/4 3/4 0; 0 0 1/2 1/2; 0 0 0 1]);
+
+%!test
+%! crv = nrbtestcrv;
+%! crv2 = nrbkntins (crv, [0.1, 0.3, 0.4, 0.5, 0.6, 0.8, 0.96 0.98]);
+%! C = basiskntins (crv.order-1,crv.knots,crv2.knots);
+%! for ii = 1:4
+%!   assert (max (abs(C*crv.coefs(ii,:)' - crv2.coefs(ii,:)')) < 1e-14 )
+%! end
+
+%!test
+%! crv = nrbtestcrv;
+%! crv2 = nrbkntins (crv, [0.50000001, 0.5000001, 0.500001, 0.50001, 0.5001]);
+%! C = basiskntins (crv.order-1,crv.knots,crv2.knots);
+%! for ii = 1:4
+%!   assert (max (abs(C*crv.coefs(ii,:)' - crv2.coefs(ii,:)')) < 1e-14 )
+%! end
+
+%!test
+%! crv = nrbtestcrv;
+%! crv2 = nrbkntins (crv, [0.1, 0.3, 0.4, 0.5, 0.6, 0.8, 0.96 0.98]);
+%! C = basiskntins (crv.order-1,crv.knots,crv2.knots);
+%! x = linspace (0, 1, 10);
+%! s = findspan (crv.number-1, crv.order-1, x, crv.knots);
+%! s2 = findspan (crv2.number-1, crv2.order-1, x, crv2.knots);
+%! N = basisfun (s, x, crv.order-1, crv.knots);
+%! N2 = basisfun (s2, x, crv2.order-1, crv2.knots);
+%! c = numbasisfun (s, x, crv.order-1, crv.knots) + 1;
+%! c2 = numbasisfun (s2, x, crv2.order-1, crv2.knots) + 1;
+%! for ii = 1:numel(x)
+%!   assert (abs(N2(ii,:) * C(c2(ii,:),c(ii,:)) - N(ii,:)) < 1e-14)
+%! end