X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fnurbs-1.3.6%2Fbspdegelev.m;fp=octave_packages%2Fnurbs-1.3.6%2Fbspdegelev.m;h=5c68c0b9b8dacf891202edfe323c822341c4ec94;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/nurbs-1.3.6/bspdegelev.m b/octave_packages/nurbs-1.3.6/bspdegelev.m new file mode 100644 index 0000000..5c68c0b --- /dev/null +++ b/octave_packages/nurbs-1.3.6/bspdegelev.m @@ -0,0 +1,296 @@ +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 . + +[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