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