1 function [ic,ik] = bspdegelev(d,c,k,t)
3 % BSPDEGELEV: Degree elevate a univariate B-Spline.
7 % [ic,ik] = bspdegelev(d,c,k,t)
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.
18 % ic - Control points of the new B-Spline.
19 % ik - Knot vector of the new B-Spline.
21 % Copyright (C) 2000 Mark Spink, 2007 Daniel Claxton
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.
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.
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/>.
38 % int bspdegelev(int d, double *c, int mc, int nc, double *k, int nk,
39 % int t, int *nh, double *ic, double *ik)
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;
49 % double **ctrl = vec2mat(c, mc, nc);
50 % ic = zeros(mc,nc*(t)); % double **ictrl = vec2mat(ic, mc, nc*(t+1));
52 n = nc - 1; % n = nc - 1;
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));
60 m = n + d + 1; % m = n + d + 1;
61 ph = d + t; % ph = d + t;
62 ph2 = floor(ph / 2); % ph2 = ph / 2;
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; %
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);
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);
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];
85 kind = ph+1; % kind = ph+1;
90 ua = k(1); % ua = k[0];
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];
95 for i=0:ph % for (i = 0; i <= ph; i++)
96 ik(i+1) = ua; % ik[i] = ua;
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];
104 % // big loop thru knot vector
105 while b < m % while (b < m) {
107 while b < m && k(b+1) == k(b+2) % while (b < m && k[b] == k[b+1])
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;
116 % // insert knot u(b) r times
117 if oldr > 0 % if (oldr > 0)
118 lbz = floor((oldr+2)/2); % lbz = (oldr+2) / 2;
123 if r > 0 % if (r > 0)
124 rbz = ph - floor((r+1)/2); % rbz = ph - (r+1)/2;
126 rbz = ph; % rbz = ph;
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);
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;
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];
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];
153 % // end of insert knot
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;
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];
169 % // end of degree elevating bezier
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;
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];
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];
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];
215 first = first - 1; % first--;
216 last = last + 1; % last++;
219 % // end of removing knot n=k[a]
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++;
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];
234 cind = cind + 1; % cind++;
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];
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];
255 for i=0:ph % for (i = 0; i <= ph; i++)
256 ik(kind+i+1) = ub; % ik[kind+i] = ub;
260 % End big while loop % // end while loop
265 % freevec2mat(ictrl);
266 % freematrix(bezalfs);
269 % freematrix(Nextbpts);
276 function b = bincoeff(n,k)
277 % Computes the binomial coefficient.
285 % Algorithm from 'Numerical Recipes in C, 2nd Edition' pg215.
287 % double bincoeff(int n, int k)
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)));
292 function f = factln(n)
294 if n <= 1, f = 0; return, end
295 f = gammaln(n+1); %log(factorial(n));