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