--- /dev/null
+function [ic,ik] = bspkntins(d,c,k,u)
+
+% BSPKNTINS: Insert knots into a B-Spline
+%
+% Calling Sequence:
+%
+% [ic,ik] = bspkntins(d,c,k,u)
+%
+% INPUT:
+%
+% d - spline degree integer
+% c - control points double matrix(mc,nc)
+% k - knot sequence double vector(nk)
+% u - new knots double vector(nu)
+%
+% OUTPUT:
+%
+% ic - new control points double matrix(mc,nc+nu)
+% ik - new knot sequence double vector(nk+nu)
+%
+% Modified version of Algorithm A5.4 from 'The NURBS BOOK' pg164.
+%
+% Copyright (C) 2000 Mark Spink, 2007 Daniel Claxton, 2010 Rafael Vazquez
+%
+% 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);
+u = sort(u);
+nu = numel(u);
+nk = numel(k);
+ %
+ % int bspkntins(int d, double *c, int mc, int nc, double *k, int nk,
+ % double *u, int nu, double *ic, double *ik)
+ % {
+ % int ierr = 0;
+ % int a, b, r, l, i, j, m, n, s, q, ind;
+ % double alfa;
+ %
+ % double **ctrl = vec2mat(c, mc, nc);
+ic = zeros(mc,nc+nu); % double **ictrl = vec2mat(ic, mc, nc+nu);
+ik = zeros(1,nk+nu);
+ %
+n = size(c,2) - 1; % n = nc - 1;
+r = length(u) - 1; % r = nu - 1;
+ %
+m = n + d + 1; % m = n + d + 1;
+a = findspan(n, d, u(1), k); % a = findspan(n, d, u[0], k);
+b = findspan(n, d, u(r+1), k); % b = findspan(n, d, u[r], k);
+b = b+1; % ++b;
+ %
+for q=0:mc-1 % for (q = 0; q < mc; q++) {
+ 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];
+ 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];
+end % }
+
+for j=0:a, ik(j+1) = k(j+1); end % for (j = 0; j <= a; j++) ik[j] = k[j];
+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];
+ %
+i = b + d - 1; % i = b + d - 1;
+s = b + d + r; % s = b + d + r;
+
+for j=r:-1:0 % for (j = r; j >= 0; j--) {
+ while u(j+1) <= k(i+1) && i > a % while (u[j] <= k[i] && i > a) {
+ for q=0:mc-1 % for (q = 0; q < mc; q++)
+ ic(q+1,s-d) = c(q+1,i-d); % ictrl[s-d-1][q] = ctrl[i-d-1][q];
+ end
+ ik(s+1) = k(i+1); % ik[s] = k[i];
+ s = s - 1; % --s;
+ i = i - 1; % --i;
+ end % }
+
+ for q=0:mc-1 % for (q = 0; q < mc; q++)
+ ic(q+1,s-d) = ic(q+1,s-d+1); % ictrl[s-d-1][q] = ictrl[s-d][q];
+ end
+
+ for l=1:d % for (l = 1; l <= d; l++) {
+ ind = s - d + l; % ind = s - d + l;
+ alfa = ik(s+l+1) - u(j+1); % alfa = ik[s+l] - u[j];
+ if abs(alfa) == 0 % if (fabs(alfa) == 0.0)
+ for q=0:mc-1 % for (q = 0; q < mc; q++)
+ ic(q+1,ind) = ic(q+1,ind+1); % ictrl[ind-1][q] = ictrl[ind][q];
+ end
+ else % else {
+ alfa = alfa/(ik(s+l+1) - k(i-d+l+1)); % alfa /= (ik[s+l] - k[i-d+l]);
+ for q=0:mc-1 % for (q = 0; q < mc; q++)
+ tmp = (1-alfa)*ic(q+1,ind+1);
+ 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];
+ end
+ end % }
+ end % }
+ %
+ ik(s+1) = u(j+1); % ik[s] = u[j];
+ s = s - 1; % --s;
+end % }
+ %
+ % freevec2mat(ctrl);
+ % freevec2mat(ictrl);
+ %
+ % return ierr;
+end % }
+
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff