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