1 |
foucault |
27 |
//--------------------------------------------------------------------------- |
2 |
|
|
|
3 |
|
|
#pragma hdrstop |
4 |
|
|
|
5 |
|
|
#include "gestionversion.h" |
6 |
|
|
|
7 |
|
|
#include "stbsplines.h" |
8 |
|
|
#include <vector> |
9 |
|
|
#include "st_gestionnaire.h" |
10 |
|
|
#include "ot_systeme.h" |
11 |
|
|
#include "constantegeo.h" |
12 |
|
|
|
13 |
|
|
#include <math.h> |
14 |
|
|
|
15 |
|
|
//--------------------------------------------------------------------------- |
16 |
|
|
#pragma package(smart_init) |
17 |
|
|
|
18 |
|
|
|
19 |
|
|
ST_B_SPLINE_SURF::ST_B_SPLINE_SURF(long LigneCourante,std::string idori,int bs_degre_u,int bs_degre_v,std::vector<int> bs_indexptsctr,std::vector<int> bs_knots_multiplicities_u,std::vector<int> bs_knots_multiplicities_v,std::vector<double> bs_knots_u,std::vector<double> bs_knots_v):ST_SURFACE(LigneCourante,idori),degre_u(bs_degre_u),degre_v(bs_degre_v),sens(1) |
20 |
|
|
{ |
21 |
francois |
283 |
int r_u=bs_knots_multiplicities_u.size(); |
22 |
|
|
for (int k=0;k<r_u;k++) |
23 |
|
|
{ |
24 |
foucault |
27 |
for (int j=0;j<bs_knots_multiplicities_u[k];j++) |
25 |
francois |
283 |
knots_u.insert(knots_u.end(),bs_knots_u[k]); |
26 |
|
|
} |
27 |
|
|
int r_v=bs_knots_multiplicities_v.size(); |
28 |
|
|
for (int k=0;k<r_v;k++) |
29 |
|
|
{ |
30 |
foucault |
27 |
for (int j=0;j<bs_knots_multiplicities_v[k];j++) |
31 |
francois |
283 |
knots_v.insert(knots_v.end(),bs_knots_v[k]); |
32 |
|
|
} |
33 |
|
|
nb_point=bs_indexptsctr.size(); |
34 |
|
|
nb_ptsctr_u=knots_u.size()-degre_u-1; |
35 |
|
|
nb_ptsctr_v=knots_v.size()-degre_v-1; |
36 |
|
|
for (int i=0;i<nb_point;i++) |
37 |
|
|
{ |
38 |
foucault |
27 |
indexptsctr.insert(indexptsctr.end(),bs_indexptsctr[i]); |
39 |
francois |
283 |
} |
40 |
|
|
umin=knots_u[0]; |
41 |
|
|
umax=knots_u[knots_u.size()-1]; |
42 |
|
|
vmin=knots_v[0]; |
43 |
|
|
vmax=knots_v[knots_v.size()-1]; |
44 |
foucault |
27 |
} |
45 |
|
|
|
46 |
|
|
ST_B_SPLINE_SURF::ST_B_SPLINE_SURF(int bs_degre_u,int bs_degre_v,std::vector<double> &bs_knots_u,std::vector<double> &bs_knots_v,std::vector<double> &bs_ptsctr,std::vector<double> &bs_poids,int sense):ST_SURFACE(),degre_u(bs_degre_u),degre_v(bs_degre_v),sens(sense) |
47 |
|
|
{ |
48 |
francois |
283 |
int r_u=bs_knots_u.size(); |
49 |
|
|
for (int k=0;k<r_u;k++) |
50 |
|
|
knots_u.insert(knots_u.end(),bs_knots_u[k]); |
51 |
|
|
int r_v=bs_knots_v.size(); |
52 |
|
|
for (int k=0;k<r_v;k++) |
53 |
foucault |
27 |
knots_v.insert(knots_v.end(),bs_knots_v[k]); |
54 |
francois |
283 |
nb_point=bs_ptsctr.size()/3; |
55 |
|
|
nb_ptsctr_u=knots_u.size()-degre_u-1; |
56 |
|
|
nb_ptsctr_v=knots_v.size()-degre_v-1; |
57 |
|
|
for (int i=0;i<nb_point;i++) |
58 |
|
|
{ |
59 |
|
|
double x=bs_ptsctr[3*i]; |
60 |
|
|
double y=bs_ptsctr[3*i+1]; |
61 |
|
|
double z=bs_ptsctr[3*i+2]; |
62 |
|
|
double w=bs_poids[i]; |
63 |
|
|
ptsctr.push_back(OT_VECTEUR_4D(w*x,w*y,w*z,w)); |
64 |
|
|
} |
65 |
|
|
periodique_u=1; |
66 |
|
|
for (int i=0; i<nb_ptsctr_v; i++) |
67 |
|
|
{ |
68 |
foucault |
27 |
double *xyz1=ptsctr[i]; |
69 |
|
|
double *xyz2=ptsctr[nb_ptsctr_v*(nb_ptsctr_u-1)+i]; |
70 |
|
|
if ((!(OPERATEUR::egal (xyz1[0],xyz2[0],1E-10))) || (!(OPERATEUR::egal (xyz1[1],xyz2[1],1E-10))) || (!(OPERATEUR::egal (xyz1[2],xyz2[2],1E-10)))) |
71 |
francois |
283 |
{ |
72 |
|
|
periodique_u=0; |
73 |
|
|
break; |
74 |
foucault |
27 |
} |
75 |
francois |
283 |
} |
76 |
|
|
if (periodique_u==1) |
77 |
|
|
{ |
78 |
foucault |
27 |
int i=knots_u.size(); |
79 |
|
|
periode_u=(knots_u[i-1]-knots_u[0]); |
80 |
francois |
283 |
} |
81 |
|
|
else periode_u=0; |
82 |
|
|
periodique_v=1; |
83 |
|
|
for (int j=0; j<nb_ptsctr_u; j++) |
84 |
|
|
{ |
85 |
foucault |
27 |
double *xyz3=ptsctr[j]; |
86 |
|
|
double *xyz4=ptsctr[(j*nb_ptsctr_v+nb_ptsctr_v-1)]; |
87 |
|
|
if ((!(OPERATEUR::egal (xyz3[0],xyz4[0],1E-10))) || (!(OPERATEUR::egal (xyz3[1],xyz4[1],1E-10))) || (!(OPERATEUR::egal (xyz3[2],xyz4[2],1E-10)))) |
88 |
francois |
283 |
{ |
89 |
|
|
periodique_v=0; |
90 |
|
|
break; |
91 |
foucault |
27 |
} |
92 |
francois |
283 |
} |
93 |
|
|
if (periodique_v==1) |
94 |
|
|
{ |
95 |
foucault |
27 |
int j=knots_v.size(); |
96 |
|
|
periode_v=(knots_v[j-1]-knots_v[0]); |
97 |
francois |
283 |
} |
98 |
|
|
else periode_v=0; |
99 |
|
|
umin=knots_u[0]; |
100 |
|
|
umax=knots_u[knots_u.size()-1]; |
101 |
|
|
vmin=knots_v[0]; |
102 |
|
|
vmax=knots_v[knots_v.size()-1]; |
103 |
foucault |
27 |
} |
104 |
|
|
|
105 |
|
|
|
106 |
|
|
ST_B_SPLINE_SURF::~ST_B_SPLINE_SURF() |
107 |
|
|
{ |
108 |
|
|
} |
109 |
|
|
|
110 |
|
|
void ST_B_SPLINE_SURF::set_periodique_u(int val) |
111 |
francois |
283 |
{ |
112 |
|
|
periodique_u=val; |
113 |
|
|
if (periodique_u==1) |
114 |
|
|
{ |
115 |
foucault |
27 |
int i=knots_u.size(); |
116 |
|
|
periode_u=(knots_u[i-1]-knots_u[0]); |
117 |
francois |
283 |
} |
118 |
foucault |
27 |
} |
119 |
|
|
|
120 |
|
|
void ST_B_SPLINE_SURF::set_periodique_v(int val) |
121 |
|
|
{ |
122 |
francois |
283 |
periodique_v=val; |
123 |
|
|
if (periodique_v==1) |
124 |
|
|
{ |
125 |
foucault |
27 |
int j=knots_v.size(); |
126 |
|
|
periode_v=(knots_v[j-1]-knots_v[0]); |
127 |
francois |
283 |
} |
128 |
foucault |
27 |
} |
129 |
|
|
|
130 |
|
|
int ST_B_SPLINE_SURF::get_intervalle(int nb_pt, int degre, double t, std::vector<double> &knots) |
131 |
|
|
{ |
132 |
francois |
283 |
int compteur = 0; |
133 |
|
|
int inter; |
134 |
|
|
if (OPERATEUR::egal(t,knots[nb_pt],1E-10)==1) return nb_pt-1; |
135 |
|
|
if (OPERATEUR::egal(t,knots[degre],1E-10)==1) return degre; |
136 |
|
|
else |
137 |
|
|
{ |
138 |
foucault |
27 |
int low=degre; |
139 |
|
|
int high=nb_pt+1; |
140 |
|
|
int mid=((low+high)/2); |
141 |
|
|
while ((t<knots[mid-1]) || (t>=knots[mid])) |
142 |
francois |
283 |
{ |
143 |
|
|
if (t<knots[mid-1]) high=mid; |
144 |
|
|
else low=mid; |
145 |
|
|
mid=(low+high)/2; |
146 |
|
|
compteur++; |
147 |
|
|
} |
148 |
foucault |
27 |
inter=mid-1; |
149 |
francois |
283 |
} |
150 |
|
|
return inter; |
151 |
foucault |
27 |
} |
152 |
|
|
|
153 |
|
|
void ST_B_SPLINE_SURF::get_valeur_fonction(int inter, double t, int degre, std::vector<double> &knots,double *grand_n) |
154 |
|
|
{ |
155 |
francois |
283 |
double saved; |
156 |
|
|
grand_n[0]=1.0; |
157 |
|
|
double *gauche=new double[degre]; |
158 |
|
|
double *droite=new double[degre]; |
159 |
|
|
for (int j=1;j<=degre;j++) |
160 |
|
|
{ |
161 |
foucault |
27 |
gauche[j-1]= t-knots[inter-j+1]; |
162 |
|
|
droite[j-1]=knots[inter+j]-t; |
163 |
|
|
saved=0.0; |
164 |
|
|
for (int r=0;r<j;r++) |
165 |
francois |
283 |
{ |
166 |
|
|
double temp=grand_n[r]/(droite[r]+ gauche[j-r-1]); |
167 |
|
|
grand_n[r]=saved+droite[r]* temp; |
168 |
|
|
saved=gauche[j-r-1]*temp; |
169 |
|
|
} |
170 |
foucault |
27 |
grand_n[j]=saved; |
171 |
francois |
283 |
} |
172 |
|
|
delete [] gauche; |
173 |
|
|
delete [] droite; |
174 |
foucault |
27 |
} |
175 |
|
|
|
176 |
|
|
|
177 |
|
|
void ST_B_SPLINE_SURF::evaluer(double *uv,double *xyz) |
178 |
|
|
{ |
179 |
|
|
|
180 |
francois |
283 |
double u = uv[0]; |
181 |
|
|
double v = uv[1]; |
182 |
foucault |
27 |
|
183 |
francois |
283 |
if (sens==-1) u=umin+umax-u; |
184 |
foucault |
27 |
|
185 |
|
|
#define P(i,j) ptsctr [ (uspan - degre_u + i) * nb_ptsctr_v + (vspan - degre_v + j) ] |
186 |
|
|
|
187 |
francois |
283 |
int uspan = get_intervalle(nb_ptsctr_u,degre_u,u, knots_u); |
188 |
foucault |
27 |
|
189 |
francois |
283 |
int vspan = get_intervalle(nb_ptsctr_v,degre_v,v, knots_v); |
190 |
foucault |
27 |
|
191 |
francois |
283 |
double Nu[20]; |
192 |
foucault |
27 |
get_valeur_fonction(uspan,u,degre_u,knots_u,Nu); |
193 |
|
|
|
194 |
|
|
double Nv[20]; |
195 |
|
|
get_valeur_fonction(vspan,v,degre_v,knots_v,Nv); |
196 |
|
|
|
197 |
francois |
283 |
OT_VECTEUR_4D temp[20]; |
198 |
foucault |
27 |
|
199 |
francois |
283 |
int l; |
200 |
|
|
for (l=0;l<=degre_v;l++) { |
201 |
|
|
temp[l] =0.0 ; |
202 |
|
|
for (int k=0;k<=degre_u;k++) { |
203 |
|
|
temp[l] += Nu[k]*P(k,l) ; |
204 |
|
|
} |
205 |
foucault |
27 |
} |
206 |
francois |
283 |
OT_VECTEUR_4D sp(0,0,0,0) ; |
207 |
|
|
for (l=0;l<=degre_v;l++) { |
208 |
|
|
sp += Nv[l]*temp[l]; |
209 |
|
|
} |
210 |
foucault |
27 |
|
211 |
francois |
283 |
// transform homogeneous coordinates to 3D coordinates |
212 |
|
|
for (int i=0; i<3; i++) |
213 |
|
|
xyz[i] = sp[i]/sp.w(); |
214 |
foucault |
27 |
|
215 |
|
|
#undef P |
216 |
|
|
} |
217 |
|
|
|
218 |
|
|
void ST_B_SPLINE_SURF::deriver_fonction(int inter,double t,int degre,int dd,std::vector<double> &knots,double *f_deriver) |
219 |
|
|
{ |
220 |
|
|
#define f_deriver(i,j) (*(f_deriver+(i)*(degre+1)+j)) |
221 |
|
|
#define grand_n(i,j) (*(grand_n+(i)*(degre+1)+j)) |
222 |
|
|
#define a(i,j) (*(a+(i)*(dd+1)+j)) |
223 |
francois |
283 |
double grand_n[256]; |
224 |
|
|
double saved; |
225 |
|
|
grand_n(0,0)=1.0; |
226 |
|
|
double gauche[16]; |
227 |
|
|
double droite[16]; |
228 |
|
|
for (int j=1;j<=degre;j++) |
229 |
|
|
{ |
230 |
foucault |
27 |
gauche[j]= t-knots[inter-j+1]; |
231 |
|
|
droite[j]=knots[inter+j]-t; |
232 |
|
|
saved=0.0; |
233 |
|
|
for (int r=0;r<j;r++) |
234 |
francois |
283 |
{ |
235 |
|
|
grand_n(j,r)=(droite[r+1]+ gauche[j-r]); |
236 |
|
|
double temp = grand_n(r,j-1)/grand_n(j,r); |
237 |
foucault |
27 |
|
238 |
francois |
283 |
grand_n(r,j)= saved+droite[r+1]*temp; |
239 |
|
|
saved=gauche[j-r]*temp; |
240 |
|
|
} |
241 |
foucault |
27 |
grand_n(j,j)=saved; |
242 |
francois |
283 |
} |
243 |
|
|
for (int j=0;j<=degre;j++) |
244 |
|
|
{ |
245 |
foucault |
27 |
f_deriver(0,j)= grand_n(j,degre); |
246 |
francois |
283 |
} |
247 |
|
|
double a[256]; |
248 |
|
|
for (int r=0;r<=degre;r++) |
249 |
|
|
{ |
250 |
|
|
int s1=0; |
251 |
|
|
int s2=1; |
252 |
foucault |
27 |
a(0,0)=1.0; |
253 |
|
|
for (int k=1;k<=dd; k++) |
254 |
francois |
283 |
{ |
255 |
|
|
double d=0.0; |
256 |
|
|
int rk=r-k; |
257 |
|
|
int pk=degre-k; |
258 |
|
|
if (r>=k) |
259 |
|
|
{ |
260 |
|
|
a(s2,0)=a(s1,0)/grand_n(pk+1,rk); |
261 |
|
|
d= a(s2,0)* grand_n(rk,pk); |
262 |
|
|
} |
263 |
|
|
int j1; |
264 |
|
|
int j2; |
265 |
|
|
if (rk>=-1) j1=1; |
266 |
|
|
else j1= -rk; |
267 |
|
|
if (r-1<=pk) j2=k-1; |
268 |
|
|
else j2=degre-r; |
269 |
|
|
for (int j=j1;j<=j2;j++) |
270 |
|
|
{ |
271 |
|
|
a(s2,j) = (a(s1,j)-a(s1,j-1))/grand_n(pk+1,rk+j); |
272 |
|
|
d+=a(s2,j)*grand_n(rk+j,pk); |
273 |
|
|
} |
274 |
|
|
if (r<=pk) |
275 |
|
|
{ |
276 |
|
|
a(s2,k) = -a(s1,k-1)/grand_n(pk+1,r); |
277 |
|
|
d+=a(s2,k)*grand_n(r,pk); |
278 |
|
|
} |
279 |
|
|
f_deriver(k,r)=d; |
280 |
|
|
int j=s1; |
281 |
|
|
s1=s2; |
282 |
|
|
s2=j; |
283 |
foucault |
27 |
} |
284 |
francois |
283 |
} |
285 |
|
|
int r=degre; |
286 |
|
|
for (int k=1;k<=dd;k++) |
287 |
|
|
{ |
288 |
|
|
for (int j=0;j<=degre;j++) |
289 |
foucault |
27 |
{ |
290 |
francois |
283 |
f_deriver(k,j)*=r; |
291 |
|
|
} |
292 |
foucault |
27 |
r*=(degre-k); |
293 |
francois |
283 |
} |
294 |
foucault |
27 |
#undef f_deriver |
295 |
|
|
#undef grand_n |
296 |
francois |
283 |
#undef a |
297 |
foucault |
27 |
} |
298 |
|
|
|
299 |
|
|
// Setup the binomial coefficients into th matrix Bin |
300 |
|
|
// Bin(i,j) = (i j) |
301 |
|
|
// The binomical coefficients are defined as follow |
302 |
|
|
// (n) n! |
303 |
|
|
// (k) = k!(n-k)! 0<=k<=n |
304 |
francois |
283 |
// and the following relationship applies |
305 |
foucault |
27 |
// (n+1) (n) ( n ) |
306 |
|
|
// ( k ) = (k) + (k-1) |
307 |
francois |
283 |
void ST_B_SPLINE_SURF::binomialCoef(double * Bin, int d) { |
308 |
|
|
int n,k; |
309 |
|
|
// Setup the first line |
310 |
|
|
Bin[0] = 1.0 ; |
311 |
|
|
for (k=d-1;k>0;--k) |
312 |
|
|
Bin[k] = 0.0 ; |
313 |
|
|
// Setup the other lines |
314 |
|
|
for (n=0;n<d-1;n++) { |
315 |
|
|
Bin[(n+1)*d+0] = 1.0 ; |
316 |
|
|
for (k=1;k<d;k++) |
317 |
|
|
if (n+1<k) |
318 |
|
|
Bin[(n+1)*d+k] = 0.0 ; |
319 |
|
|
else |
320 |
|
|
Bin[(n+1)*d+k] = Bin[n*d+k] + Bin[n*d+k-1] ; |
321 |
|
|
} |
322 |
foucault |
27 |
} |
323 |
|
|
|
324 |
|
|
void ST_B_SPLINE_SURF::deriver_bs_kieme(int nb_ptsctr_u,int degre_u,std::vector<double> &knots_u,int nb_ptsctr_v,int degre_v,std::vector<double> &knots_v,double u,double v,int d,OT_VECTEUR_4D * skl) |
325 |
|
|
{ |
326 |
|
|
#define skl(i,j) skl[(i)*(d+1)+j] |
327 |
|
|
#define Nu(i,j) Nu[(i)*(degre_u+1)+j] |
328 |
|
|
#define Nv(i,j) Nv[(i)*(degre_v+1)+j] |
329 |
|
|
#define P(i,j) ptsctr[(i)*(nb_ptsctr_v)+j] |
330 |
|
|
|
331 |
francois |
283 |
int i,j,k,l,s,r,dd; |
332 |
foucault |
27 |
|
333 |
francois |
283 |
int du = std::min(d,degre_u); |
334 |
|
|
for ( k=degre_u+1;k<=d;k++) |
335 |
|
|
{ |
336 |
|
|
for (int l=0;l<=d-k;k++) |
337 |
foucault |
27 |
{ |
338 |
francois |
283 |
skl(k,l)=0.0; |
339 |
foucault |
27 |
} |
340 |
francois |
283 |
} |
341 |
|
|
int dv = std::min(d,degre_v); |
342 |
|
|
for ( l=degre_v+1;l<=d;l++) |
343 |
|
|
{ |
344 |
|
|
for ( k=0;k<=d-l;k++) |
345 |
foucault |
27 |
{ |
346 |
francois |
283 |
skl(k,l)=0.0; |
347 |
foucault |
27 |
} |
348 |
francois |
283 |
} |
349 |
foucault |
27 |
|
350 |
francois |
283 |
int inter_u = get_intervalle(nb_ptsctr_u,degre_u,u, knots_u); |
351 |
|
|
double Nu [ 256 ]; |
352 |
|
|
deriver_fonction(inter_u,u,degre_u,du,knots_u,Nu); |
353 |
foucault |
27 |
|
354 |
francois |
283 |
int inter_v = get_intervalle(nb_ptsctr_v,degre_v,v, knots_v); |
355 |
|
|
double Nv [ 256 ]; |
356 |
|
|
deriver_fonction(inter_v,v,degre_v,dv,knots_v,Nv); |
357 |
|
|
|
358 |
|
|
OT_VECTEUR_4D temp [ 16 ]; |
359 |
|
|
for ( k=0; k<=du;k++) |
360 |
|
|
{ |
361 |
|
|
for ( s=0; s<=degre_v;s++) |
362 |
foucault |
27 |
{ |
363 |
francois |
283 |
temp[s]=0.0; |
364 |
|
|
for ( r=0; r<=degre_u;r++) |
365 |
|
|
temp[s]+=P((inter_u-degre_u+r),(inter_v-degre_v+s))*Nu(k,r); |
366 |
foucault |
27 |
} |
367 |
|
|
dd = std::min(d-k,dv); |
368 |
francois |
283 |
for ( l=0; l<=dd;l++) |
369 |
foucault |
27 |
{ |
370 |
francois |
283 |
skl(k,l) = 0.0; |
371 |
|
|
for ( s=0; s<=degre_v;s++) |
372 |
|
|
skl(k,l) += temp[s]*Nv(l,s); |
373 |
foucault |
27 |
} |
374 |
francois |
283 |
} |
375 |
foucault |
27 |
|
376 |
|
|
#undef skl |
377 |
|
|
#undef f_deriver_u |
378 |
|
|
#undef f_deriver_v |
379 |
|
|
#undef P |
380 |
|
|
} |
381 |
|
|
|
382 |
|
|
void ST_B_SPLINE_SURF::deriver_kieme(double *uv,int d,OT_VECTEUR_4D *skl) |
383 |
|
|
{ |
384 |
|
|
#define skl(i,j) (*(skl+(i)*(d+1)+j)) |
385 |
|
|
#define ders(i,j) (*(ders+(i)*(d+1)+j)) |
386 |
francois |
283 |
double u = uv[0]; |
387 |
|
|
double v = uv[1]; |
388 |
|
|
int i,j,k,l; |
389 |
|
|
OT_VECTEUR_4D pv,pv2; |
390 |
foucault |
27 |
|
391 |
francois |
283 |
OT_VECTEUR_4D ders [256]; |
392 |
foucault |
27 |
|
393 |
francois |
283 |
double Bin[256]; |
394 |
|
|
int dbin=d+1; |
395 |
|
|
binomialCoef(Bin, dbin); |
396 |
foucault |
27 |
#define Bin(i,j) (*(Bin+(i)*(dbin)+j)) |
397 |
|
|
|
398 |
francois |
283 |
deriver_bs_kieme(nb_ptsctr_u,degre_u,knots_u,nb_ptsctr_v,degre_v,knots_v,u,v,d,ders); |
399 |
foucault |
27 |
|
400 |
francois |
283 |
for (k=0;k<=d;++k) |
401 |
foucault |
27 |
{ |
402 |
francois |
283 |
for (l=0;l<=d-k;++l) |
403 |
|
|
{ |
404 |
|
|
pv = ders(k,l); |
405 |
|
|
for (j=1;j<=l;j++) |
406 |
|
|
pv -= Bin(l,j)*ders(0,j).w()*skl(k,l-j) ; |
407 |
|
|
for (i=1;i<=k;i++) |
408 |
|
|
{ |
409 |
|
|
pv -= Bin(k,i)*ders(i,0).w()*skl(k-i,l) ; |
410 |
|
|
pv2 = 0.0 ; |
411 |
|
|
for (j=1;j<=l;j++) |
412 |
|
|
pv2 += Bin(l,j)*ders(i,j).w()*skl(k-i,l-j) ; |
413 |
|
|
pv -= Bin(k,i)*pv2 ; |
414 |
|
|
} |
415 |
|
|
skl(k,l) = pv; |
416 |
|
|
skl(k,l) /= ders[0].w(); |
417 |
|
|
} |
418 |
foucault |
27 |
} |
419 |
|
|
#undef skl |
420 |
|
|
#undef ders |
421 |
francois |
283 |
#undef Bin |
422 |
foucault |
27 |
} |
423 |
|
|
|
424 |
|
|
|
425 |
|
|
void ST_B_SPLINE_SURF::deriver(double *uv,double *xyzdu, double *xyzdv) |
426 |
|
|
{ |
427 |
francois |
283 |
OT_VECTEUR_4D skl[4]; |
428 |
|
|
double uvl[3]; |
429 |
|
|
uvl[0]=uv[0]; |
430 |
|
|
uvl[1]=uv[1]; |
431 |
|
|
if (sens==-1) uvl[0]=umin+umax-uvl[0]; |
432 |
|
|
deriver_kieme(uvl,1,skl); |
433 |
foucault |
27 |
#define skl(i,j) skl[(i)*(2)+j] |
434 |
|
|
|
435 |
francois |
283 |
for (int i=0; i<3; i++) |
436 |
|
|
xyzdu[i] = sens * skl(1,0)[i]; |
437 |
|
|
|
438 |
|
|
for (int i=0; i<3; i++) |
439 |
|
|
xyzdv[i]=skl(0,1)[i]; |
440 |
|
|
|
441 |
|
|
#undef skl |
442 |
foucault |
27 |
} |
443 |
|
|
|
444 |
|
|
void ST_B_SPLINE_SURF::deriver_seconde(double *uv,double* xyzduu,double* xyzduv,double* xyzdvv,double *xyz, double *xyzdu, double *xyzdv) |
445 |
francois |
283 |
{ |
446 |
foucault |
27 |
#define skl(i,j) skl[(i)*(3)+j] |
447 |
francois |
283 |
OT_VECTEUR_4D skl[9]; |
448 |
foucault |
27 |
|
449 |
francois |
283 |
double uvl[3]; |
450 |
|
|
uvl[0]=uv[0]; |
451 |
|
|
uvl[1]=uv[1]; |
452 |
|
|
if (sens==-1) uvl[0]=umin+umax-uvl[0]; |
453 |
|
|
deriver_kieme(uvl,2,skl); |
454 |
foucault |
27 |
|
455 |
francois |
283 |
for (int i=0; i<3; i++) |
456 |
|
|
xyzduu[i]=skl(2,0)[i]; |
457 |
foucault |
27 |
|
458 |
francois |
283 |
for (int i=0; i<3; i++) |
459 |
|
|
xyzduv[i] = sens * skl(1,1)[i]; |
460 |
foucault |
27 |
|
461 |
francois |
283 |
for (int i=0; i<3; i++) |
462 |
|
|
xyzdvv[i]=skl(0,2)[i]; |
463 |
foucault |
27 |
|
464 |
francois |
283 |
if ((xyzdu!=NULL) && (xyzdv!=NULL ) ) |
465 |
|
|
{ |
466 |
|
|
for (int i=0; i<3; i++) |
467 |
|
|
xyzdu[i] = sens * skl(1,0)[i]; |
468 |
foucault |
27 |
|
469 |
francois |
283 |
for (int i=0; i<3; i++) |
470 |
|
|
xyzdv[i]=skl(0,1)[i]; |
471 |
|
|
} |
472 |
|
|
if (xyz!=NULL) |
473 |
|
|
{ |
474 |
|
|
for (int i=0; i<3; i++) |
475 |
|
|
xyz[i]=skl(0,0)[i]; |
476 |
|
|
} |
477 |
foucault |
27 |
|
478 |
francois |
283 |
unsigned kk, kkk; |
479 |
|
|
double LIMITE = 1E99; |
480 |
|
|
for (kk = 0; kk<3; kk++) |
481 |
foucault |
27 |
if ( (xyzduu[kk] > LIMITE ) || (xyzduv[kk] > LIMITE ) || (xyzdvv[kk] > LIMITE ) ) |
482 |
francois |
283 |
break; |
483 |
|
|
if (kk < 3) |
484 |
|
|
{ |
485 |
foucault |
27 |
deriver_seconde(uv,xyzduu,xyzduv,xyzdvv,xyz, xyzdu, xyzdv); |
486 |
|
|
printf("pb"); |
487 |
francois |
283 |
} |
488 |
foucault |
27 |
|
489 |
|
|
#undef skl |
490 |
|
|
} |
491 |
|
|
|
492 |
|
|
void ST_B_SPLINE_SURF::inverser(double *uv,double *xyz,double precision) |
493 |
|
|
{ |
494 |
francois |
283 |
int code=0; |
495 |
|
|
int num_point=(nb_ptsctr_u+nb_ptsctr_v)/2; |
496 |
|
|
do |
497 |
|
|
{ |
498 |
foucault |
27 |
code=inverser2(uv,xyz,num_point,precision); |
499 |
|
|
num_point=num_point*2; |
500 |
francois |
283 |
} |
501 |
|
|
while (code==0&&num_point<1500); |
502 |
|
|
} |
503 |
foucault |
27 |
|
504 |
|
|
|
505 |
|
|
int ST_B_SPLINE_SURF::inverser2(double *uv,double *xyz,int num_test, double precision) |
506 |
|
|
{ |
507 |
francois |
283 |
double uvi[2]; |
508 |
|
|
double xyz1[3]; |
509 |
|
|
double xyzdu[3]; |
510 |
|
|
double xyzdv[3]; |
511 |
|
|
double xyzduu[3]; |
512 |
|
|
double xyzduv[3]; |
513 |
|
|
double xyzdvv[3]; |
514 |
|
|
OT_VECTEUR_3D Pt(xyz[0],xyz[1],xyz[2]); |
515 |
|
|
double a[4]; |
516 |
|
|
double delta[2]; |
517 |
|
|
double b[9]; |
518 |
|
|
double delta_u; |
519 |
|
|
double delta_v; |
520 |
|
|
double ui; |
521 |
|
|
double vi; |
522 |
foucault |
27 |
|
523 |
francois |
283 |
double distance_ref=1e308; |
524 |
|
|
int refu; |
525 |
|
|
int refv; |
526 |
foucault |
27 |
|
527 |
francois |
283 |
for (int ii=0;ii<num_test+1;ii++) |
528 |
|
|
{ |
529 |
foucault |
27 |
double u=umin+ii*1./num_test*(umax-umin); |
530 |
|
|
uvi[0]=u; |
531 |
francois |
283 |
for (int jj=0;jj<num_test+1;jj++) |
532 |
|
|
{ |
533 |
|
|
double v=vmin+jj*1./num_test*(vmax-vmin); |
534 |
|
|
uvi[1]=v; |
535 |
|
|
evaluer(uvi,xyz1); |
536 |
|
|
OT_VECTEUR_3D S(xyz1[0],xyz1[1],xyz1[2]); |
537 |
|
|
OT_VECTEUR_3D Distance = S-Pt; |
538 |
|
|
double longueur=Distance.get_longueur2(); |
539 |
|
|
if (longueur<distance_ref) |
540 |
|
|
{ |
541 |
|
|
distance_ref=longueur; |
542 |
|
|
refu=ii; |
543 |
|
|
refv=jj; |
544 |
|
|
} |
545 |
foucault |
27 |
} |
546 |
francois |
283 |
} |
547 |
foucault |
27 |
|
548 |
francois |
283 |
double uii=umin+refu*1./num_test*(umax-umin); |
549 |
|
|
double vii=vmin+refv*1./num_test*(vmax-vmin); |
550 |
|
|
int compteur = 0; |
551 |
|
|
do |
552 |
|
|
{ |
553 |
foucault |
27 |
compteur++; |
554 |
|
|
ui=uii; |
555 |
|
|
vi=vii; |
556 |
|
|
uvi[0]=ui; |
557 |
|
|
uvi[1]=vi; |
558 |
|
|
deriver_seconde(uvi,xyzduu,xyzduv,xyzdvv,xyz1,xyzdu,xyzdv); |
559 |
|
|
OT_VECTEUR_3D S(xyz1[0],xyz1[1],xyz1[2]); |
560 |
|
|
OT_VECTEUR_3D Su(xyzdu[0],xyzdu[1],xyzdu[2]); |
561 |
|
|
OT_VECTEUR_3D Sv(xyzdv[0],xyzdv[1],xyzdv[2]); |
562 |
|
|
OT_VECTEUR_3D Suu(xyzduu[0],xyzduu[1],xyzduu[2]); |
563 |
|
|
OT_VECTEUR_3D Suv(xyzduv[0],xyzduv[1],xyzduv[2]); |
564 |
|
|
OT_VECTEUR_3D Svv(xyzdvv[0],xyzdvv[1],xyzdvv[2]); |
565 |
|
|
OT_VECTEUR_3D Distance = S-Pt; |
566 |
|
|
a[0]=Su.get_longueur2()+Distance*Suu; |
567 |
|
|
a[1]=Su*Sv+Distance*Suv; |
568 |
|
|
a[2]=Su*Sv+Distance*Suv; |
569 |
|
|
a[3]=Sv.get_longueur2()+Distance*Svv; |
570 |
francois |
283 |
b[0]=Distance*Su; |
571 |
|
|
b[0]=-b[0]; |
572 |
|
|
b[1]=Distance*Sv; |
573 |
|
|
b[1]=-b[1]; |
574 |
foucault |
27 |
double denominateur_delta=(a[0]*a[3]-a[2]*a[1]); |
575 |
francois |
35 |
if (fabs(denominateur_delta) < ( (fabs(a[0])+fabs(a[1])+fabs(a[2])+fabs(a[3]))*1e-12 ) ) |
576 |
foucault |
27 |
return 0; |
577 |
|
|
delta_u=(b[0]*a[3]-b[1]*a[1])/denominateur_delta; |
578 |
|
|
delta_v=(a[0]*b[1]-a[2]*b[0])/denominateur_delta; |
579 |
francois |
35 |
if (fabs(delta_u)>umax-umin) return 0; |
580 |
|
|
if (fabs(delta_v)>vmax-vmin) return 0; |
581 |
foucault |
27 |
if (Su.get_longueur2()>1E-14) |
582 |
|
|
uii=ui+delta_u; |
583 |
|
|
else |
584 |
|
|
uii=ui; |
585 |
|
|
if (Sv.get_longueur2()>1E-14) |
586 |
|
|
vii=vi+delta_v; |
587 |
|
|
else |
588 |
|
|
vii=vi; |
589 |
|
|
if (periodique_u==0) |
590 |
francois |
283 |
{ |
591 |
|
|
if (uii<umin) uii=umin; |
592 |
|
|
if (uii>umax) uii=umax; |
593 |
|
|
} |
594 |
foucault |
27 |
if (periodique_v==0) |
595 |
francois |
283 |
{ |
596 |
|
|
if (vii<vmin) vii=vmin; |
597 |
|
|
if (vii>vmax) vii=vmax; |
598 |
|
|
} |
599 |
foucault |
27 |
if (periodique_u==1) |
600 |
francois |
283 |
{ |
601 |
|
|
if (uii<umin) uii=umax - (umin-uii); |
602 |
|
|
if (uii>umax) uii=umin + (uii-umax); |
603 |
|
|
} |
604 |
foucault |
27 |
if (periodique_v==1) |
605 |
francois |
283 |
{ |
606 |
|
|
if (vii<vmin) vii=vmax - (vmin-vii); |
607 |
|
|
if (vii>vmax) vii=vmin + (vii-vmax); |
608 |
|
|
} |
609 |
foucault |
27 |
delta_u=uii-ui; |
610 |
|
|
delta_v=vii-vi; |
611 |
|
|
if (compteur>500) return 0; |
612 |
francois |
283 |
} |
613 |
foucault |
27 |
|
614 |
francois |
283 |
while ((fabs(delta_u)>precision)||(fabs(delta_v)>precision)); |
615 |
|
|
uv[0]=uii; |
616 |
|
|
uv[1]=vii; |
617 |
|
|
return 1; |
618 |
foucault |
27 |
} |
619 |
|
|
|
620 |
|
|
int ST_B_SPLINE_SURF::est_periodique_u(void) |
621 |
|
|
{ |
622 |
francois |
283 |
return periodique_u; |
623 |
foucault |
27 |
} |
624 |
|
|
int ST_B_SPLINE_SURF::est_periodique_v(void) |
625 |
|
|
{ |
626 |
francois |
283 |
return periodique_v; |
627 |
foucault |
27 |
} |
628 |
|
|
double ST_B_SPLINE_SURF::get_periode_u(void) |
629 |
|
|
{ |
630 |
francois |
283 |
return periode_u; |
631 |
foucault |
27 |
} |
632 |
|
|
double ST_B_SPLINE_SURF::get_periode_v(void) |
633 |
|
|
{ |
634 |
francois |
283 |
return periode_v; |
635 |
foucault |
27 |
} |
636 |
|
|
double ST_B_SPLINE_SURF::get_umin(void) |
637 |
|
|
{ |
638 |
francois |
283 |
return umin; |
639 |
foucault |
27 |
} |
640 |
|
|
double ST_B_SPLINE_SURF::get_umax(void) |
641 |
|
|
{ |
642 |
francois |
283 |
return umax; |
643 |
foucault |
27 |
} |
644 |
|
|
double ST_B_SPLINE_SURF::get_vmin(void) |
645 |
|
|
{ |
646 |
francois |
283 |
return vmin; |
647 |
foucault |
27 |
} |
648 |
|
|
double ST_B_SPLINE_SURF::get_vmax(void) |
649 |
|
|
{ |
650 |
francois |
283 |
return vmax; |
651 |
foucault |
27 |
} |
652 |
|
|
|
653 |
|
|
int ST_B_SPLINE_SURF::get_type_geometrique(TPL_LISTE_ENTITE<double> ¶m) |
654 |
|
|
{ |
655 |
francois |
283 |
for (int i=0;i<nb_ptsctr_u-(degre_u+1);i++) |
656 |
|
|
{ |
657 |
|
|
param.ajouter(knots_u[i]); |
658 |
|
|
} |
659 |
|
|
for (int i=0;i<nb_ptsctr_v-(degre_v+1);i++) |
660 |
|
|
{ |
661 |
|
|
param.ajouter(knots_v[i]); |
662 |
|
|
} |
663 |
|
|
double xyz[3]; |
664 |
|
|
for (int i=0;i<nb_point;i++) |
665 |
|
|
{ |
666 |
|
|
param.ajouter(ptsctr[i].x()/ptsctr[i].w()); |
667 |
|
|
param.ajouter(ptsctr[i].y()/ptsctr[i].w()); |
668 |
|
|
param.ajouter(ptsctr[i].z()/ptsctr[i].w()); |
669 |
|
|
} |
670 |
|
|
for (int i=0;i<nb_point;i++) |
671 |
|
|
{ |
672 |
|
|
param.ajouter(ptsctr[i].w()); |
673 |
|
|
} |
674 |
|
|
param.ajouter(degre_u); |
675 |
|
|
param.ajouter(degre_v); |
676 |
|
|
return MGCo_BSPLINES; |
677 |
foucault |
27 |
} |
678 |
|
|
|
679 |
|
|
void ST_B_SPLINE_SURF::initialiser(class ST_GESTIONNAIRE* gest) |
680 |
|
|
{ |
681 |
francois |
283 |
for (int i=0;i<nb_point;i++) |
682 |
|
|
{ |
683 |
foucault |
27 |
ST_POINT* stpoint=gest->lst_point.getid(indexptsctr[i]); |
684 |
|
|
double xyz[3]; |
685 |
|
|
stpoint->evaluer(xyz); |
686 |
|
|
double w=1.0; |
687 |
|
|
ptsctr.push_back(OT_VECTEUR_4D(w*xyz[0],w*xyz[1],w*xyz[2],w)); |
688 |
francois |
283 |
} |
689 |
foucault |
27 |
|
690 |
francois |
283 |
periodique_u=1; |
691 |
|
|
for (int i=0; i<nb_ptsctr_v; i++) |
692 |
|
|
{ |
693 |
foucault |
27 |
ST_POINT* point1_u=gest->lst_point.getid(indexptsctr[i]); |
694 |
|
|
ST_POINT* point2_u=gest->lst_point.getid(indexptsctr[nb_ptsctr_v*(nb_ptsctr_u-1)+i]); |
695 |
|
|
double xyz1[3]; |
696 |
|
|
double xyz2[3]; |
697 |
|
|
point1_u->evaluer(xyz1); |
698 |
|
|
point2_u->evaluer(xyz2); |
699 |
|
|
if ((!(OPERATEUR::egal (xyz1[0],xyz2[0],1E-10))) || (!(OPERATEUR::egal (xyz1[1],xyz2[1],1E-10))) || (!(OPERATEUR::egal (xyz1[2],xyz2[2],1E-10)))) |
700 |
francois |
283 |
{ |
701 |
|
|
periodique_u=0; |
702 |
foucault |
27 |
} |
703 |
francois |
283 |
} |
704 |
|
|
if (periodique_u==1) |
705 |
|
|
{ |
706 |
foucault |
27 |
int i=knots_u.size(); |
707 |
|
|
periode_u=(knots_u[i-1]-knots_u[0]); |
708 |
francois |
283 |
} |
709 |
|
|
else periode_u=0; |
710 |
|
|
periodique_v=1; |
711 |
|
|
for (int j=0; j<nb_ptsctr_u; j++) |
712 |
|
|
{ |
713 |
foucault |
27 |
ST_POINT* point1_v=gest->lst_point.getid(indexptsctr[j*nb_ptsctr_v]); |
714 |
|
|
ST_POINT* point2_v=gest->lst_point.getid(indexptsctr[j*nb_ptsctr_v+nb_ptsctr_v-1]); |
715 |
|
|
double xyz3[3]; |
716 |
|
|
double xyz4[3]; |
717 |
|
|
point1_v->evaluer(xyz3); |
718 |
|
|
point2_v->evaluer(xyz4); |
719 |
|
|
if ((!(OPERATEUR::egal (xyz3[0],xyz4[0],1E-10))) || (!(OPERATEUR::egal (xyz3[1],xyz4[1],1E-10))) || (!(OPERATEUR::egal (xyz3[2],xyz4[2],1E-10)))) |
720 |
francois |
283 |
{ |
721 |
|
|
periodique_v=0; |
722 |
foucault |
27 |
} |
723 |
francois |
283 |
} |
724 |
|
|
if (periodique_v==1) |
725 |
|
|
{ |
726 |
foucault |
27 |
int j=knots_v.size(); |
727 |
|
|
periode_v=(knots_v[j-1]-knots_v[0]); |
728 |
francois |
283 |
} |
729 |
|
|
else periode_v=0; |
730 |
foucault |
27 |
} |
731 |
|
|
|
732 |
|
|
|
733 |
|
|
|
734 |
|
|
void ST_B_SPLINE_SURF::est_util(ST_GESTIONNAIRE* gest) |
735 |
|
|
{ |
736 |
francois |
283 |
util=true; |
737 |
|
|
for (int i=0;i<nb_point;i++) |
738 |
foucault |
27 |
gest->lst_point.getid(indexptsctr[i])->est_util(gest); |
739 |
|
|
} |
740 |
|
|
|
741 |
|
|
|
742 |
|
|
|
743 |
|
|
void ST_B_SPLINE_SURF::get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> ¶m) |
744 |
|
|
{ |
745 |
|
|
|
746 |
|
|
|
747 |
|
|
// The first parameter indicate the code access |
748 |
francois |
283 |
param.ajouter(2); |
749 |
foucault |
27 |
// The follewing two parameters of the list indicate the orders of the net points |
750 |
|
|
|
751 |
francois |
283 |
param.ajouter(degre_u+1); |
752 |
|
|
param.ajouter(degre_v+1); |
753 |
foucault |
27 |
|
754 |
|
|
// The follewing two parameters indicate the number of rows and colons of the control points |
755 |
|
|
// respectively to the two parameters directions |
756 |
|
|
|
757 |
francois |
283 |
param.ajouter(nb_ptsctr_u); |
758 |
|
|
param.ajouter(nb_ptsctr_v); |
759 |
foucault |
27 |
|
760 |
|
|
// this present the knot vector in the u-direction |
761 |
|
|
|
762 |
francois |
283 |
for (unsigned int i=0;i<knots_u.size();i++) |
763 |
|
|
{ |
764 |
|
|
param.ajouter(knots_u[i]); |
765 |
|
|
} |
766 |
foucault |
27 |
|
767 |
|
|
//This present the knot vector in the v-direction |
768 |
|
|
|
769 |
francois |
283 |
for (unsigned int j=0;j<knots_v.size();j++) |
770 |
|
|
{ |
771 |
|
|
param.ajouter(knots_v[j]); |
772 |
|
|
} |
773 |
foucault |
27 |
|
774 |
francois |
283 |
for (int v=0;v<nb_ptsctr_v;v++) |
775 |
|
|
{ |
776 |
|
|
for (int u=0;u<nb_ptsctr_u;u++) |
777 |
|
|
{ |
778 |
|
|
int pt=u*nb_ptsctr_v+v; |
779 |
|
|
double w=ptsctr[pt].w(); |
780 |
|
|
double inv_w=1/w; |
781 |
|
|
double x=ptsctr[pt].x()*inv_w; |
782 |
|
|
double y=ptsctr[pt].y()*inv_w; |
783 |
|
|
double z=ptsctr[pt].z()*inv_w; |
784 |
|
|
param.ajouter(x); |
785 |
|
|
param.ajouter(y); |
786 |
|
|
param.ajouter(z); |
787 |
|
|
param.ajouter(w); |
788 |
|
|
} |
789 |
foucault |
27 |
|
790 |
francois |
283 |
} |
791 |
souaissa |
57 |
|
792 |
francois |
283 |
/* |
793 |
|
|
for(unsigned int pt=0;pt<ptsctr.size();pt++) |
794 |
|
|
{ |
795 |
foucault |
27 |
|
796 |
francois |
283 |
double w=ptsctr[pt].w(); |
797 |
|
|
double inv_w=1/w; |
798 |
|
|
double x=ptsctr[pt].x()*inv_w; |
799 |
|
|
double y=ptsctr[pt].y()*inv_w; |
800 |
|
|
double z=ptsctr[pt].z()*inv_w; |
801 |
|
|
param.ajouter(x); |
802 |
|
|
param.ajouter(y); |
803 |
|
|
param.ajouter(z); |
804 |
|
|
param.ajouter(w); |
805 |
foucault |
27 |
|
806 |
francois |
283 |
// OT_VECTEUR_4D ctrpoint=ptsctr[pt]; |
807 |
|
|
// for(int i=0;i<4;i++) |
808 |
|
|
// { |
809 |
|
|
// param.ajouter(ctrpoint[i]); |
810 |
|
|
// } |
811 |
foucault |
27 |
|
812 |
|
|
|
813 |
|
|
|
814 |
francois |
283 |
} */ |
815 |
foucault |
27 |
|
816 |
francois |
283 |
indx_premier_ptctr=5+knots_u.size()+knots_v.size(); |
817 |
foucault |
27 |
|
818 |
|
|
|
819 |
|
|
|
820 |
|
|
|
821 |
souaissa |
57 |
|
822 |
francois |
283 |
|
823 |
|
|
|
824 |
|
|
|
825 |
foucault |
27 |
} |