octave_packages/nurbs-1.3.6/bspdegelev.m

   1 function [ic,ik] = bspdegelev(d,c,k,t)
   2
   3 % BSPDEGELEV:  Degree elevate a univariate B-Spline.
   4 %
   5 % Calling Sequence:
   6 %
   7 %   [ic,ik] = bspdegelev(d,c,k,t)
   8 %
   9 %   INPUT:
  10 %
  11 %   d - Degree of the B-Spline.
  12 %   c - Control points, matrix of size (dim,nc).
  13 %   k - Knot sequence, row vector of size nk.
  14 %   t - Raise the B-Spline degree t times.
  15 %
  16 %   OUTPUT:
  17 %
  18 %   ic - Control points of the new B-Spline.
  19 %   ik - Knot vector of the new B-Spline.
  20 %
  21 %    Copyright (C) 2000 Mark Spink, 2007 Daniel Claxton
  22 %
  23 %    This program is free software: you can redistribute it and/or modify
  24 %    it under the terms of the GNU General Public License as published by
  25 %    the Free Software Foundation, either version 2 of the License, or
  26 %    (at your option) any later version.
  27
  28 %    This program is distributed in the hope that it will be useful,
  29 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  30 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  31 %    GNU General Public License for more details.
  32 %
  33 %    You should have received a copy of the GNU General Public License
  34 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  35
  36 [mc,nc] = size(c);
  37                                                           %
  38                                                           % int bspdegelev(int d, double *c, int mc, int nc, double *k, int nk,
  39                                                           %                int t, int *nh, double *ic, double *ik)
  40                                                           % {
  41                                                           %   int row,col;
  42                                                           %
  43                                                           %   int ierr = 0;
  44                                                           %   int i, j, q, s, m, ph, ph2, mpi, mh, r, a, b, cind, oldr, mul;
  45                                                           %   int n, lbz, rbz, save, tr, kj, first, kind, last, bet, ii;
  46                                                           %   double inv, ua, ub, numer, den, alf, gam;
  47                                                           %   double **bezalfs, **bpts, **ebpts, **Nextbpts, *alfs;
  48                                                           %
  49                                                           %   double **ctrl  = vec2mat(c, mc, nc);
  50 % ic = zeros(mc,nc*(t));                                  %   double **ictrl = vec2mat(ic, mc, nc*(t+1));
  51                                                           %
  52 n = nc - 1;                                               %   n = nc - 1;
  53                                                           %
  54 bezalfs =  zeros(d+1,d+t+1);                              %   bezalfs = matrix(d+1,d+t+1);
  55 bpts = zeros(mc,d+1);                                     %   bpts = matrix(mc,d+1);
  56 ebpts = zeros(mc,d+t+1);                                  %   ebpts = matrix(mc,d+t+1);
  57 Nextbpts = zeros(mc,d+1);                                 %   Nextbpts = matrix(mc,d+1);
  58 alfs = zeros(d,1);                                        %   alfs = (double *) mxMalloc(d*sizeof(double));
  59                                                           %
  60 m = n + d + 1;                                            %   m = n + d + 1;
  61 ph = d + t;                                               %   ph = d + t;
  62 ph2 = floor(ph / 2);                                      %   ph2 = ph / 2;
  63                                                           %
  64                                                           %   // compute bezier degree elevation coefficeients
  65 bezalfs(1,1) = 1;                                         %   bezalfs[0][0] = bezalfs[ph][d] = 1.0;
  66 bezalfs(d+1,ph+1) = 1;                                    %
  67
  68 for i=1:ph2                                               %   for (i = 1; i <= ph2; i++) {
  69    inv = 1/bincoeff(ph,i);                                %     inv = 1.0 / bincoeff(ph,i);
  70    mpi = min(d,i);                                        %     mpi = min(d,i);
  71                                                           %
  72    for j=max(0,i-t):mpi                                   %     for (j = max(0,i-t); j <= mpi; j++)
  73        bezalfs(j+1,i+1) = inv*bincoeff(d,j)*bincoeff(t,i-j);  %       bezalfs[i][j] = inv * bincoeff(d,j) * bincoeff(t,i-j);
  74    end
  75 end                                                       %   }
  76                                                           %
  77 for i=ph2+1:ph-1                                          %   for (i = ph2+1; i <= ph-1; i++) {
  78    mpi = min(d,i);                                        %     mpi = min(d, i);
  79    for j=max(0,i-t):mpi                                   %     for (j = max(0,i-t); j <= mpi; j++)
  80        bezalfs(j+1,i+1) = bezalfs(d-j+1,ph-i+1);          %       bezalfs[i][j] = bezalfs[ph-i][d-j];
  81    end
  82 end                                                       %   }
  83                                                           %
  84 mh = ph;                                                  %   mh = ph;
  85 kind = ph+1;                                              %   kind = ph+1;
  86 r = -1;                                                   %   r = -1;
  87 a = d;                                                    %   a = d;
  88 b = d+1;                                                  %   b = d+1;
  89 cind = 1;                                                 %   cind = 1;
  90 ua = k(1);                                                %   ua = k[0];
  91                                                           %
  92 for ii=0:mc-1                                             %   for (ii = 0; ii < mc; ii++)
  93    ic(ii+1,1) = c(ii+1,1);                                %     ictrl[0][ii] = ctrl[0][ii];
  94 end                                                       %
  95 for i=0:ph                                                %   for (i = 0; i <= ph; i++)
  96    ik(i+1) = ua;                                          %     ik[i] = ua;
  97 end                                                       %
  98                                                           %   // initialise first bezier seg
  99 for i=0:d                                                 %   for (i = 0; i <= d; i++)
 100    for ii=0:mc-1                                          %     for (ii = 0; ii < mc; ii++)
 101       bpts(ii+1,i+1) = c(ii+1,i+1);                       %       bpts[i][ii] = ctrl[i][ii];
 102    end
 103 end                                                       %
 104                                                           %   // big loop thru knot vector
 105 while b < m                                               %   while (b < m)  {
 106    i = b;                                                 %     i = b;
 107    while b < m && k(b+1) == k(b+2)                        %     while (b < m && k[b] == k[b+1])
 108       b = b + 1;                                          %       b++;
 109    end                                                    %
 110    mul = b - i + 1;                                       %     mul = b - i + 1;
 111    mh = mh + mul + t;                                     %     mh += mul + t;
 112    ub = k(b+1);                                           %     ub = k[b];
 113    oldr = r;                                              %     oldr = r;
 114    r = d - mul;                                           %     r = d - mul;
 115                                                           %
 116                                                           %     // insert knot u(b) r times
 117    if oldr > 0                                            %     if (oldr > 0)
 118       lbz = floor((oldr+2)/2);                            %       lbz = (oldr+2) / 2;
 119    else                                                   %     else
 120       lbz = 1;                                            %       lbz = 1;
 121    end                                                    %
 122
 123    if r > 0                                               %     if (r > 0)
 124       rbz = ph - floor((r+1)/2);                          %       rbz = ph - (r+1)/2;
 125    else                                                   %     else
 126       rbz = ph;                                           %       rbz = ph;
 127    end                                                    %
 128
 129    if r > 0                                               %     if (r > 0) {
 130                                                           %       // insert knot to get bezier segment
 131       numer = ub - ua;                                    %       numer = ub - ua;
 132       for q=d:-1:mul+1                                    %       for (q = d; q > mul; q--)
 133          alfs(q-mul) = numer / (k(a+q+1)-ua);             %         alfs[q-mul-1] = numer / (k[a+q]-ua);
 134       end
 135
 136       for j=1:r                                           %       for (j = 1; j <= r; j++)  {
 137          save = r - j;                                    %         save = r - j;
 138          s = mul + j;                                     %         s = mul + j;
 139                                                           %
 140          for q=d:-1:s                                     %         for (q = d; q >= s; q--)
 141             for ii=0:mc-1                                 %           for (ii = 0; ii < mc; ii++)
 142                tmp1 = alfs(q-s+1)*bpts(ii+1,q+1);
 143                tmp2 = (1-alfs(q-s+1))*bpts(ii+1,q);
 144                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];
 145             end
 146          end                                              %
 147
 148          for ii=0:mc-1                                    %         for (ii = 0; ii < mc; ii++)
 149             Nextbpts(ii+1,save+1) = bpts(ii+1,d+1);       %           Nextbpts[save][ii] = bpts[d][ii];
 150          end
 151       end                                                 %       }
 152    end                                                    %     }
 153                                                           %     // end of insert knot
 154                                                           %
 155                                                           %     // degree elevate bezier
 156    for i=lbz:ph                                           %     for (i = lbz; i <= ph; i++)  {
 157       for ii=0:mc-1                                       %       for (ii = 0; ii < mc; ii++)
 158          ebpts(ii+1,i+1) = 0;                             %         ebpts[i][ii] = 0.0;
 159       end
 160       mpi = min(d, i);                                    %       mpi = min(d, i);
 161       for j=max(0,i-t):mpi                                %       for (j = max(0,i-t); j <= mpi; j++)
 162          for ii=0:mc-1                                    %         for (ii = 0; ii < mc; ii++)
 163             tmp1 = ebpts(ii+1,i+1);
 164             tmp2 = bezalfs(j+1,i+1)*bpts(ii+1,j+1);
 165             ebpts(ii+1,i+1) = tmp1 + tmp2;                %           ebpts[i][ii] = ebpts[i][ii] + bezalfs[i][j]*bpts[j][ii];
 166          end
 167       end
 168    end                                                    %     }
 169                                                           %     // end of degree elevating bezier
 170                                                           %
 171    if oldr > 1                                            %     if (oldr > 1)  {
 172                                                           %       // must remove knot u=k[a] oldr times
 173       first = kind - 2;                                                    %       first = kind - 2;
 174       last = kind;                                        %       last = kind;
 175       den = ub - ua;                                      %       den = ub - ua;
 176       bet = floor((ub-ik(kind)) / den);                   %       bet = (ub-ik[kind-1]) / den;
 177                                                           %
 178                                                           %       // knot removal loop
 179       for tr=1:oldr-1                                     %       for (tr = 1; tr < oldr; tr++)  {
 180          i = first;                                       %         i = first;
 181          j = last;                                        %         j = last;
 182          kj = j - kind + 1;                               %         kj = j - kind + 1;
 183          while j-i > tr                                   %         while (j - i > tr)  {
 184                                                           %           // loop and compute the new control points
 185                                                           %           // for one removal step
 186             if i < cind                                   %           if (i < cind)  {
 187                alf = (ub-ik(i+1))/(ua-ik(i+1));           %             alf = (ub-ik[i])/(ua-ik[i]);
 188                for ii=0:mc-1                              %             for (ii = 0; ii < mc; ii++)
 189                   tmp1 = alf*ic(ii+1,i+1);
 190                   tmp2 = (1-alf)*ic(ii+1,i);
 191                   ic(ii+1,i+1) = tmp1 + tmp2;             %               ictrl[i][ii] = alf * ictrl[i][ii] + (1.0-alf) * ictrl[i-1][ii];
 192                end
 193             end                                           %           }
 194             if j >= lbz                                   %           if (j >= lbz)  {
 195                if j-tr <= kind-ph+oldr                    %             if (j-tr <= kind-ph+oldr) {
 196                   gam = (ub-ik(j-tr+1)) / den;            %               gam = (ub-ik[j-tr]) / den;
 197                   for ii=0:mc-1                           %               for (ii = 0; ii < mc; ii++)
 198                      tmp1 = gam*ebpts(ii+1,kj+1);
 199                      tmp2 = (1-gam)*ebpts(ii+1,kj+2);
 200                      ebpts(ii+1,kj+1) = tmp1 + tmp2;      %                 ebpts[kj][ii] = gam*ebpts[kj][ii] + (1.0-gam)*ebpts[kj+1][ii];
 201                   end                                     %             }
 202                else                                       %             else  {
 203                   for ii=0:mc-1                           %               for (ii = 0; ii < mc; ii++)
 204                      tmp1 = bet*ebpts(ii+1,kj+1);
 205                      tmp2 = (1-bet)*ebpts(ii+1,kj+2);
 206                      ebpts(ii+1,kj+1) = tmp1 + tmp2;      %                 ebpts[kj][ii] = bet*ebpts[kj][ii] + (1.0-bet)*ebpts[kj+1][ii];
 207                   end
 208                end                                        %             }
 209             end                                           %           }
 210             i = i + 1;                                    %           i++;
 211             j = j - 1;                                    %           j--;
 212             kj = kj - 1;                                  %           kj--;
 213          end                                              %         }
 214                                                           %
 215          first = first - 1;                               %         first--;
 216          last = last + 1;                                 %         last++;
 217       end                                                 %       }
 218    end                                                    %     }
 219                                                           %     // end of removing knot n=k[a]
 220                                                           %
 221                                                           %     // load the knot ua
 222    if a ~= d                                              %     if (a != d)
 223       for i=0:ph-oldr-1                                   %       for (i = 0; i < ph-oldr; i++)  {
 224          ik(kind+1) = ua;                                 %         ik[kind] = ua;
 225          kind = kind + 1;                                 %         kind++;
 226       end
 227    end                                                    %       }
 228                                                           %
 229                                                           %     // load ctrl pts into ic
 230       for j=lbz:rbz                                       %     for (j = lbz; j <= rbz; j++)  {
 231          for ii=0:mc-1                                    %       for (ii = 0; ii < mc; ii++)
 232             ic(ii+1,cind+1) = ebpts(ii+1,j+1);            %         ictrl[cind][ii] = ebpts[j][ii];
 233          end
 234          cind = cind + 1;                                 %       cind++;
 235       end                                                 %     }
 236                                                           %
 237       if b < m                                            %     if (b < m)  {
 238                                                           %       // setup for next pass thru loop
 239          for j=0:r-1                                      %       for (j = 0; j < r; j++)
 240             for ii=0:mc-1                                 %         for (ii = 0; ii < mc; ii++)
 241                bpts(ii+1,j+1) = Nextbpts(ii+1,j+1);       %           bpts[j][ii] = Nextbpts[j][ii];
 242             end
 243          end
 244          for j=r:d                                        %       for (j = r; j <= d; j++)
 245             for ii=0:mc-1                                 %         for (ii = 0; ii < mc; ii++)
 246                bpts(ii+1,j+1) = c(ii+1,b-d+j+1);          %           bpts[j][ii] = ctrl[b-d+j][ii];
 247             end
 248          end
 249          a = b;                                           %       a = b;
 250          b = b+1;                                         %       b++;
 251          ua = ub;                                         %       ua = ub;
 252                                                           %     }
 253       else                                                %     else
 254                                                           %       // end knot
 255          for i=0:ph                                       %       for (i = 0; i <= ph; i++)
 256             ik(kind+i+1) = ub;                            %         ik[kind+i] = ub;
 257          end
 258       end
 259 end                                                       %   }
 260 % End big while loop                                      %   // end while loop
 261                                                           %
 262                                                           %   *nh = mh - ph - 1;
 263                                                           %
 264                                                           %   freevec2mat(ctrl);
 265                                                           %   freevec2mat(ictrl);
 266                                                           %   freematrix(bezalfs);
 267                                                           %   freematrix(bpts);
 268                                                           %   freematrix(ebpts);
 269                                                           %   freematrix(Nextbpts);
 270                                                           %   mxFree(alfs);
 271                                                           %
 272                                                           %   return(ierr);
 273 end                                                       % }
 274
 275
 276 function b = bincoeff(n,k)
 277 %  Computes the binomial coefficient.
 278 %
 279 %      ( n )      n!
 280 %      (   ) = --------
 281 %      ( k )   k!(n-k)!
 282 %
 283 %  b = bincoeff(n,k)
 284 %
 285 %  Algorithm from 'Numerical Recipes in C, 2nd Edition' pg215.
 286
 287                                                           % double bincoeff(int n, int k)
 288                                                           % {
 289 b = floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));      %   return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));
 290 end                                                       % }
 291
 292 function f = factln(n)
 293 % computes ln(n!)
 294 if n <= 1, f = 0; return, end
 295 f = gammaln(n+1); %log(factorial(n));
 296 end