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