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