Mercurial > forge

comparison extra/nurbs/inst/basiskntins.m @ 12641:e3998369a32e octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added function to compute subdivision coefficients
author rafavzqz
date Fri, 19 Jun 2015 16:24:09 +0000
parents
children 7623d14dd29c
comparison
equal deleted inserted replaced
12640:de98e4cb9248 12641:e3998369a32e
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 %
16 % S = basisfun_kntins (deg, t, u);
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