--- /dev/null
+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
\ No newline at end of file
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff