octave_packages/nurbs-1.3.6/bspkntins.m

   1 function [ic,ik] = bspkntins(d,c,k,u)
   2
   3 % BSPKNTINS:  Insert knots into a B-Spline
   4 %
   5 % Calling Sequence:
   6 %
   7 %   [ic,ik] = bspkntins(d,c,k,u)
   8 %
   9 %  INPUT:
  10 %
  11 %    d - spline degree             integer
  12 %    c - control points            double  matrix(mc,nc)
  13 %    k - knot sequence             double  vector(nk)
  14 %    u - new knots                 double  vector(nu)
  15 %
  16 %  OUTPUT:
  17 %
  18 %    ic - new control points double  matrix(mc,nc+nu)
  19 %    ik - new knot sequence  double  vector(nk+nu)
  20 %
  21 %  Modified version of Algorithm A5.4 from 'The NURBS BOOK' pg164.
  22 %
  23 %    Copyright (C) 2000 Mark Spink, 2007 Daniel Claxton, 2010 Rafael Vazquez
  24 %
  25 %    This program is free software: you can redistribute it and/or modify
  26 %    it under the terms of the GNU General Public License as published by
  27 %    the Free Software Foundation, either version 2 of the License, or
  28 %    (at your option) any later version.
  29
  30 %    This program is distributed in the hope that it will be useful,
  31 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  32 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  33 %    GNU General Public License for more details.
  34 %
  35 %    You should have received a copy of the GNU General Public License
  36 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  37
  38 [mc,nc] = size(c);
  39 u  = sort(u);
  40 nu = numel(u);
  41 nk = numel(k);
  42                                                      %
  43                                                      % int bspkntins(int d, double *c, int mc, int nc, double *k, int nk,
  44                                                      %               double *u, int nu, double *ic, double *ik)
  45                                                      % {
  46                                                      %   int ierr = 0;
  47                                                      %   int a, b, r, l, i, j, m, n, s, q, ind;
  48                                                      %   double alfa;
  49                                                      %
  50                                                      %   double **ctrl  = vec2mat(c, mc, nc);
  51 ic = zeros(mc,nc+nu);                                %   double **ictrl = vec2mat(ic, mc, nc+nu);
  52 ik = zeros(1,nk+nu);
  53                                                      %
  54 n = size(c,2) - 1;                                   %   n = nc - 1;
  55 r = length(u) - 1;                                   %   r = nu - 1;
  56                                                      %
  57 m = n + d + 1;                                       %   m = n + d + 1;
  58 a = findspan(n, d, u(1), k);                         %   a = findspan(n, d, u[0], k);
  59 b = findspan(n, d, u(r+1), k);                       %   b = findspan(n, d, u[r], k);
  60 b = b+1;                                             %   ++b;
  61                                                      %
  62 for q=0:mc-1                                         %   for (q = 0; q < mc; q++)  {
  63    for j=0:a-d, ic(q+1,j+1) = c(q+1,j+1); end        %     for (j = 0; j <= a-d; j++) ictrl[j][q] = ctrl[j][q];
  64    for j=b-1:n, ic(q+1,j+r+2) = c(q+1,j+1); end      %     for (j = b-1; j <= n; j++) ictrl[j+r+1][q] = ctrl[j][q];
  65 end                                                  %   }
  66
  67 for j=0:a, ik(j+1) = k(j+1); end                     %   for (j = 0; j <= a; j++)   ik[j] = k[j];
  68 for j=b+d:m, ik(j+r+2) = k(j+1); end                 %   for (j = b+d; j <= m; j++) ik[j+r+1] = k[j];
  69                                                      %
  70 i = b + d - 1;                                       %   i = b + d - 1;
  71 s = b + d + r;                                       %   s = b + d + r;
  72
  73 for j=r:-1:0                                         %   for (j = r; j >= 0; j--) {
  74    while u(j+1) <= k(i+1) && i > a                   %     while (u[j] <= k[i] && i > a) {
  75        for q=0:mc-1                                  %       for (q = 0; q < mc; q++)
  76            ic(q+1,s-d) = c(q+1,i-d);                 %         ictrl[s-d-1][q] = ctrl[i-d-1][q];
  77        end
  78        ik(s+1) = k(i+1);                             %       ik[s] = k[i];
  79        s = s - 1;                                    %       --s;
  80        i = i - 1;                                    %       --i;
  81    end                                               %     }
  82
  83    for q=0:mc-1                                      %     for (q = 0; q < mc; q++)
  84        ic(q+1,s-d) = ic(q+1,s-d+1);                  %       ictrl[s-d-1][q] = ictrl[s-d][q];
  85    end
  86
  87    for l=1:d                                         %     for (l = 1; l <= d; l++)  {
  88        ind = s - d + l;                              %       ind = s - d + l;
  89        alfa = ik(s+l+1) - u(j+1);                    %       alfa = ik[s+l] - u[j];
  90        if abs(alfa) == 0                             %       if (fabs(alfa) == 0.0)
  91            for q=0:mc-1                              %         for (q = 0; q < mc; q++)
  92                ic(q+1,ind) = ic(q+1,ind+1);          %           ictrl[ind-1][q] = ictrl[ind][q];
  93            end
  94        else                                          %       else  {
  95            alfa = alfa/(ik(s+l+1) - k(i-d+l+1));     %         alfa /= (ik[s+l] - k[i-d+l]);
  96            for q=0:mc-1                              %         for (q = 0; q < mc; q++)
  97                tmp = (1-alfa)*ic(q+1,ind+1);
  98                ic(q+1,ind) = alfa*ic(q+1,ind) + tmp; %           ictrl[ind-1][q] = alfa*ictrl[ind-1][q]+(1.0-alfa)*ictrl[ind][q];
  99            end
 100        end                                           %       }
 101    end                                               %     }
 102    %
 103    ik(s+1) = u(j+1);                                 %     ik[s] = u[j];
 104    s = s - 1;                                        %     --s;
 105 end                                                  %   }
 106                                                      %
 107                                                      %   freevec2mat(ctrl);
 108                                                      %   freevec2mat(ictrl);
 109                                                      %
 110                                                      %   return ierr;
 111 end                                                  % }
 112