function [ic,ik] = bspdegelev(d,c,k,t)

% BSPDEGELEV:  Degree elevate a univariate B-Spline. 
% 
% Calling Sequence:
% 
%   [ic,ik] = bspdegelev(d,c,k,t)
% 
%   INPUT:
% 
%   d - Degree of the B-Spline.
%   c - Control points, matrix of size (dim,nc).
%   k - Knot sequence, row vector of size nk.
%   t - Raise the B-Spline degree t times.
% 
%   OUTPUT:
%
%   ic - Control points of the new B-Spline. 
%   ik - Knot vector of the new B-Spline.
% 
%    Copyright (C) 2000 Mark Spink, 2007 Daniel Claxton
%
%    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 2 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/>.

[mc,nc] = size(c);
                                                          %
                                                          % int bspdegelev(int d, double *c, int mc, int nc, double *k, int nk,
                                                          %                int t, int *nh, double *ic, double *ik)
                                                          % {
                                                          %   int row,col;
                                                          %
                                                          %   int ierr = 0;
                                                          %   int i, j, q, s, m, ph, ph2, mpi, mh, r, a, b, cind, oldr, mul;
                                                          %   int n, lbz, rbz, save, tr, kj, first, kind, last, bet, ii;
                                                          %   double inv, ua, ub, numer, den, alf, gam;
                                                          %   double **bezalfs, **bpts, **ebpts, **Nextbpts, *alfs;
                                                          %
                                                          %   double **ctrl  = vec2mat(c, mc, nc);
% ic = zeros(mc,nc*(t));                                  %   double **ictrl = vec2mat(ic, mc, nc*(t+1));
                                                          %
n = nc - 1;                                               %   n = nc - 1;
                                                          %
bezalfs =  zeros(d+1,d+t+1);                              %   bezalfs = matrix(d+1,d+t+1);
bpts = zeros(mc,d+1);                                     %   bpts = matrix(mc,d+1);
ebpts = zeros(mc,d+t+1);                                  %   ebpts = matrix(mc,d+t+1);
Nextbpts = zeros(mc,d+1);                                 %   Nextbpts = matrix(mc,d+1);
alfs = zeros(d,1);                                        %   alfs = (double *) mxMalloc(d*sizeof(double));
                                                          %
m = n + d + 1;                                            %   m = n + d + 1;
ph = d + t;                                               %   ph = d + t;
ph2 = floor(ph / 2);                                      %   ph2 = ph / 2;
                                                          %
                                                          %   // compute bezier degree elevation coefficeients
bezalfs(1,1) = 1;                                         %   bezalfs[0][0] = bezalfs[ph][d] = 1.0;
bezalfs(d+1,ph+1) = 1;                                    %

for i=1:ph2                                               %   for (i = 1; i <= ph2; i++) {
   inv = 1/bincoeff(ph,i);                                %     inv = 1.0 / bincoeff(ph,i);
   mpi = min(d,i);                                        %     mpi = min(d,i);
                                                          %
   for j=max(0,i-t):mpi                                   %     for (j = max(0,i-t); j <= mpi; j++)
       bezalfs(j+1,i+1) = inv*bincoeff(d,j)*bincoeff(t,i-j);  %       bezalfs[i][j] = inv * bincoeff(d,j) * bincoeff(t,i-j);
   end                                                       
end                                                       %   }
                                                          %
for i=ph2+1:ph-1                                          %   for (i = ph2+1; i <= ph-1; i++) {
   mpi = min(d,i);                                        %     mpi = min(d, i);
   for j=max(0,i-t):mpi                                   %     for (j = max(0,i-t); j <= mpi; j++)
       bezalfs(j+1,i+1) = bezalfs(d-j+1,ph-i+1);          %       bezalfs[i][j] = bezalfs[ph-i][d-j];
   end                                                       
end                                                       %   }
                                                          %
mh = ph;                                                  %   mh = ph;      
kind = ph+1;                                              %   kind = ph+1;
r = -1;                                                   %   r = -1;
a = d;                                                    %   a = d;
b = d+1;                                                  %   b = d+1;
cind = 1;                                                 %   cind = 1;
ua = k(1);                                                %   ua = k[0]; 
                                                          %
for ii=0:mc-1                                             %   for (ii = 0; ii < mc; ii++)
   ic(ii+1,1) = c(ii+1,1);                                %     ictrl[0][ii] = ctrl[0][ii];
end                                                       %
for i=0:ph                                                %   for (i = 0; i <= ph; i++)
   ik(i+1) = ua;                                          %     ik[i] = ua;
end                                                       %
                                                          %   // initialise first bezier seg
for i=0:d                                                 %   for (i = 0; i <= d; i++)
   for ii=0:mc-1                                          %     for (ii = 0; ii < mc; ii++)
      bpts(ii+1,i+1) = c(ii+1,i+1);                       %       bpts[i][ii] = ctrl[i][ii];
   end                                                       
end                                                       %
                                                          %   // big loop thru knot vector
while b < m                                               %   while (b < m)  {
   i = b;                                                 %     i = b;
   while b < m && k(b+1) == k(b+2)                        %     while (b < m && k[b] == k[b+1])
      b = b + 1;                                          %       b++;
   end                                                    %
   mul = b - i + 1;                                       %     mul = b - i + 1;
   mh = mh + mul + t;                                     %     mh += mul + t;
   ub = k(b+1);                                           %     ub = k[b];
   oldr = r;                                              %     oldr = r;
   r = d - mul;                                           %     r = d - mul;
                                                          %
                                                          %     // insert knot u(b) r times
   if oldr > 0                                            %     if (oldr > 0)
      lbz = floor((oldr+2)/2);                            %       lbz = (oldr+2) / 2;
   else                                                   %     else
      lbz = 1;                                            %       lbz = 1;
   end                                                    %
   
   if r > 0                                               %     if (r > 0)
      rbz = ph - floor((r+1)/2);                          %       rbz = ph - (r+1)/2;
   else                                                   %     else
      rbz = ph;                                           %       rbz = ph;
   end                                                    %
   
   if r > 0                                               %     if (r > 0) {
                                                          %       // insert knot to get bezier segment
      numer = ub - ua;                                    %       numer = ub - ua;
      for q=d:-1:mul+1                                    %       for (q = d; q > mul; q--)
         alfs(q-mul) = numer / (k(a+q+1)-ua);             %         alfs[q-mul-1] = numer / (k[a+q]-ua);
      end                                           
      
      for j=1:r                                           %       for (j = 1; j <= r; j++)  {
         save = r - j;                                    %         save = r - j;
         s = mul + j;                                     %         s = mul + j;
                                                          %
         for q=d:-1:s                                     %         for (q = d; q >= s; q--)
            for ii=0:mc-1                                 %           for (ii = 0; ii < mc; ii++)
               tmp1 = alfs(q-s+1)*bpts(ii+1,q+1); 
               tmp2 = (1-alfs(q-s+1))*bpts(ii+1,q); 
               bpts(ii+1,q+1) = tmp1 + tmp2;              %             bpts[q][ii] = alfs[q-s]*bpts[q][ii]+(1.0-alfs[q-s])*bpts[q-1][ii];
            end                                              
         end                                              %
         
         for ii=0:mc-1                                    %         for (ii = 0; ii < mc; ii++)
            Nextbpts(ii+1,save+1) = bpts(ii+1,d+1);       %           Nextbpts[save][ii] = bpts[d][ii];
         end                                                 
      end                                                 %       }
   end                                                    %     }
                                                          %     // end of insert knot
                                                          %
                                                          %     // degree elevate bezier
   for i=lbz:ph                                           %     for (i = lbz; i <= ph; i++)  {
      for ii=0:mc-1                                       %       for (ii = 0; ii < mc; ii++)
         ebpts(ii+1,i+1) = 0;                             %         ebpts[i][ii] = 0.0;
      end                                                    
      mpi = min(d, i);                                    %       mpi = min(d, i);
      for j=max(0,i-t):mpi                                %       for (j = max(0,i-t); j <= mpi; j++)
         for ii=0:mc-1                                    %         for (ii = 0; ii < mc; ii++)
            tmp1 = ebpts(ii+1,i+1); 
            tmp2 = bezalfs(j+1,i+1)*bpts(ii+1,j+1);
            ebpts(ii+1,i+1) = tmp1 + tmp2;                %           ebpts[i][ii] = ebpts[i][ii] + bezalfs[i][j]*bpts[j][ii];
         end                                                 
      end                                                    
   end                                                    %     }
                                                          %     // end of degree elevating bezier
                                                          %
   if oldr > 1                                            %     if (oldr > 1)  {
                                                          %       // must remove knot u=k[a] oldr times
      first = kind - 2;                                                    %       first = kind - 2;
      last = kind;                                        %       last = kind;
      den = ub - ua;                                      %       den = ub - ua;
      bet = floor((ub-ik(kind)) / den);                   %       bet = (ub-ik[kind-1]) / den;
                                                          %
                                                          %       // knot removal loop
      for tr=1:oldr-1                                     %       for (tr = 1; tr < oldr; tr++)  {
         i = first;                                       %         i = first;
         j = last;                                        %         j = last;
         kj = j - kind + 1;                               %         kj = j - kind + 1;
         while j-i > tr                                   %         while (j - i > tr)  {
                                                          %           // loop and compute the new control points
                                                          %           // for one removal step
            if i < cind                                   %           if (i < cind)  {
               alf = (ub-ik(i+1))/(ua-ik(i+1));           %             alf = (ub-ik[i])/(ua-ik[i]);
               for ii=0:mc-1                              %             for (ii = 0; ii < mc; ii++)
                  tmp1 = alf*ic(ii+1,i+1); 
                  tmp2 = (1-alf)*ic(ii+1,i); 
                  ic(ii+1,i+1) = tmp1 + tmp2;             %               ictrl[i][ii] = alf * ictrl[i][ii] + (1.0-alf) * ictrl[i-1][ii];
               end                                           
            end                                           %           }
            if j >= lbz                                   %           if (j >= lbz)  {
               if j-tr <= kind-ph+oldr                    %             if (j-tr <= kind-ph+oldr) {
                  gam = (ub-ik(j-tr+1)) / den;            %               gam = (ub-ik[j-tr]) / den;
                  for ii=0:mc-1                           %               for (ii = 0; ii < mc; ii++)
                     tmp1 = gam*ebpts(ii+1,kj+1); 
                     tmp2 = (1-gam)*ebpts(ii+1,kj+2); 
                     ebpts(ii+1,kj+1) = tmp1 + tmp2;      %                 ebpts[kj][ii] = gam*ebpts[kj][ii] + (1.0-gam)*ebpts[kj+1][ii];
                  end                                     %             }
               else                                       %             else  {
                  for ii=0:mc-1                           %               for (ii = 0; ii < mc; ii++)
                     tmp1 = bet*ebpts(ii+1,kj+1);                                     
                     tmp2 = (1-bet)*ebpts(ii+1,kj+2);                                     
                     ebpts(ii+1,kj+1) = tmp1 + tmp2;      %                 ebpts[kj][ii] = bet*ebpts[kj][ii] + (1.0-bet)*ebpts[kj+1][ii];
                  end                                        
               end                                        %             }
            end                                           %           }
            i = i + 1;                                    %           i++;
            j = j - 1;                                    %           j--;
            kj = kj - 1;                                  %           kj--;
         end                                              %         }
                                                          %
         first = first - 1;                               %         first--;
         last = last + 1;                                 %         last++;
      end                                                 %       }
   end                                                    %     }
                                                          %     // end of removing knot n=k[a]
                                                          %
                                                          %     // load the knot ua
   if a ~= d                                              %     if (a != d)
      for i=0:ph-oldr-1                                   %       for (i = 0; i < ph-oldr; i++)  {
         ik(kind+1) = ua;                                 %         ik[kind] = ua;
         kind = kind + 1;                                 %         kind++;
      end
   end                                                    %       }
                                                          %
                                                          %     // load ctrl pts into ic
      for j=lbz:rbz                                       %     for (j = lbz; j <= rbz; j++)  {
         for ii=0:mc-1                                    %       for (ii = 0; ii < mc; ii++)
            ic(ii+1,cind+1) = ebpts(ii+1,j+1);            %         ictrl[cind][ii] = ebpts[j][ii];
         end                                                 
         cind = cind + 1;                                 %       cind++;
      end                                                 %     }
                                                          %
      if b < m                                            %     if (b < m)  {
                                                          %       // setup for next pass thru loop
         for j=0:r-1                                      %       for (j = 0; j < r; j++)
            for ii=0:mc-1                                 %         for (ii = 0; ii < mc; ii++)
               bpts(ii+1,j+1) = Nextbpts(ii+1,j+1);       %           bpts[j][ii] = Nextbpts[j][ii];
            end                                           
         end                                              
         for j=r:d                                        %       for (j = r; j <= d; j++)
            for ii=0:mc-1                                 %         for (ii = 0; ii < mc; ii++)
               bpts(ii+1,j+1) = c(ii+1,b-d+j+1);          %           bpts[j][ii] = ctrl[b-d+j][ii];
            end                                              
         end                                                 
         a = b;                                           %       a = b;
         b = b+1;                                         %       b++;
         ua = ub;                                         %       ua = ub;
                                                          %     }
      else                                                %     else
                                                          %       // end knot
         for i=0:ph                                       %       for (i = 0; i <= ph; i++)
            ik(kind+i+1) = ub;                            %         ik[kind+i] = ub;
         end                                                 
      end                                                    
end                                                       %   }
% End big while loop                                      %   // end while loop
                                                          %
                                                          %   *nh = mh - ph - 1;
                                                          %
                                                          %   freevec2mat(ctrl);
                                                          %   freevec2mat(ictrl);
                                                          %   freematrix(bezalfs);
                                                          %   freematrix(bpts);
                                                          %   freematrix(ebpts);
                                                          %   freematrix(Nextbpts);
                                                          %   mxFree(alfs);
                                                          %
                                                          %   return(ierr);
end                                                       % }

                                                          
function b = bincoeff(n,k)
%  Computes the binomial coefficient.
%
%      ( n )      n!
%      (   ) = --------
%      ( k )   k!(n-k)!
%
%  b = bincoeff(n,k)
%
%  Algorithm from 'Numerical Recipes in C, 2nd Edition' pg215.

                                                          % double bincoeff(int n, int k)
                                                          % {
b = floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));      %   return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));
end                                                       % }

function f = factln(n)
% computes ln(n!)
if n <= 1, f = 0; return, end
f = gammaln(n+1); %log(factorial(n));
end