0027232: Configuration - fix mblen missing building issue on Android
[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-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17
18 #include <BSplCLib.hxx>
19 #include <Geom2d_BSplineCurve.hxx>
20 #include <Geom_BSplineCurve.hxx>
21 #include <gp_Pnt2d.hxx>
22 #include <Hermit.hxx>
23 #include <PLib.hxx>
24 #include <Precision.hxx>
25 #include <Standard_Boolean.hxx>
26 #include <Standard_DimensionError.hxx>
27 #include <Standard_Real.hxx>
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>
33
34 #include <algorithm>
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 static 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
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);
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
76 static 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
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);
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
111 static 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
128 static 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
158 static 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
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();
215
216   //sort of the array of knots
217   std::sort (knots.begin(), knots.end());
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())
231         Standard_DimensionError::Raise("Hermit Impossible Tolerance");       
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       }
290     }
291   } // end of the loop
292   
293   if (!SignDenom(Polesinit)) //invertion of the polynome sign
294   {
295     for (index=0;index<=3;index++)
296       Polesinit(index).SetCoord(0.0,-Polesinit(index).Y());
297   }
298
299   // loop of positivity
300   if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0))
301   {
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   
315   if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0))
316   {
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
379 static 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
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 
407   {
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     }
419   }
420   else
421   {
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();
439
440   //sort of the array of knots
441   std::sort (knots.begin(), knots.end());
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
448   // loop to check the tolerances on poles
449   if (TolPoles!=0.0)
450   {
451     Polemax(Polesinit,min,max);
452     Standard_Real  Polemin=Polesinit(min).Y();
453     Standard_Real  Polemax=Polesinit(max).Y();
454     if (((Polemax)>=((1/TolPoles)*Polemin))||((Polemin==0.0)&&(Polemax>=(1/TolPoles))))
455     {
456       if (Polesinit(0).Y()>=(1/TolPoles)*Polesinit(3).Y()||Polesinit(0).Y()<=TolPoles*Polesinit(3).Y())
457         Standard_DimensionError::Raise("Hermit Impossible Tolerance");       
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       }
523     }
524   } // end of the loop
525   
526   if (!SignDenom(Polesinit)) // invertion of the polynome sign
527   {
528     for (index=0;index<=3;index++)
529       Polesinit(index).SetCoord(0.0,-Polesinit(index).Y());
530   }
531
532   // boucle de positivite
533   if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0))
534   {
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   
548   if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0))
549   {
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  
611 static 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
628 static 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
647 Handle(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
717     Standard_DimensionError::Raise("Hermit Impossible Tolerance");
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
737 Handle(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
807     Standard_DimensionError::Raise("Hermit Impossible Tolerance");
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
827 void 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
899     Standard_DimensionError::Raise("Hermit Impossible Tolerance");
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