b311480e |
1 | // Created on: 1997-01-15 |
2 | // Created by: Stagiaire Francois DUMONT |
3 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 |
4 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
5 | // |
973c2be1 |
6 | // This file is part of Open CASCADE Technology software library. |
b311480e |
7 | // |
973c2be1 |
8 | // This library is free software; you can redistribute it and / or modify it |
9 | // under the terms of the GNU Lesser General Public version 2.1 as published |
10 | // by the Free Software Foundation, with special exception defined in the file |
11 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
12 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
13 | // |
973c2be1 |
14 | // Alternatively, this file may be used under the terms of Open CASCADE |
15 | // commercial license or contractual agreement. |
7fd59977 |
16 | |
17 | #include <Hermit.ixx> |
18 | #include <Geom_BSplineCurve.hxx> |
19 | #include <Geom2d_BSplineCurve.hxx> |
20 | #include <BSplCLib.hxx> |
21 | #include <TColStd_Array1OfReal.hxx> |
22 | #include <TColStd_HArray1OfReal.hxx> |
23 | #include <TColStd_Array1OfInteger.hxx> |
24 | #include <TColStd_HArray1OfInteger.hxx> |
25 | #include <PLib.hxx> |
26 | #include <Standard_Boolean.hxx> |
27 | #include <gp_Pnt2d.hxx> |
28 | #include <TColgp_Array1OfPnt2d.hxx> |
29 | #include <Standard_DimensionError.hxx> |
30 | #include <Standard_Real.hxx> |
31 | #include <TCollection_CompareOfReal.hxx> |
32 | #include <SortTools_QuickSortOfReal.hxx> |
33 | #include <Precision.hxx> |
34 | |
35 | //======================================================================= |
36 | //function : HermiteCoeff |
37 | //purpose : calculate the Hermite coefficients of degree 3 from BS and |
38 | // store them in TAB(4 coefficients) |
39 | //======================================================================= |
40 | |
41 | static void HermiteCoeff(const Handle(Geom_BSplineCurve)& BS, |
42 | TColStd_Array1OfReal& TAB) |
43 | |
44 | { |
45 | TColStd_Array1OfReal Knots(1,BS->NbKnots()); |
46 | TColStd_Array1OfReal Weights(1,BS->NbPoles()); |
47 | TColStd_Array1OfInteger Mults(1,BS->NbKnots()); |
48 | Standard_Integer Degree,Index0,Index1; // denominateur value for u=0 & u=1 |
49 | Standard_Real Denom0,Denom1, // denominator value for u=0 & u=1 |
50 | Deriv0,Deriv1 ; // derivative denominator value for u=0 & 1 |
51 | Standard_Boolean Periodic; |
52 | |
53 | BS->Knots(Knots); |
54 | BSplCLib::Reparametrize(0.0,1.0,Knots); //affinity on the nodal vector |
55 | BS->Weights(Weights); |
56 | BS->Multiplicities(Mults); |
57 | Degree = BS->Degree(); |
58 | Periodic = BS->IsPeriodic(); |
59 | Index0 = BS->FirstUKnotIndex(); |
60 | Index1 = BS->LastUKnotIndex()-1; |
61 | |
62 | BSplCLib::D1(0.0,Index0,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom0,Deriv0); |
63 | BSplCLib::D1(1.0,Index1,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom1,Deriv1); |
64 | TAB(0) = 1/Denom0; //Hermit coefficients |
65 | TAB(1) = -Deriv0/(Denom0*Denom0); |
66 | TAB(2) = -Deriv1/(Denom1*Denom1); |
67 | TAB(3) = 1/Denom1; |
68 | |
69 | } |
70 | |
71 | //======================================================================= |
72 | //function : HermiteCoeff |
73 | //purpose : calculate the Hermite coefficients of degree 3 from BS and |
74 | // store them in TAB(4 coefficients) |
75 | //======================================================================= |
76 | |
77 | static void HermiteCoeff(const Handle(Geom2d_BSplineCurve)& BS, |
78 | TColStd_Array1OfReal& TAB) |
79 | |
80 | { |
81 | TColStd_Array1OfReal Knots(1,BS->NbKnots()); |
82 | TColStd_Array1OfReal Weights(1,BS->NbPoles()); |
83 | TColStd_Array1OfInteger Mults(1,BS->NbKnots()); |
84 | Standard_Integer Degree,Index0,Index1; |
85 | Standard_Real Denom0,Denom1, // denominateur value for u=0 & u=1 |
86 | Deriv0,Deriv1 ; // denominator value for u=0 & u=1 |
87 | Standard_Boolean Periodic; // derivative denominatur value for u=0 & 1 |
88 | |
89 | BS->Knots(Knots); |
90 | BSplCLib::Reparametrize(0.0,1.0,Knots); //affinity on the nodal vector |
91 | BS->Weights(Weights); |
92 | BS->Multiplicities(Mults); |
93 | Degree = BS->Degree(); |
94 | Periodic = BS->IsPeriodic(); |
95 | Index0 = BS->FirstUKnotIndex(); |
96 | Index1 = BS->LastUKnotIndex()-1; |
97 | |
98 | BSplCLib::D1(0.0,Index0,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom0,Deriv0); |
99 | BSplCLib::D1(1.0,Index1,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom1,Deriv1); |
100 | TAB(0) = 1/Denom0; //Hermit coefficients |
101 | TAB(1) = -Deriv0/(Denom0*Denom0); |
102 | TAB(2) = -Deriv1/(Denom1*Denom1); |
103 | TAB(3) = 1/Denom1; |
104 | |
105 | } |
106 | |
107 | //======================================================================= |
108 | //function : SignDenom |
109 | //purpose : give the sign of Herm(0) True=Positive |
110 | //======================================================================= |
111 | |
112 | static Standard_Boolean SignDenom(const TColgp_Array1OfPnt2d& Poles) |
113 | |
114 | { |
115 | Standard_Boolean Result; |
116 | |
117 | if (Poles(0).Y()<0) |
118 | Result=Standard_False; |
119 | else Result=Standard_True; |
120 | return Result; |
121 | } |
122 | |
123 | //======================================================================= |
124 | //function : Polemax |
125 | //purpose : give the max and the min of the Poles (by their index) |
126 | //======================================================================= |
127 | |
128 | |
129 | static void Polemax(const TColgp_Array1OfPnt2d& Poles, |
130 | Standard_Integer& min, |
131 | Standard_Integer& max) |
132 | |
133 | { |
134 | // Standard_Integer i,index=0; |
135 | Standard_Integer i; |
136 | Standard_Real Max,Min; //intermediate value of max and min ordinates |
137 | min=0;max=0; //initialisation of the indices |
138 | |
139 | Min=Poles(0).Y(); //initialisation of the intermediate value |
140 | Max=Poles(0).Y(); |
141 | for (i=1;i<=(Poles.Length()-1);i++){ |
142 | if (Poles(i).Y()<Min){ |
143 | Min=Poles(i).Y(); |
144 | min=i; |
145 | } |
146 | if (Poles(i).Y()>Max){ |
147 | Max=Poles(i).Y(); |
148 | max=i; |
149 | } |
150 | } |
151 | } |
152 | |
153 | //======================================================================= |
154 | //function : PolyTest |
155 | //purpose : give the knots U4 and U5 to insert to a(u) |
156 | //======================================================================= |
157 | |
158 | |
159 | static void PolyTest(const TColStd_Array1OfReal& Herm, |
160 | const Handle(Geom_BSplineCurve)& BS, |
161 | Standard_Real& U4, |
162 | Standard_Real& U5, |
163 | Standard_Integer& boucle, |
164 | const Standard_Real TolPoles, |
165 | // const Standard_Real TolKnots, |
166 | const Standard_Real , |
167 | const Standard_Real Ux, |
168 | const Standard_Real Uy) |
169 | |
170 | { |
171 | Standard_Integer index,i, |
172 | I1=0,I2=0,I3=0,I4=0; //knots index |
173 | TColgp_Array1OfPnt2d Polesinit(0,3) ; |
174 | Handle(TColStd_HArray1OfReal) Knots; //array of the BSpline knots + the ones inserted |
175 | Standard_Integer cas=0,mark=0,dercas=0, //loop marks |
176 | min,max; //Pole min and max indices |
177 | Standard_Real Us1,Us2,a; //bounderies value of the knots to be inserted |
178 | // gp_Pnt2d P ; |
179 | TCollection_CompareOfReal Comp; |
180 | |
181 | U4=0.0;U5=1.0; //default value |
182 | if (Ux!=1.0){ |
183 | BS->LocateU(Ux,0.0,I1,I2); //localization of the inserted knots |
184 | if (Uy!=0.0) |
185 | BS->LocateU(Uy,0.0,I3,I4); |
186 | } |
187 | |
188 | if (I1==I2) //definition and filling of the |
189 | if((I3==I4)||(I3==0)){ //array of knots |
190 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()); |
191 | for (i=1;i<=BS->NbKnots();i++) |
192 | Knots->SetValue(i,BS->Knot(i)); |
193 | } |
194 | else{ |
195 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1); |
196 | for (i=1;i<=BS->NbKnots();i++) |
197 | Knots->SetValue(i,BS->Knot(i)); |
198 | Knots->SetValue(BS->NbKnots()+1,Uy); |
199 | } |
200 | else{ |
201 | if((I3==I4)||(I3==0)){ |
202 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1); |
203 | for (i=1;i<=BS->NbKnots();i++) |
204 | Knots->SetValue(i,BS->Knot(i)); |
205 | Knots->SetValue(BS->NbKnots()+1,Ux); |
206 | } |
207 | else{ |
208 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+2); |
209 | for (i=1;i<=BS->NbKnots();i++) |
210 | Knots->SetValue(i,BS->Knot(i)); |
211 | Knots->SetValue(BS->NbKnots()+1,Ux); |
212 | Knots->SetValue(BS->NbKnots()+2,Uy); |
213 | } |
214 | } |
215 | |
216 | TColStd_Array1OfReal knots(1,Knots->Length()); |
217 | knots=Knots->ChangeArray1(); |
218 | SortTools_QuickSortOfReal::Sort(knots,Comp); //sort of the array of knots |
219 | |
220 | Polesinit(0).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form |
221 | Polesinit(1).SetCoord(0.0,Herm(0)+Herm(1)/3.0); |
222 | Polesinit(2).SetCoord(0.0,Herm(3)-Herm(2)/3.0); |
223 | Polesinit(3).SetCoord(0.0,Herm(3)); |
224 | |
225 | //loop to check the tolerances on poles |
226 | if (TolPoles!=0.0){ |
227 | Polemax(Polesinit,min,max); |
228 | Standard_Real Polemin=Polesinit(min).Y(); |
229 | Standard_Real Polemax=Polesinit(max).Y(); |
230 | if (((Polemax)>=((1/TolPoles)*Polemin))||((Polemin==0.0)&&(Polemax>=(1/TolPoles)))){ |
231 | if (Polesinit(0).Y()>=(1/TolPoles)*Polesinit(3).Y()||Polesinit(0).Y()<=TolPoles*Polesinit(3).Y()) |
eafb234b |
232 | Standard_DimensionError::Raise("Hermit Impossible Tolerance"); |
233 | if ((max==0)||(max==3)) |
234 | { |
235 | for (i=0;i<=3;i++) |
236 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax)); |
237 | } |
238 | else if ((max==1)||(max==2)) { |
239 | if ((min==0)||(min==3)) |
240 | { |
241 | for (i=0;i<=3;i++) |
242 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-(1/TolPoles)*Polemin)); |
243 | } |
244 | else{ |
245 | if ((TolPoles*Polemax<Polesinit(0).Y())&&(TolPoles*Polemax<Polesinit(3).Y())){ |
246 | for (i=0;i<=3;i++) |
247 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax)); |
248 | mark=1; |
249 | } |
250 | if ((1/TolPoles*Polemin>Polesinit(0).Y())&&(1/TolPoles*Polemin>Polesinit(3).Y())&&(mark==0)){ |
251 | for (i=0;i<=3;i++) |
252 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-1/TolPoles*Polemin)); |
253 | mark=1; |
254 | } |
255 | if (mark==0){ |
256 | Standard_Real Pole0,Pole3; |
257 | Pole0=Polesinit(0).Y(); |
258 | Pole3=Polesinit(3).Y(); |
259 | if (Pole0<3){ |
260 | a=Log10(Pole3/Pole0); |
261 | if (boucle==2) |
262 | { |
263 | for (i=0;i<=3;i++) |
264 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); |
265 | } |
266 | if (boucle==1) |
267 | { |
268 | for (i=0;i<=3;i++) |
269 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles)))))); |
270 | dercas=1; |
271 | } |
272 | } |
273 | if (Pole0>Pole3) |
274 | { |
275 | a=Log10(Pole0/Pole3); |
276 | if (boucle==2) |
277 | { |
278 | for (i=0;i<=3;i++) |
279 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); |
280 | } |
281 | if (boucle==1) |
282 | { |
283 | for (i=0;i<=3;i++) |
284 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles)))))); |
285 | dercas=1; |
286 | } |
287 | } |
288 | } |
289 | } |
290 | } |
7fd59977 |
291 | } |
eafb234b |
292 | } // end of the loop |
7fd59977 |
293 | |
eafb234b |
294 | if (!SignDenom(Polesinit)) //invertion of the polynome sign |
295 | { |
7fd59977 |
296 | for (index=0;index<=3;index++) |
297 | Polesinit(index).SetCoord(0.0,-Polesinit(index).Y()); |
eafb234b |
298 | } |
299 | |
300 | // loop of positivity |
301 | if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0)) |
302 | { |
7fd59977 |
303 | Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y()); |
304 | if (boucle==2) |
305 | Us1=Us1*knots(2); |
306 | if (boucle==1) |
307 | if (Ux!=0.0) |
308 | Us1=Us1*Ux; |
309 | BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1); |
310 | if (I1<2) |
311 | U4=Us1; |
312 | else |
313 | U4=knots(I1); |
314 | } |
315 | |
eafb234b |
316 | if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0)) |
317 | { |
7fd59977 |
318 | Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y()); |
319 | if (boucle==2) |
320 | Us2=knots(knots.Length()-1)+Us2*(1-knots(knots.Length()-1)); |
321 | if (boucle==1) |
322 | if (Ux!=0.0) |
323 | Us2=Uy+Us2*(1-Uy); |
324 | BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I1,Us2); |
325 | if (I1>=(knots.Length()-1)) |
326 | U5=Us2; |
327 | else |
328 | U5=knots(I1+1); |
329 | } |
330 | |
331 | if (dercas==1) |
332 | boucle++; |
333 | |
334 | if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()<0.0)){ |
335 | Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y()); |
336 | Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y()); |
337 | if (boucle!=0) |
338 | if (Ux!=0.0){ |
339 | Us1=Us1*Ux; |
340 | Us2=Uy+Us2*(1-Uy); |
341 | } |
342 | if (Us2<=Us1){ |
343 | BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1); |
344 | if (knots(I1)>=Us2) //insertion of one knot for the two poles |
345 | U4=knots(I1); |
346 | else{ |
347 | if (I1>=2){ //insertion to the left and |
348 | U4=knots(I1); //to the right without a new knot |
349 | BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2); |
350 | if (I3<(BS->NbKnots()-1)){ |
351 | U5=knots(I3+1); |
352 | cas=1; |
353 | } |
354 | } |
355 | if(cas==0) //insertion of only one new knot |
356 | U4=(Us1+Us2)/2; |
357 | } |
358 | } |
359 | else{ //insertion of two knots |
360 | BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1); |
361 | if (I1>=2) |
362 | U4=knots(I1); |
363 | else |
364 | U4=Us1; |
365 | BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2); |
366 | if (I3<(BS->NbKnots()-1)) |
367 | U5=knots(I3+1); |
368 | else |
369 | U5=Us2; |
370 | } |
371 | } |
372 | } |
373 | |
374 | //======================================================================= |
375 | //function : PolyTest |
376 | //purpose : give the knots U4 and U5 to insert to a(u) |
377 | //======================================================================= |
378 | |
379 | |
380 | static void PolyTest(const TColStd_Array1OfReal& Herm, |
381 | const Handle(Geom2d_BSplineCurve)& BS, |
382 | Standard_Real& U4, |
383 | Standard_Real& U5, |
384 | Standard_Integer& boucle, |
385 | const Standard_Real TolPoles, |
386 | // const Standard_Real TolKnots, |
387 | const Standard_Real , |
388 | const Standard_Real Ux, |
389 | const Standard_Real Uy) |
390 | |
391 | { |
392 | Standard_Integer index,i, |
393 | I1=0,I2=0,I3=0,I4=0; //knots index |
394 | TColgp_Array1OfPnt2d Polesinit(0,3) ; |
395 | Handle(TColStd_HArray1OfReal) Knots; //array of the BSpline knots + the ones inserted |
396 | Standard_Integer cas=0,mark=0,dercas=0, //loop marks |
397 | min,max; //Pole min and max indices |
398 | Standard_Real Us1,Us2,a; //bounderies value of the knots to be inserted |
399 | // gp_Pnt2d P ; |
400 | TCollection_CompareOfReal Comp; |
401 | |
402 | U4=0.0;U5=1.0; //default value |
403 | if (Ux!=1.0){ |
404 | BS->LocateU(Ux,0.0,I1,I2); //localization of the inserted knots |
405 | if (Uy!=0.0) |
406 | BS->LocateU(Uy,0.0,I3,I4); |
407 | } |
408 | |
409 | if (I1==I2) //definition and filling of the |
eafb234b |
410 | { |
7fd59977 |
411 | if((I3==I4)||(I3==0)){ //array of knots |
412 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()); |
413 | for (i=1;i<=BS->NbKnots();i++) |
414 | Knots->SetValue(i,BS->Knot(i)); |
415 | } |
416 | else{ |
417 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1); |
418 | for (i=1;i<=BS->NbKnots();i++) |
419 | Knots->SetValue(i,BS->Knot(i)); |
420 | Knots->SetValue(BS->NbKnots()+1,Uy); |
421 | } |
eafb234b |
422 | } |
423 | else |
424 | { |
7fd59977 |
425 | if((I3==I4)||(I3==0)){ |
426 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1); |
427 | for (i=1;i<=BS->NbKnots();i++) |
428 | Knots->SetValue(i,BS->Knot(i)); |
429 | Knots->SetValue(BS->NbKnots()+1,Ux); |
430 | } |
431 | else{ |
432 | Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+2); |
433 | for (i=1;i<=BS->NbKnots();i++) |
434 | Knots->SetValue(i,BS->Knot(i)); |
435 | Knots->SetValue(BS->NbKnots()+1,Ux); |
436 | Knots->SetValue(BS->NbKnots()+2,Uy); |
437 | } |
438 | } |
439 | |
440 | TColStd_Array1OfReal knots(1,Knots->Length()); |
441 | knots=Knots->ChangeArray1(); |
442 | SortTools_QuickSortOfReal::Sort(knots,Comp); //sort of the array of knots |
443 | |
444 | Polesinit(0).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form |
445 | Polesinit(1).SetCoord(0.0,Herm(0)+Herm(1)/3.0); |
446 | Polesinit(2).SetCoord(0.0,Herm(3)-Herm(2)/3.0); |
447 | Polesinit(3).SetCoord(0.0,Herm(3)); |
448 | |
eafb234b |
449 | // loop to check the tolerances on poles |
450 | if (TolPoles!=0.0) |
451 | { |
7fd59977 |
452 | Polemax(Polesinit,min,max); |
453 | Standard_Real Polemin=Polesinit(min).Y(); |
454 | Standard_Real Polemax=Polesinit(max).Y(); |
eafb234b |
455 | if (((Polemax)>=((1/TolPoles)*Polemin))||((Polemin==0.0)&&(Polemax>=(1/TolPoles)))) |
456 | { |
7fd59977 |
457 | if (Polesinit(0).Y()>=(1/TolPoles)*Polesinit(3).Y()||Polesinit(0).Y()<=TolPoles*Polesinit(3).Y()) |
eafb234b |
458 | Standard_DimensionError::Raise("Hermit Impossible Tolerance"); |
459 | if ((max==0)||(max==3)) |
460 | { |
461 | for (i=0;i<=3;i++) |
462 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax)); |
463 | } |
464 | else if ((max==1)||(max==2)) |
465 | { |
466 | if ((min==0)||(min==3)) |
467 | { |
468 | for (i=0;i<=3;i++) |
469 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-(1/TolPoles)*Polemin)); |
470 | } |
471 | else |
472 | { |
473 | if ((TolPoles*Polemax<Polesinit(0).Y())&&(TolPoles*Polemax<Polesinit(3).Y())) |
474 | { |
475 | for (i=0;i<=3;i++) |
476 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax)); |
477 | mark=1; |
478 | } |
479 | |
480 | if ((1/TolPoles*Polemin>Polesinit(0).Y())&&(1/TolPoles*Polemin>Polesinit(3).Y())&&(mark==0)) |
481 | { |
482 | for (i=0;i<=3;i++) |
483 | Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-1/TolPoles*Polemin)); |
484 | mark=1; |
485 | } |
486 | if (mark==0) |
487 | { |
488 | Standard_Real Pole0,Pole3; |
489 | Pole0=Polesinit(0).Y(); |
490 | Pole3=Polesinit(3).Y(); |
491 | if (Pole0<3) |
492 | { |
493 | a=Log10(Pole3/Pole0); |
494 | if (boucle==2) |
495 | { |
496 | for (i=0;i<=3;i++) |
497 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); |
498 | } |
499 | if (boucle==1) |
500 | { |
501 | for (i=0;i<=3;i++) |
502 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles)))))); |
503 | dercas=1; |
504 | } |
505 | } |
506 | if (Pole0>Pole3) |
507 | { |
508 | a=Log10(Pole0/Pole3); |
509 | if (boucle==2) |
510 | { |
511 | for (i=0;i<=3;i++) |
512 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); |
513 | } |
514 | else if (boucle==1) |
515 | { |
516 | for (i=0;i<=3;i++) |
517 | Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles)))))); |
518 | dercas=1; |
519 | } |
520 | } |
521 | } |
522 | } |
523 | } |
7fd59977 |
524 | } |
eafb234b |
525 | } // end of the loop |
7fd59977 |
526 | |
eafb234b |
527 | if (!SignDenom(Polesinit)) // invertion of the polynome sign |
528 | { |
7fd59977 |
529 | for (index=0;index<=3;index++) |
530 | Polesinit(index).SetCoord(0.0,-Polesinit(index).Y()); |
eafb234b |
531 | } |
532 | |
533 | // boucle de positivite |
534 | if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0)) |
535 | { |
7fd59977 |
536 | Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y()); |
537 | if (boucle==2) |
538 | Us1=Us1*knots(2); |
539 | if (boucle==1) |
540 | if (Ux!=0.0) |
541 | Us1=Us1*Ux; |
542 | BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1); |
543 | if (I1<2) |
544 | U4=Us1; |
545 | else |
546 | U4=knots(I1); |
547 | } |
548 | |
eafb234b |
549 | if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0)) |
550 | { |
7fd59977 |
551 | Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y()); |
552 | if (boucle==2) |
553 | Us2=knots(knots.Length()-1)+Us2*(1-knots(knots.Length()-1)); |
554 | if (boucle==1) |
555 | if (Ux!=0.0) |
556 | Us2=Uy+Us2*(1-Uy); |
557 | BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I1,Us2); |
558 | if (I1>=(knots.Length()-1)) |
559 | U5=Us2; |
560 | else |
561 | U5=knots(I1+1); |
562 | } |
563 | |
564 | if (dercas==1) |
565 | boucle++; |
566 | |
567 | if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()<0.0)){ |
568 | Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y()); |
569 | Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y()); |
570 | if (boucle!=0) |
571 | if (Ux!=0.0){ |
572 | Us1=Us1*Ux; |
573 | Us2=Uy+Us2*(1-Uy); |
574 | } |
575 | if (Us2<=Us1){ |
576 | BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1); |
577 | if (knots(I1)>=Us2) //insertion of one knot for the two poles |
578 | U4=knots(I1); |
579 | else{ |
580 | if (I1>=2){ //insertion to the left and |
581 | U4=knots(I1); //to the right without a new knot |
582 | BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2); |
583 | if (I3<(BS->NbKnots()-1)){ |
584 | U5=knots(I3+1); |
585 | cas=1; |
586 | } |
587 | } |
588 | if(cas==0) //insertion of only one new knot |
589 | U4=(Us1+Us2)/2; |
590 | } |
591 | } |
592 | else{ //insertion of two knots |
593 | BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1); |
594 | if (I1>=2) |
595 | U4=knots(I1); |
596 | else |
597 | U4=Us1; |
598 | BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2); |
599 | if (I3<(BS->NbKnots()-1)) |
600 | U5=knots(I3+1); |
601 | else |
602 | U5=Us2; |
603 | } |
604 | } |
605 | } |
606 | |
607 | //======================================================================= |
608 | //function : InsertKnots |
609 | //purpose : insert the knots in BS knot sequence if they are not null. |
610 | //======================================================================= |
611 | |
612 | static void InsertKnots(Handle(Geom2d_BSplineCurve)& BS, |
613 | const Standard_Real U4, |
614 | const Standard_Real U5) |
615 | |
616 | { |
617 | if (U4!=0.0) //insertion of :0 knot if U4=0 |
618 | BS->InsertKnot(U4); // 1 knot if U4=U5 |
619 | if ((U5!=1.0)&&(U5!=U4)) // 2 knots otherwise |
620 | BS->InsertKnot(U5); |
621 | |
622 | } |
623 | |
624 | //======================================================================= |
625 | //function : MovePoles |
626 | //purpose : move the poles above the x axis |
627 | //======================================================================= |
628 | |
629 | static void MovePoles(Handle(Geom2d_BSplineCurve)& BS) |
630 | |
631 | { |
632 | gp_Pnt2d P ; |
633 | // Standard_Integer i,index; |
634 | Standard_Integer i; |
635 | |
636 | for (i=3;i<=(BS->NbPoles()-2);i++){ |
637 | P.SetCoord(1,(BS->Pole(i).Coord(1))); //raising of the no constrained poles to |
638 | P.SetCoord(2,(BS->Pole(1).Coord(2))); //the first pole level |
639 | BS->SetPole(i,P); |
640 | } |
641 | } |
642 | |
643 | //======================================================================= |
644 | //function : Solution |
645 | //purpose : |
646 | //======================================================================= |
647 | |
648 | Handle(Geom2d_BSplineCurve) Hermit::Solution(const Handle(Geom_BSplineCurve)& BS, |
649 | const Standard_Real TolPoles, |
650 | const Standard_Real TolKnots) |
651 | |
652 | { |
653 | TColStd_Array1OfReal Herm(0,3); |
654 | Standard_Real Upos1=0.0, Upos2=1.0, //positivity knots |
655 | Ux=0.0, Uy=1.0, |
656 | Utol1=0.0, Utol2=1.0, //tolerance knots |
657 | Uint1=0.0, Uint2=1.0; //tolerance knots for the first loop |
658 | Standard_Integer boucle=1; //loop mark |
659 | TColStd_Array1OfReal Knots(1,2); |
660 | TColStd_Array1OfInteger Multiplicities(1,2); |
661 | TColgp_Array1OfPnt2d Poles(1,4); |
662 | Standard_Integer zeroboucle = 0 ; |
663 | HermiteCoeff(BS,Herm); //computation of the Hermite coefficient |
664 | |
665 | Poles(1).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form |
666 | Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0); |
667 | Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0); |
668 | Poles(4).SetCoord(0.0,Herm(3)); |
669 | Knots(1)=0.0; |
670 | Knots(2)=1.0; |
671 | Multiplicities(1)=4; |
672 | Multiplicities(2)=4; |
673 | |
674 | Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic |
675 | Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//BSpline without modif |
676 | |
677 | PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots |
678 | InsertKnots(BS2,Upos1,Upos2); //and insertion |
679 | |
680 | if (Upos1!=0.0) |
681 | if (Upos2!=1.0){ |
682 | Ux=Min(Upos1,Upos2); |
683 | Uy=Max(Upos1,Upos2); |
684 | } |
685 | else{ |
686 | Ux=Upos1; |
687 | Uy=Upos1; |
688 | } |
689 | else{ |
690 | Ux=Upos2; |
691 | Uy=Upos2; |
692 | } |
693 | |
694 | Herm(0)=BS2->Pole(1).Y(); //computation of the Hermite coefficient on the |
695 | Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); //positive BSpline |
696 | Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y()); |
697 | Herm(3)=BS2->Pole(BS2->NbPoles()).Y(); |
698 | |
699 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy); //computation of the tolerance knots |
700 | InsertKnots(BS2,Utol1,Utol2); //and insertion |
701 | |
702 | if (boucle==2){ //insertion of two knots |
703 | Herm(0)=BS2->Pole(1).Y(); |
704 | Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); |
705 | Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y()); |
706 | Herm(3)=BS2->Pole(BS2->NbPoles()).Y(); |
707 | if (Utol1==0.0){ |
708 | Uint2=Utol2; |
709 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0); |
710 | } |
711 | else{ |
712 | Uint1=Utol1; |
713 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0); |
714 | } |
715 | InsertKnots(BS2,Utol1,Utol2); |
716 | } |
717 | if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance |
718 | Standard_DimensionError::Raise("Hermit Impossible Tolerance"); |
719 | else{ |
720 | if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0)) //test on the final inserted knots |
721 | InsertKnots(BS1,BS2->Knot(2),1.0); |
722 | else{ |
723 | if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0)) |
724 | InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),1.0); |
725 | else |
726 | InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),BS2->Knot(2)); |
727 | } |
728 | MovePoles(BS1); //relocation of the no-contrained knots |
729 | } |
730 | return BS1; |
731 | } |
732 | |
733 | |
734 | //======================================================================= |
735 | // Solution |
736 | //======================================================================= |
737 | |
738 | Handle(Geom2d_BSplineCurve) Hermit::Solution(const Handle(Geom2d_BSplineCurve)& BS, |
739 | const Standard_Real TolPoles, |
740 | const Standard_Real TolKnots) |
741 | |
742 | { |
743 | TColStd_Array1OfReal Herm(0,3); |
744 | Standard_Real Upos1=0.0, Upos2=1.0, //positivity knots |
745 | Ux=0.0, Uy=1.0, |
746 | Utol1=0.0, Utol2=1.0, //tolerance knots |
747 | Uint1=0.0, Uint2=1.0; //tolerance knots for the first loop |
748 | Standard_Integer boucle=1; //loop mark |
749 | TColStd_Array1OfReal Knots(1,2); |
750 | TColStd_Array1OfInteger Multiplicities(1,2); |
751 | TColgp_Array1OfPnt2d Poles(1,4); |
752 | Standard_Integer zeroboucle = 0 ; |
753 | HermiteCoeff(BS,Herm); //computation of the Hermite coefficient |
754 | |
755 | Poles(1).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form |
756 | Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0); |
757 | Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0); |
758 | Poles(4).SetCoord(0.0,Herm(3)); |
759 | Knots(1)=0.0; |
760 | Knots(2)=1.0; |
761 | Multiplicities(1)=4; |
762 | Multiplicities(2)=4; |
763 | |
764 | Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic |
765 | Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//BSpline without modif |
766 | |
767 | PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots |
768 | InsertKnots(BS2,Upos1,Upos2); //and insertion |
769 | |
770 | if (Upos1!=0.0) |
771 | if (Upos2!=1.0){ |
772 | Ux=Min(Upos1,Upos2); |
773 | Uy=Max(Upos1,Upos2); |
774 | } |
775 | else{ |
776 | Ux=Upos1; |
777 | Uy=Upos1; |
778 | } |
779 | else{ |
780 | Ux=Upos2; |
781 | Uy=Upos2; |
782 | } |
783 | |
784 | Herm(0)=BS2->Pole(1).Y(); //computation of the Hermite coefficient on the |
785 | Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); //positive BSpline |
786 | Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y()); |
787 | Herm(3)=BS2->Pole(BS2->NbPoles()).Y(); |
788 | |
789 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy); //computation of the tolerance knots |
790 | InsertKnots(BS2,Utol1,Utol2); //and insertion |
791 | |
792 | if (boucle==2){ //insertion of two knots |
793 | Herm(0)=BS2->Pole(1).Y(); |
794 | Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); |
795 | Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y()); |
796 | Herm(3)=BS2->Pole(BS2->NbPoles()).Y(); |
797 | if (Utol1==0.0){ |
798 | Uint2=Utol2; |
799 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0); |
800 | } |
801 | else{ |
802 | Uint1=Utol1; |
803 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0); |
804 | } |
805 | InsertKnots(BS2,Utol1,Utol2); |
806 | } |
807 | if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance |
808 | Standard_DimensionError::Raise("Hermit Impossible Tolerance"); |
809 | else{ |
810 | if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0)) //test on the final inserted knots |
811 | InsertKnots(BS1,BS2->Knot(2),1.0); |
812 | else{ |
813 | if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0)) |
814 | InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),1.0); |
815 | else |
816 | InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),BS2->Knot(2)); |
817 | } |
818 | MovePoles(BS1); //relocation of the no-contrained knots |
819 | } |
820 | return BS1; |
821 | } |
822 | |
823 | //======================================================================= |
824 | //function : Solutionbis |
825 | //purpose : |
826 | //======================================================================= |
827 | |
828 | void Hermit::Solutionbis(const Handle(Geom_BSplineCurve)& BS, |
829 | Standard_Real & Knotmin, |
830 | Standard_Real & Knotmax, |
831 | const Standard_Real TolPoles, |
832 | const Standard_Real TolKnots) |
833 | |
834 | { |
835 | TColStd_Array1OfReal Herm(0,3); |
836 | Standard_Real Upos1=0.0, Upos2=1.0, //positivity knots |
837 | Ux=0.0, Uy=1.0, |
838 | Utol1=0.0, Utol2=1.0, //tolerance knots |
839 | Uint1=0.0, Uint2=1.0; //tolerance knots for the first loop |
840 | Standard_Integer boucle=1; //loop mark |
841 | TColStd_Array1OfReal Knots(1,2); |
842 | TColStd_Array1OfInteger Multiplicities(1,2); |
843 | TColgp_Array1OfPnt2d Poles(1,4); |
844 | Standard_Integer zeroboucle = 0 ; |
845 | HermiteCoeff(BS,Herm); //computation of the Hermite coefficient |
846 | |
847 | Poles(1).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form |
848 | Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0); |
849 | Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0); |
850 | Poles(4).SetCoord(0.0,Herm(3)); |
851 | Knots(1)=0.0; |
852 | Knots(2)=1.0; |
853 | Multiplicities(1)=4; |
854 | Multiplicities(2)=4; |
855 | |
856 | Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic |
857 | //BSpline without modif |
858 | |
859 | PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots |
860 | InsertKnots(BS2,Upos1,Upos2); //and insertion |
861 | |
862 | if (Upos1!=0.0) |
863 | if (Upos2!=1.0){ |
864 | Ux=Min(Upos1,Upos2); |
865 | Uy=Max(Upos1,Upos2); |
866 | } |
867 | else{ |
868 | Ux=Upos1; |
869 | Uy=Upos1; |
870 | } |
871 | else{ |
872 | Ux=Upos2; |
873 | Uy=Upos2; |
874 | } |
875 | |
876 | Herm(0)=BS2->Pole(1).Y(); //computation of the Hermite coefficient on the |
877 | Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); //positive BSpline |
878 | Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y()); |
879 | Herm(3)=BS2->Pole(BS2->NbPoles()).Y(); |
880 | |
881 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy); //computation of the tolerance knots |
882 | InsertKnots(BS2,Utol1,Utol2); //and insertion |
883 | |
884 | if (boucle==2){ //insertion of two knots |
885 | Herm(0)=BS2->Pole(1).Y(); |
886 | Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); |
887 | Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y()); |
888 | Herm(3)=BS2->Pole(BS2->NbPoles()).Y(); |
889 | if (Utol1==0.0){ |
890 | Uint2=Utol2; |
891 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0); |
892 | } |
893 | else{ |
894 | Uint1=Utol1; |
895 | PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0); |
896 | } |
897 | InsertKnots(BS2,Utol1,Utol2); |
898 | } |
899 | if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance |
900 | Standard_DimensionError::Raise("Hermit Impossible Tolerance"); |
901 | else{ |
902 | if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0)) //test on the final inserted knots |
903 | Knotmin=BS2->Knot(2); |
904 | else{ |
905 | if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0)) |
906 | Knotmax=BS2->Knot(BS2->NbKnots()-1); |
907 | else{ |
908 | Knotmin=BS2->Knot(2); |
909 | Knotmax=BS2->Knot(BS2->NbKnots()-1); |
910 | } |
911 | } |
912 | } |
913 | } |
914 | |
915 | |
916 | |