c4881e2a9eff8cc4b78b579f6453031e08705c08
[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         for (i=0;i<=3;i++)
241           Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
242       if ((max==1)||(max==2))
243         if ((min==0)||(min==3))                                             
244           for (i=0;i<=3;i++)
245             Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-(1/TolPoles)*Polemin));
246         else{                                                                  
247           if ((TolPoles*Polemax<Polesinit(0).Y())&&(TolPoles*Polemax<Polesinit(3).Y())){
248             for (i=0;i<=3;i++)                                             
249               Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
250             mark=1;
251           }
252           if ((1/TolPoles*Polemin>Polesinit(0).Y())&&(1/TolPoles*Polemin>Polesinit(3).Y())&&(mark==0)){
253             for (i=0;i<=3;i++)                                             
254               Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-1/TolPoles*Polemin));
255             mark=1;
256           }
257           if (mark==0){
258             Standard_Real Pole0,Pole3;
259             Pole0=Polesinit(0).Y();
260             Pole3=Polesinit(3).Y();
261             if (Pole0<3){                         
262               a=Log10(Pole3/Pole0);                      
263               if (boucle==2)                      
264                 for (i=0;i<=3;i++)                                                    
265                   Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); 
266               if (boucle==1){
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               a=Log10(Pole0/Pole3);                           
274               if (boucle==2)                      
275                 for (i=0;i<=3;i++)                                                    
276                   Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); 
277               if (boucle==1){
278                 for (i=0;i<=3;i++)                                                    
279                   Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles)))))); 
280                 dercas=1;
281               }
282             }
283           }
284         }
285     }
286   }                          //end of the loop
287   
288   if (!SignDenom(Polesinit))                       //invertion of the polynome sign
289     for (index=0;index<=3;index++)
290       Polesinit(index).SetCoord(0.0,-Polesinit(index).Y());
291   
292                                 //loop of positivity
293   if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0)){              
294     Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
295     if (boucle==2)
296       Us1=Us1*knots(2);
297     if (boucle==1)
298       if (Ux!=0.0)
299         Us1=Us1*Ux;
300     BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
301     if (I1<2)
302       U4=Us1;
303     else
304       U4=knots(I1);
305   }
306   
307   if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0)){
308     Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
309     if (boucle==2)
310       Us2=knots(knots.Length()-1)+Us2*(1-knots(knots.Length()-1));
311     if (boucle==1)
312       if (Ux!=0.0)
313         Us2=Uy+Us2*(1-Uy);
314     BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I1,Us2);
315     if (I1>=(knots.Length()-1))
316       U5=Us2;
317     else
318       U5=knots(I1+1);
319   }
320   
321   if (dercas==1)
322     boucle++;
323   
324   if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()<0.0)){
325     Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
326     Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
327     if (boucle!=0)
328       if (Ux!=0.0){
329         Us1=Us1*Ux;
330         Us2=Uy+Us2*(1-Uy);
331       }
332     if (Us2<=Us1){                                   
333       BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
334       if (knots(I1)>=Us2)                                    //insertion of one knot for the two poles
335         U4=knots(I1);
336       else{
337         if (I1>=2){                                          //insertion to the left and
338           U4=knots(I1);                                      //to the right without a new knot
339           BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2);
340           if (I3<(BS->NbKnots()-1)){
341             U5=knots(I3+1);
342             cas=1;
343           }
344         }
345         if(cas==0)                                          //insertion of only one new knot
346           U4=(Us1+Us2)/2;
347       }
348     }
349     else{                                                      //insertion of two knots
350       BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
351       if (I1>=2)                                                 
352         U4=knots(I1);
353       else
354         U4=Us1;
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       else
359         U5=Us2;
360     }
361   }
362
363
364 //=======================================================================
365 //function : PolyTest
366 //purpose  : give the knots U4 and U5 to insert to a(u)
367 //=======================================================================
368    
369
370 static void PolyTest(const TColStd_Array1OfReal&        Herm,
371                      const Handle(Geom2d_BSplineCurve)& BS,
372                      Standard_Real&                     U4,
373                      Standard_Real&                     U5,
374                      Standard_Integer&                  boucle,
375                      const Standard_Real                TolPoles,
376 //                   const Standard_Real                TolKnots,
377                      const Standard_Real                ,
378                      const Standard_Real                Ux, 
379                      const Standard_Real                Uy)
380
381 {
382   Standard_Integer               index,i,
383                                  I1=0,I2=0,I3=0,I4=0;    //knots index
384   TColgp_Array1OfPnt2d           Polesinit(0,3) ;
385   Handle(TColStd_HArray1OfReal)  Knots;                  //array of the BSpline knots + the ones inserted
386   Standard_Integer               cas=0,mark=0,dercas=0,  //loop marks
387                                  min,max;                //Pole min and max indices
388   Standard_Real                  Us1,Us2,a;              //bounderies value of the knots to be inserted
389 //  gp_Pnt2d                       P ;
390   TCollection_CompareOfReal      Comp;
391   
392   U4=0.0;U5=1.0;                                         //default value
393   if (Ux!=1.0){
394     BS->LocateU(Ux,0.0,I1,I2);                          //localization of the inserted knots
395     if (Uy!=0.0)
396       BS->LocateU(Uy,0.0,I3,I4);
397   }
398
399   if (I1==I2)                                            //definition and filling of the 
400     if((I3==I4)||(I3==0)){                               //array of knots
401       Knots=new TColStd_HArray1OfReal(1,BS->NbKnots());
402       for (i=1;i<=BS->NbKnots();i++)
403         Knots->SetValue(i,BS->Knot(i));
404     }
405     else{
406       Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1);
407       for (i=1;i<=BS->NbKnots();i++)
408         Knots->SetValue(i,BS->Knot(i));
409       Knots->SetValue(BS->NbKnots()+1,Uy);
410     }
411   else{
412     if((I3==I4)||(I3==0)){
413       Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1);
414       for (i=1;i<=BS->NbKnots();i++)
415         Knots->SetValue(i,BS->Knot(i));
416       Knots->SetValue(BS->NbKnots()+1,Ux);
417     }
418     else{
419       Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+2);
420       for (i=1;i<=BS->NbKnots();i++)
421         Knots->SetValue(i,BS->Knot(i));
422       Knots->SetValue(BS->NbKnots()+1,Ux);
423       Knots->SetValue(BS->NbKnots()+2,Uy);
424     }
425   }
426
427   TColStd_Array1OfReal knots(1,Knots->Length());
428   knots=Knots->ChangeArray1();
429   SortTools_QuickSortOfReal::Sort(knots,Comp);     //sort of the array of knots
430
431   Polesinit(0).SetCoord(0.0,Herm(0));              //poles of the Hermite polynome in the BSpline form
432   Polesinit(1).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
433   Polesinit(2).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
434   Polesinit(3).SetCoord(0.0,Herm(3));
435
436                               //loop to check the tolerances on poles
437   if (TolPoles!=0.0){
438     Polemax(Polesinit,min,max);
439     Standard_Real  Polemin=Polesinit(min).Y();
440     Standard_Real  Polemax=Polesinit(max).Y();
441     if (((Polemax)>=((1/TolPoles)*Polemin))||((Polemin==0.0)&&(Polemax>=(1/TolPoles)))){
442       if (Polesinit(0).Y()>=(1/TolPoles)*Polesinit(3).Y()||Polesinit(0).Y()<=TolPoles*Polesinit(3).Y())
443         Standard_DimensionError::Raise("Hermit Impossible Tolerance");       
444       if ((max==0)||(max==3))                                                
445         for (i=0;i<=3;i++)
446           Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
447       if ((max==1)||(max==2))
448         if ((min==0)||(min==3))                                             
449           for (i=0;i<=3;i++)
450             Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-(1/TolPoles)*Polemin));
451         else{                                                                  
452           if ((TolPoles*Polemax<Polesinit(0).Y())&&(TolPoles*Polemax<Polesinit(3).Y())){
453             for (i=0;i<=3;i++)                                             
454               Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
455             mark=1;
456           }
457           if ((1/TolPoles*Polemin>Polesinit(0).Y())&&(1/TolPoles*Polemin>Polesinit(3).Y())&&(mark==0)){
458             for (i=0;i<=3;i++)                                             
459               Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-1/TolPoles*Polemin));
460             mark=1;
461           }
462           if (mark==0){
463             Standard_Real Pole0,Pole3;
464             Pole0=Polesinit(0).Y();
465             Pole3=Polesinit(3).Y();
466             if (Pole0<3){                         
467               a=Log10(Pole3/Pole0);                      
468               if (boucle==2)                      
469                 for (i=0;i<=3;i++)                                                    
470                   Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); 
471               if (boucle==1){
472                 for (i=0;i<=3;i++)                                                    
473                   Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles)))))); 
474                 dercas=1;
475               }
476             }
477             if (Pole0>Pole3){                         
478               a=Log10(Pole0/Pole3);                           
479               if (boucle==2)                      
480                 for (i=0;i<=3;i++)                                                    
481                   Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0))))); 
482               if (boucle==1){
483                 for (i=0;i<=3;i++)                                                    
484                   Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles)))))); 
485                 dercas=1;
486               }
487             }
488           }
489         }
490     }
491   }                          //end of the loop
492   
493   if (!SignDenom(Polesinit))                       //invertion of the polynome sign
494     for (index=0;index<=3;index++)
495       Polesinit(index).SetCoord(0.0,-Polesinit(index).Y());
496   
497                                 //boucle de positivite
498   if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0)){
499     Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
500     if (boucle==2)
501       Us1=Us1*knots(2);
502     if (boucle==1)
503       if (Ux!=0.0)
504         Us1=Us1*Ux;
505     BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
506     if (I1<2)
507       U4=Us1;
508     else
509       U4=knots(I1);
510   }
511   
512   if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0)){
513     Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
514     if (boucle==2)
515       Us2=knots(knots.Length()-1)+Us2*(1-knots(knots.Length()-1));
516     if (boucle==1)
517       if (Ux!=0.0)
518         Us2=Uy+Us2*(1-Uy);
519     BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I1,Us2);
520     if (I1>=(knots.Length()-1))
521       U5=Us2;
522     else
523       U5=knots(I1+1);
524   }
525   
526   if (dercas==1)
527     boucle++;
528   
529   if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()<0.0)){
530     Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
531     Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
532     if (boucle!=0)
533       if (Ux!=0.0){
534         Us1=Us1*Ux;
535         Us2=Uy+Us2*(1-Uy);
536       }
537     if (Us2<=Us1){                                   
538       BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
539       if (knots(I1)>=Us2)                                    //insertion of one knot for the two poles
540         U4=knots(I1);
541       else{
542         if (I1>=2){                                          //insertion to the left and
543           U4=knots(I1);                                    //to the right without a new knot
544           BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2);
545           if (I3<(BS->NbKnots()-1)){
546             U5=knots(I3+1);
547             cas=1;
548           }
549         }
550         if(cas==0)                                          //insertion of only one new knot
551           U4=(Us1+Us2)/2;
552       }
553     }
554     else{                                                      //insertion of two knots
555       BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
556       if (I1>=2)                                                 
557         U4=knots(I1);
558       else
559         U4=Us1;
560       BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2);
561       if (I3<(BS->NbKnots()-1))
562         U5=knots(I3+1);
563       else
564         U5=Us2;
565     }
566   }
567
568                                
569 //=======================================================================
570 //function : InsertKnots
571 //purpose  : insert the knots in BS knot sequence if they are not null. 
572 //======================================================================= 
573  
574 static void InsertKnots(Handle(Geom2d_BSplineCurve)& BS,
575                         const Standard_Real          U4,
576                         const Standard_Real          U5)
577
578 {
579   if (U4!=0.0)                                  //insertion of :0 knot if U4=0
580     BS->InsertKnot(U4);                         //              1 knot if U4=U5 
581   if ((U5!=1.0)&&(U5!=U4))                      //              2 knots otherwise
582     BS->InsertKnot(U5);
583   
584 }
585
586 //=======================================================================
587 //function : MovePoles
588 //purpose  : move the poles above the x axis
589 //======================================================================= 
590
591 static void MovePoles(Handle(Geom2d_BSplineCurve)& BS)
592
593 {
594   gp_Pnt2d  P ;
595 //  Standard_Integer i,index; 
596   Standard_Integer i; 
597
598   for (i=3;i<=(BS->NbPoles()-2);i++){
599     P.SetCoord(1,(BS->Pole(i).Coord(1)));           //raising of the no constrained poles to
600     P.SetCoord(2,(BS->Pole(1).Coord(2)));           //the first pole level          
601     BS->SetPole(i,P);
602   }
603 }
604
605 //=======================================================================
606 //function : Solution
607 //purpose  :
608 //======================================================================= 
609
610 Handle(Geom2d_BSplineCurve) Hermit::Solution(const Handle(Geom_BSplineCurve)& BS,
611                                              const Standard_Real              TolPoles,
612                                              const Standard_Real              TolKnots)
613      
614 {
615   TColStd_Array1OfReal       Herm(0,3);
616   Standard_Real              Upos1=0.0, Upos2=1.0,   //positivity knots
617                              Ux=0.0,    Uy=1.0,
618                              Utol1=0.0, Utol2=1.0,   //tolerance knots
619                              Uint1=0.0, Uint2=1.0;   //tolerance knots for the first loop
620   Standard_Integer           boucle=1;               //loop mark
621   TColStd_Array1OfReal       Knots(1,2);             
622   TColStd_Array1OfInteger    Multiplicities(1,2);
623   TColgp_Array1OfPnt2d       Poles(1,4);
624   Standard_Integer           zeroboucle = 0 ;
625   HermiteCoeff(BS,Herm);                             //computation of the Hermite coefficient
626
627   Poles(1).SetCoord(0.0,Herm(0));                    //poles of the Hermite polynome in the BSpline form
628   Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
629   Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
630   Poles(4).SetCoord(0.0,Herm(3));
631   Knots(1)=0.0;
632   Knots(2)=1.0;
633   Multiplicities(1)=4;
634   Multiplicities(2)=4;
635
636   Handle(Geom2d_BSplineCurve)  BS1=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic
637   Handle(Geom2d_BSplineCurve)  BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//BSpline without modif
638
639   PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots
640   InsertKnots(BS2,Upos1,Upos2);                      //and insertion
641     
642   if (Upos1!=0.0)
643     if (Upos2!=1.0){
644       Ux=Min(Upos1,Upos2);
645       Uy=Max(Upos1,Upos2);
646     }
647     else{
648       Ux=Upos1;
649       Uy=Upos1;
650     }
651   else{
652     Ux=Upos2;
653     Uy=Upos2;
654   }
655
656   Herm(0)=BS2->Pole(1).Y();                           //computation of the Hermite coefficient on the
657   Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());      //positive BSpline
658   Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
659   Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
660   
661   PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy);          //computation of the tolerance knots
662   InsertKnots(BS2,Utol1,Utol2);                                          //and insertion
663   
664   if (boucle==2){                                     //insertion of two knots
665     Herm(0)=BS2->Pole(1).Y();
666     Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());
667     Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
668     Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
669     if (Utol1==0.0){
670       Uint2=Utol2;
671       PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0);
672     }
673     else{
674       Uint1=Utol1;
675       PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0);
676     }
677     InsertKnots(BS2,Utol1,Utol2);
678   }
679   if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance
680     Standard_DimensionError::Raise("Hermit Impossible Tolerance");
681   else{
682     if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0))    //test on the final inserted knots
683       InsertKnots(BS1,BS2->Knot(2),1.0);
684     else{
685       if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0))
686         InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),1.0);
687       else
688         InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),BS2->Knot(2));
689     }
690     MovePoles(BS1);                                  //relocation of the no-contrained knots
691   }
692   return BS1;
693 }
694
695
696 //=======================================================================
697 // Solution
698 //======================================================================= 
699
700 Handle(Geom2d_BSplineCurve) Hermit::Solution(const Handle(Geom2d_BSplineCurve)& BS,
701                                              const Standard_Real                TolPoles,
702                                              const Standard_Real                TolKnots)
703
704 {
705   TColStd_Array1OfReal       Herm(0,3);
706   Standard_Real              Upos1=0.0, Upos2=1.0,   //positivity knots
707                              Ux=0.0,    Uy=1.0, 
708                              Utol1=0.0, Utol2=1.0,   //tolerance knots
709                              Uint1=0.0, Uint2=1.0;   //tolerance knots for the first loop
710   Standard_Integer           boucle=1;               //loop mark
711   TColStd_Array1OfReal       Knots(1,2);             
712   TColStd_Array1OfInteger    Multiplicities(1,2);
713   TColgp_Array1OfPnt2d       Poles(1,4);
714   Standard_Integer           zeroboucle = 0 ;
715   HermiteCoeff(BS,Herm);                             //computation of the Hermite coefficient
716
717   Poles(1).SetCoord(0.0,Herm(0));                    //poles of the Hermite polynome in the BSpline form
718   Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
719   Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
720   Poles(4).SetCoord(0.0,Herm(3));
721   Knots(1)=0.0;
722   Knots(2)=1.0;
723   Multiplicities(1)=4;
724   Multiplicities(2)=4;
725
726   Handle(Geom2d_BSplineCurve)  BS1=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic
727   Handle(Geom2d_BSplineCurve)  BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//BSpline without modif
728
729   PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots
730   InsertKnots(BS2,Upos1,Upos2);                      //and insertion
731     
732   if (Upos1!=0.0)
733     if (Upos2!=1.0){
734       Ux=Min(Upos1,Upos2);
735       Uy=Max(Upos1,Upos2);
736     }
737     else{
738       Ux=Upos1;
739       Uy=Upos1;
740     }
741   else{
742     Ux=Upos2;
743     Uy=Upos2;
744   }
745   
746   Herm(0)=BS2->Pole(1).Y();                           //computation of the Hermite coefficient on the
747   Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());      //positive BSpline
748   Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
749   Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
750   
751   PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy);          //computation of the tolerance knots
752   InsertKnots(BS2,Utol1,Utol2);                                          //and insertion
753   
754   if (boucle==2){                                      //insertion of two knots
755     Herm(0)=BS2->Pole(1).Y();
756     Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());
757     Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
758     Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
759     if (Utol1==0.0){
760       Uint2=Utol2;
761       PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0);
762     }
763     else{
764       Uint1=Utol1;
765       PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0);
766     }
767     InsertKnots(BS2,Utol1,Utol2);
768   }
769   if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance
770     Standard_DimensionError::Raise("Hermit Impossible Tolerance");
771   else{
772     if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0))    //test on the final inserted knots
773       InsertKnots(BS1,BS2->Knot(2),1.0);
774     else{
775       if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0))
776         InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),1.0);
777       else
778         InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),BS2->Knot(2));
779     }
780     MovePoles(BS1);                          //relocation of the no-contrained knots
781   }
782   return BS1;
783 }
784
785 //=======================================================================
786 //function : Solutionbis
787 //purpose  :
788 //======================================================================= 
789
790 void Hermit::Solutionbis(const Handle(Geom_BSplineCurve)& BS,
791                          Standard_Real &                  Knotmin,
792                          Standard_Real &                  Knotmax,
793                          const Standard_Real              TolPoles,
794                          const Standard_Real              TolKnots)
795      
796 {
797   TColStd_Array1OfReal       Herm(0,3);
798   Standard_Real              Upos1=0.0, Upos2=1.0,   //positivity knots
799                              Ux=0.0,    Uy=1.0, 
800                              Utol1=0.0, Utol2=1.0,   //tolerance knots
801                              Uint1=0.0, Uint2=1.0;   //tolerance knots for the first loop
802   Standard_Integer           boucle=1;               //loop mark
803   TColStd_Array1OfReal       Knots(1,2);             
804   TColStd_Array1OfInteger    Multiplicities(1,2);
805   TColgp_Array1OfPnt2d       Poles(1,4);
806   Standard_Integer           zeroboucle = 0 ;
807   HermiteCoeff(BS,Herm);                             //computation of the Hermite coefficient
808   
809   Poles(1).SetCoord(0.0,Herm(0));                    //poles of the Hermite polynome in the BSpline form
810   Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
811   Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
812   Poles(4).SetCoord(0.0,Herm(3));
813   Knots(1)=0.0;
814   Knots(2)=1.0;
815   Multiplicities(1)=4;
816   Multiplicities(2)=4;
817   
818   Handle(Geom2d_BSplineCurve)  BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic
819                                                                                          //BSpline without modif
820   
821   PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots
822   InsertKnots(BS2,Upos1,Upos2);                      //and insertion
823   
824   if (Upos1!=0.0)
825     if (Upos2!=1.0){
826       Ux=Min(Upos1,Upos2);
827       Uy=Max(Upos1,Upos2);
828     }
829     else{
830       Ux=Upos1;
831       Uy=Upos1;
832     }
833   else{
834     Ux=Upos2;
835     Uy=Upos2;
836   }
837   
838   Herm(0)=BS2->Pole(1).Y();                           //computation of the Hermite coefficient on the
839   Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());      //positive BSpline
840   Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
841   Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
842   
843   PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy);          //computation of the tolerance knots
844   InsertKnots(BS2,Utol1,Utol2);                                          //and insertion
845   
846   if (boucle==2){                                      //insertion of two knots
847     Herm(0)=BS2->Pole(1).Y();
848     Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());
849     Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
850     Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
851     if (Utol1==0.0){
852       Uint2=Utol2;
853       PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0);
854     }
855     else{
856       Uint1=Utol1;
857       PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0);
858     }
859     InsertKnots(BS2,Utol1,Utol2);
860   }
861   if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance
862     Standard_DimensionError::Raise("Hermit Impossible Tolerance");
863   else{
864     if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0))    //test on the final inserted knots
865       Knotmin=BS2->Knot(2);
866     else{
867       if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0))
868         Knotmax=BS2->Knot(BS2->NbKnots()-1);
869       else{
870         Knotmin=BS2->Knot(2);
871         Knotmax=BS2->Knot(BS2->NbKnots()-1);
872       }
873     }  
874   }
875 }
876
877
878