|
12641
|
1 function S = basiskntins (deg,t,u) |
|
|
2 |
|
|
3 % Compute the coefficient matrix for non-uniform B-splines subdivision. |
|
|
4 % |
|
|
5 % This represents the B-spline basis given by a coarse knot vector |
|
|
6 % in terms of the B-spline basis of a finer knot vector. |
|
|
7 % |
|
|
8 % The function is implemented for the univariate case. It is based on |
|
|
9 % the paper: |
|
|
10 % |
|
|
11 % G. Casciola, L. Romani, A general matrix representation for non-uniform |
|
|
12 % B-spline subdivision with boundary control, ALMA-DL, University of Bologna (2007) |
|
|
13 % |
|
|
14 % Calling Sequence: |
|
|
15 % |
|
12668
|
16 % S = basiskntins (deg, t, u); |
|
12641
|
17 % |
|
|
18 % INPUT: |
|
|
19 % |
|
|
20 % deg - degree of the first knot vector |
|
|
21 % t - coarse knot vector |
|
|
22 % u - fine knot vector |
|
|
23 % |
|
|
24 % OUTPUT: |
|
|
25 % |
|
|
26 % B - Value of the basis functions at the points |
|
|
27 % size(B)=[numel(u),(p+1)] for curves |
|
|
28 % or [numel(u)*numel(v), (p+1)*(q+1)] for surfaces |
|
|
29 % |
|
|
30 % N - Indices of the basis functions that are nonvanishing at each |
|
|
31 % point. size(N) == size(B) |
|
|
32 % |
|
|
33 % Copyright (C) 2015 Rafael Vazquez |
|
|
34 % |
|
|
35 % This program is free software: you can redistribute it and/or modify |
|
|
36 % it under the terms of the GNU General Public License as published by |
|
|
37 % the Free Software Foundation, either version 3 of the License, or |
|
|
38 % (at your option) any later version. |
|
|
39 |
|
|
40 % This program is distributed in the hope that it will be useful, |
|
|
41 % but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
42 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
43 % GNU General Public License for more details. |
|
|
44 % |
|
|
45 % You should have received a copy of the GNU General Public License |
|
|
46 % along with this program. If not, see <http://www.gnu.org/licenses/>. |
|
|
47 |
|
|
48 nt = length(t); |
|
|
49 nu = length(u); |
|
|
50 S = sparse (nu-deg-1,nt-deg-1); |
|
|
51 [t_mult,t_single,nt_s] = knot_mult(deg,t); |
|
|
52 [u_mult,u_single,nu_s] = knot_mult(deg,u); |
|
|
53 st = deg+1; |
|
|
54 su = deg+1; |
|
|
55 row = 1; |
|
|
56 col = 1; |
|
|
57 Sl = bs2bs(deg,t,u,st,su); |
|
|
58 S(row:deg+row,col:deg+col) = Sl; |
|
|
59 t_single(nt+1) = t(nt-deg); |
|
|
60 i = 1; |
|
|
61 |
|
|
62 for j=1:nu_s |
|
|
63 if (u_single(j) == t_single(i)) |
|
|
64 st = st+t_mult(i); |
|
|
65 col = col+t_mult(i); |
|
|
66 i = i+1; |
|
|
67 end |
|
|
68 su = su+u_mult(j); |
|
|
69 row = row+u_mult(j); |
|
|
70 Sl = bs2bs(deg,t,u,st,su); |
|
|
71 S(row:deg+row,col:deg+col) = Sl; |
|
|
72 end |
|
|
73 |
|
|
74 end |
|
|
75 |
|
|
76 |
|
|
77 function [t_mult,t_single,nt_s] = knot_mult(d,t) |
|
|
78 epsilon = 1e-14 * (t(end) - t(1)); |
|
|
79 |
|
|
80 nt = length(t); |
|
|
81 nt_s = 0; |
|
|
82 m = 1; |
|
|
83 for i = d+2:nt-d-1 |
|
|
84 if ((t(i+1) - t(i)) > epsilon) |
|
|
85 nt_s = nt_s+1; |
|
|
86 t_mult(nt_s) = m; |
|
|
87 t_single(nt_s) = t(i); |
|
|
88 m=1; |
|
|
89 else |
|
|
90 m = m+1; |
|
|
91 end |
|
|
92 end |
|
|
93 t_single(nt_s+1)=t(nt-d); |
|
|
94 t_mult(nt_s+1)=0; |
|
|
95 |
|
|
96 end |
|
|
97 |
|
|
98 function S = bs2bs(d,t,u,k,l) |
|
|
99 |
|
|
100 S = zeros(d+1); |
|
|
101 S(1,:) = bs2bs_first_row(d,t,u,k,l); |
|
|
102 |
|
|
103 for ir=1:d |
|
|
104 S(ir+1,:) = bs2bs_i_row(d,t,u,k,l,ir,S(ir,:)); |
|
|
105 end |
|
|
106 end |
|
|
107 |
|
|
108 function S = bs2bs_first_row(d,t,u,k,l) |
|
|
109 |
|
|
110 S = eye(1,d+1); |
|
|
111 for h=1:d |
|
|
112 beta_2=0.0; |
|
|
113 uu=u(l+1-h); |
|
|
114 for j=h:-1:1 |
|
|
115 d1=uu-t(k+j-h); |
|
|
116 d2=t(k+j)-uu; |
|
|
117 beta_1=S(j)/(d2+d1); |
|
|
118 S(j+1)=d1*beta_1+beta_2; |
|
|
119 beta_2=d2*beta_1; |
|
|
120 end |
|
|
121 S(1)=beta_2; |
|
|
122 end |
|
|
123 end |
|
|
124 |
|
|
125 function Si = bs2bs_i_row(d,t,u,k,l,ir,S) |
|
|
126 |
|
|
127 Si(1) = S(1)*(t(k+1)-u(l+ir))/(t(k+1)-u(l+ir-d)); |
|
|
128 |
|
|
129 for j=1:d |
|
|
130 den=t(k+j+1)-u(l+ir-d); |
|
|
131 fact=(t(k+j+1)-t(k+j-d))/(t(k+j)-t(k+j-d-1)); |
|
|
132 Si(j+1)=(fact*(S(j)*(u(l+ir)-t(k+j-d-1))-Si(j) * ... |
|
|
133 (u(l+ir-d)-t(k+j-d-1)))+S(j+1)*(t(k+j+1)-u(l+ir)))/den; |
|
|
134 end |
|
|
135 end |
|
|
136 |
|
|
137 %!test |
|
|
138 %! knt1 = [0 0 0 1/2 1 1 1]; |
|
|
139 %! knt2 = [0 0 0 1/4 1/2 3/4 1 1 1]; |
|
|
140 %! C = basiskntins (2, knt1, knt2); |
|
|
141 %! assert (full(C), [1 0 0 0; 1/2 1/2 0 0; 0 3/4 1/4 0; 0 1/4 3/4 0; 0 0 1/2 1/2; 0 0 0 1]); |
|
|
142 |
|
|
143 %!test |
|
|
144 %! crv = nrbtestcrv; |
|
|
145 %! crv2 = nrbkntins (crv, [0.1, 0.3, 0.4, 0.5, 0.6, 0.8, 0.96 0.98]); |
|
|
146 %! C = basiskntins (crv.order-1,crv.knots,crv2.knots); |
|
|
147 %! for ii = 1:4 |
|
|
148 %! assert (max (abs(C*crv.coefs(ii,:)' - crv2.coefs(ii,:)')) < 1e-14 ) |
|
|
149 %! end |
|
|
150 |
|
|
151 %!test |
|
|
152 %! crv = nrbtestcrv; |
|
|
153 %! crv2 = nrbkntins (crv, [0.50000001, 0.5000001, 0.500001, 0.50001, 0.5001]); |
|
|
154 %! C = basiskntins (crv.order-1,crv.knots,crv2.knots); |
|
|
155 %! for ii = 1:4 |
|
|
156 %! assert (max (abs(C*crv.coefs(ii,:)' - crv2.coefs(ii,:)')) < 1e-14 ) |
|
|
157 %! end |
|
|
158 |
|
|
159 %!test |
|
|
160 %! crv = nrbtestcrv; |
|
|
161 %! crv2 = nrbkntins (crv, [0.1, 0.3, 0.4, 0.5, 0.6, 0.8, 0.96 0.98]); |
|
|
162 %! C = basiskntins (crv.order-1,crv.knots,crv2.knots); |
|
|
163 %! x = linspace (0, 1, 10); |
|
|
164 %! s = findspan (crv.number-1, crv.order-1, x, crv.knots); |
|
|
165 %! s2 = findspan (crv2.number-1, crv2.order-1, x, crv2.knots); |
|
|
166 %! N = basisfun (s, x, crv.order-1, crv.knots); |
|
|
167 %! N2 = basisfun (s2, x, crv2.order-1, crv2.knots); |
|
|
168 %! c = numbasisfun (s, x, crv.order-1, crv.knots) + 1; |
|
|
169 %! c2 = numbasisfun (s2, x, crv2.order-1, crv2.knots) + 1; |
|
|
170 %! for ii = 1:numel(x) |
|
|
171 %! assert (abs(N2(ii,:) * C(c2(ii,:),c(ii,:)) - N(ii,:)) < 1e-14) |
|
|
172 %! end |
