src/ProjLib/ProjLib_ComputeApprox.cxx

   1 // Created on: 1993-09-07
   2 // Created by: Bruno DUMORTIER
   3 // Copyright (c) 1993-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 // modified by NIZHNY-OFV  Thu Jan 20 11:04:19 2005
  18
  19 #include <ProjLib_ComputeApprox.hxx>
  20
  21 #include <GeomAbs_SurfaceType.hxx>
  22 #include <GeomAbs_CurveType.hxx>
  23 #include <AppCont_Function.hxx>
  24 #include <Convert_CompBezierCurves2dToBSplineCurve2d.hxx>
  25 #include <ElSLib.hxx>
  26 #include <ElCLib.hxx>
  27 #include <BSplCLib.hxx>
  28 #include <Standard_NoSuchObject.hxx>
  29 #include <Geom_UndefinedDerivative.hxx>
  30 #include <gp.hxx>
  31 #include <gp_Trsf.hxx>
  32 #include <Precision.hxx>
  33 #include <Approx_FitAndDivide2d.hxx>
  34 #include <AppParCurves_MultiCurve.hxx>
  35 #include <Adaptor3d_HCurve.hxx>
  36 #include <Adaptor3d_HSurface.hxx>
  37 #include <TColgp_Array1OfPnt2d.hxx>
  38 #include <TColgp_Array1OfPnt.hxx>
  39 #include <TColStd_Array1OfReal.hxx>
  40 #include <TColStd_Array1OfInteger.hxx>
  41 #include <Geom_BSplineCurve.hxx>
  42 #include <Geom_BezierCurve.hxx>
  43 #include <Geom2d_BSplineCurve.hxx>
  44 #include <Geom2d_BezierCurve.hxx>
  45
  46 //#define DRAW
  47 #ifdef DRAW
  48 #include <DrawTrSurf.hxx>
  49 #endif
  50 #ifdef OCCT_DEBUG
  51 //static Standard_Boolean AffichValue = Standard_False;
  52 #endif
  53
  54 //=======================================================================
  55 //function : IsEqual
  56 //purpose  :
  57 //=======================================================================
  58 // OFV:
  59 static inline Standard_Boolean IsEqual(Standard_Real Check,Standard_Real With,Standard_Real Toler)
  60 {
  61   return ((Abs(Check - With) < Toler) ? Standard_True : Standard_False);
  62 }
  63
  64
  65 //=======================================================================
  66 //function : Value
  67 //purpose  :
  68 //=======================================================================
  69
  70 static gp_Pnt2d Function_Value(const Standard_Real U,
  71                                const Handle(Adaptor3d_HCurve)&   myCurve,
  72                                const Handle(Adaptor3d_HSurface)& mySurface,
  73                                const Standard_Real U1,
  74                                const Standard_Real U2,
  75                                const Standard_Real V1,
  76                                const Standard_Real V2,
  77                                const Standard_Boolean UCouture,
  78                                const Standard_Boolean VCouture )
  79 {
  80   Standard_Real S = 0., T = 0.;
  81
  82   gp_Pnt P3d = myCurve->Value(U);
  83   GeomAbs_SurfaceType SType = mySurface->GetType();
  84
  85   switch ( SType ) {
  86
  87   case GeomAbs_Plane:
  88     {
  89       gp_Pln Plane = mySurface->Plane();
  90       ElSLib::Parameters( Plane, P3d, S, T);
  91       break;
  92     }
  93   case GeomAbs_Cylinder:
  94     {
  95       gp_Cylinder Cylinder = mySurface->Cylinder();
  96       ElSLib::Parameters( Cylinder, P3d, S, T);
  97       break;
  98     }
  99   case GeomAbs_Cone:
 100     {
 101       gp_Cone Cone = mySurface->Cone();
 102       ElSLib::Parameters( Cone, P3d, S, T);
 103       break;
 104     }
 105   case GeomAbs_Sphere:
 106     {
 107       gp_Sphere Sphere = mySurface->Sphere();
 108       ElSLib::Parameters(Sphere, P3d, S, T);
 109       break;
 110     }
 111   case GeomAbs_Torus:
 112     {
 113       gp_Torus Torus = mySurface->Torus();
 114       ElSLib::Parameters( Torus, P3d, S, T);
 115       break;
 116     }
 117   default:
 118     Standard_NoSuchObject::Raise("ProjLib_ComputeApprox::Value");
 119   }
 120
 121   if ( UCouture) {
 122     if(S < U1 || S > U2)
 123     {
 124         S = ElCLib::InPeriod(S, U1, U2);
 125     }
 126   }
 127
 128   if ( VCouture) {
 129     if(SType == GeomAbs_Sphere) {
 130       if ( Abs( S - U1 ) > M_PI ) {
 131               T = M_PI - T;
 132               S = M_PI + S;
 133       }
 134       if(S > U1 || S < U2)
 135         S = ElCLib::InPeriod(S, U1, U2);
 136     }
 137     if(T < V1 || T > V2)
 138       T = ElCLib::InPeriod(T, V1, V2);
 139   }
 140
 141   return gp_Pnt2d(S, T);
 142 }
 143 //=======================================================================
 144 //function : D1
 145 //purpose  :
 146 //=======================================================================
 147 static Standard_Boolean Function_D1( const Standard_Real U,
 148                                     gp_Pnt2d&            P,
 149                                     gp_Vec2d&            D,
 150                                     const Handle(Adaptor3d_HCurve)&   myCurve,
 151                                     const Handle(Adaptor3d_HSurface)& mySurface,
 152                                     const Standard_Real U1,
 153                                     const Standard_Real U2,
 154                                     const Standard_Real V1,
 155                                     const Standard_Real V2,
 156                                     const Standard_Boolean UCouture,
 157                                     const Standard_Boolean VCouture )
 158 {
 159   gp_Pnt P3d;
 160   Standard_Real dU, dV;
 161
 162   P = Function_Value(U,myCurve,mySurface,U1,U2,V1,V2,UCouture,VCouture);
 163
 164   GeomAbs_SurfaceType Type = mySurface->GetType();
 165
 166   switch ( Type) {
 167   case GeomAbs_Plane:
 168   case GeomAbs_Cone:
 169   case GeomAbs_Cylinder:
 170   case GeomAbs_Sphere:
 171   case GeomAbs_Torus:
 172     {
 173       gp_Vec D1U, D1V;
 174       gp_Vec T;
 175       myCurve->D1(U,P3d,T);
 176       mySurface->D1(P.X(),P.Y(),P3d,D1U,D1V);
 177
 178       dU = T.Dot(D1U);
 179       dV = T.Dot(D1V);
 180       Standard_Real Nu = D1U.SquareMagnitude();
 181       Standard_Real Nv = D1V.SquareMagnitude();
 182
 183       if ( Nu < Epsilon(1.) || Nv < Epsilon(1.))
 184               return Standard_False;
 185
 186       dU /= Nu;
 187       dV /= Nv;
 188       D = gp_Vec2d( dU, dV);
 189     }
 190     break;
 191
 192   default:
 193     return Standard_False;
 194   }
 195
 196   return Standard_True;
 197 }
 198
 199 //=======================================================================
 200 //function : Function_SetUVBounds
 201 //purpose  :
 202 //=======================================================================
 203 static void Function_SetUVBounds(Standard_Real& myU1,
 204                                  Standard_Real& myU2,
 205                                  Standard_Real& myV1,
 206                                  Standard_Real& myV2,
 207                                  Standard_Boolean& UCouture,
 208                                  Standard_Boolean& VCouture,
 209                                  const Handle(Adaptor3d_HCurve)&   myCurve,
 210                                  const Handle(Adaptor3d_HSurface)& mySurface)
 211 {
 212   Standard_Real W1, W2, W;
 213   gp_Pnt P1, P2, P;
 214   //
 215   W1 = myCurve->FirstParameter();
 216   W2 = myCurve->LastParameter ();
 217   W  = 0.5*(W1+W2);
 218   // on ouvre l`intervalle
 219   // W1 += 1.0e-9;
 220   // W2 -= 1.0e-9;
 221   P1 = myCurve->Value(W1);
 222   P2 = myCurve->Value(W2);
 223   P  = myCurve->Value(W);
 224
 225   switch ( mySurface->GetType()) {
 226
 227     case GeomAbs_Cone:    {
 228       Standard_Real tol = Epsilon(1.);
 229       Standard_Real ptol = Precision::PConfusion();
 230       gp_Cone Cone = mySurface->Cone();
 231       VCouture = Standard_False;
 232       //Calculation of cone parameters for P == ConeApex often produces wrong
 233       //values of U
 234       gp_Pnt ConeApex = Cone.Apex();
 235       if(ConeApex.XYZ().IsEqual(P1.XYZ(), tol))
 236       {
 237         W1 += ptol;
 238         P1 = myCurve->Value(W1);
 239       }
 240       if(ConeApex.XYZ().IsEqual(P2.XYZ(), tol))
 241       {
 242         W2 -= ptol;
 243         P2 = myCurve->Value(W2);
 244       }
 245       if(ConeApex.XYZ().IsEqual(P.XYZ(), tol))
 246       {
 247         W += ptol;
 248         P = myCurve->Value(W);
 249       }
 250
 251       switch( myCurve->GetType() ){
 252       case GeomAbs_Parabola:
 253       case GeomAbs_Hyperbola:
 254       case GeomAbs_Ellipse:{
 255         Standard_Real U1, U2, V1, V2, U , V;
 256         ElSLib::Parameters( Cone, P1, U1, V1);
 257         ElSLib::Parameters( Cone, P2, U2, V2);
 258         ElSLib::Parameters( Cone, P , U , V );
 259         myU1 = Min(U1,U2);
 260         myU2 = Max(U1,U2);
 261         if (  ( U1 < U  &&  U < U2 ) && !myCurve->IsClosed() ) {
 262           UCouture = Standard_False;
 263         }
 264         else {
 265           UCouture = Standard_True;
 266           myU2 = myU1 + 2*M_PI;
 267         }
 268
 269       }
 270       break;
 271       default:    {
 272         Standard_Real U1, V1, U , V, Delta = 0., d = 0., pmin = W1, pmax = W1, dmax = 0., Uf, Ul;
 273         ElSLib::Parameters( Cone, P1, U1, V1);
 274         ElSLib::Parameters( Cone, P2, Ul, V1);
 275         myU1 = U1; myU2 = U1; Uf = U1;
 276         Standard_Real Step = .1;
 277         Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
 278         if(myCurve->GetType() == GeomAbs_Line)
 279         {
 280           nbp = 3;
 281         }
 282         nbp = Max(nbp, 3);
 283         Step = (W2 - W1) / (nbp - 1);
 284         Standard_Boolean isclandper = (!(myCurve->IsClosed()) && !(myCurve->IsPeriodic()));
 285         Standard_Boolean isFirst = Standard_True;
 286         for(Standard_Real par = W1 + Step; par <= W2; par += Step)
 287         {
 288           if(!isclandper) par += Step;
 289           P = myCurve->Value(par);
 290           ElSLib::Parameters( Cone, P, U, V);
 291           U += Delta;
 292           d = U - U1;
 293           if(d > M_PI)
 294           {
 295             if( ( (IsEqual(U,(2*M_PI),1.e-10) && (U1 >= 0. && U1 <= M_PI)) &&
 296               (IsEqual(U,Ul,1.e-10) && !IsEqual(Uf,0.,1.e-10)) ) && isclandper )
 297               U = 0.0;
 298             else
 299             {
 300               // Protection against first-last point on seam.
 301               if (isFirst)
 302                 U1 = 2*M_PI;
 303               else if (par + Step >= W2)
 304                 U = 0.0;
 305               else
 306                 Delta -= 2*M_PI;
 307             }
 308             U += Delta;
 309             d = U - U1;
 310           }
 311           else if(d < -M_PI)
 312           {
 313             if( ( (IsEqual(U,0.,1.e-10) && (U1 >= M_PI && U1 <= (2*M_PI))) &&
 314               (IsEqual(U,Ul,1.e-10) && !IsEqual(Uf,(2*M_PI),1.e-10)) ) && isclandper )
 315               U = 2*M_PI;
 316             else
 317             {
 318               // Protection against first-last point on seam.
 319               if (isFirst)
 320                 U1 = 0.0;
 321               else if (par + Step >= W2)
 322                 U = 2*M_PI;
 323               else
 324                 Delta += 2*M_PI;
 325             }
 326             U += Delta;
 327             d = U - U1;
 328           }
 329           dmax = Max(dmax, Abs(d));
 330           if(U < myU1) {myU1 = U; pmin = par;}
 331           if(U > myU2) {myU2 = U; pmax = par;}
 332           U1 = U;
 333           isFirst = Standard_False;
 334         } // for(Standard_Real par = W1 + Step; par <= W2; par += Step)
 335
 336         if(!(Abs(pmin - W1) <= Precision::PConfusion() || Abs(pmin - W2) <= Precision::PConfusion()) ) myU1 -= dmax*.5;
 337         if(!(Abs(pmax - W1) <= Precision::PConfusion() || Abs(pmax - W2) <= Precision::PConfusion()) ) myU2 += dmax*.5;
 338
 339         if((myU1 >=0. && myU1 <= 2*M_PI) && (myU2 >=0. && myU2 <= 2*M_PI) ) UCouture = Standard_False;
 340         else{
 341           U = ( myU1 + myU2 ) /2.;
 342           myU1 = U - M_PI;
 343           myU2 = U + M_PI;
 344           UCouture = Standard_True;
 345         }
 346       }
 347       break;
 348     }// switch curve type
 349   }// case Cone
 350   break;
 351
 352   case GeomAbs_Cylinder:    {
 353     gp_Cylinder Cylinder = mySurface->Cylinder();
 354     VCouture = Standard_False;
 355
 356     if (myCurve->GetType() == GeomAbs_Ellipse) {
 357
 358       Standard_Real U1, U2, V1, V2, U , V;
 359       ElSLib::Parameters( Cylinder, P1, U1, V1);
 360       ElSLib::Parameters( Cylinder, P2, U2, V2);
 361       ElSLib::Parameters( Cylinder, P , U , V );
 362       myU1 = Min(U1,U2);
 363       myU2 = Max(U1,U2);
 364
 365       if ( !myCurve->IsClosed()) {
 366         if ( myU1 < U && U < myU2) {
 367           U = ( myU1 + myU2 ) /2.;
 368           myU1 = U - M_PI;
 369           myU2 = U + M_PI;
 370         }
 371         else {
 372           U = ( myU1 + myU2 ) /2.;
 373           if ( myU1 < U) {
 374             myU1 = U - 2*M_PI;
 375             myU2 = U;
 376           }
 377           else {
 378             myU1 = U;
 379             myU2 = U + 2*M_PI;
 380           }
 381         }
 382         UCouture = Standard_True;
 383       }
 384       else {
 385         gp_Vec D1U, D1V;
 386         gp_Vec T;
 387         gp_Pnt P3d;
 388         myCurve->D1(W1,P3d,T);
 389         mySurface->D1(U1,U2,P3d,D1U,D1V);
 390         Standard_Real dU = T.Dot(D1U);
 391
 392         UCouture = Standard_True;
 393         if ( dU > 0.) {
 394           myU2 = myU1 + 2*M_PI;
 395         }
 396         else {
 397           myU2 = myU1;
 398           myU1 -= 2*M_PI;
 399         }
 400       }
 401     }
 402     else {
 403       Standard_Real U1, V1, U , V;
 404       ElSLib::Parameters( Cylinder, P1, U1, V1);
 405       Standard_Real Step = .1, Delta = 0.;
 406       Standard_Real eps = M_PI, dmax = 0., d = 0.;
 407       Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
 408       nbp = Max(nbp, 3);
 409       Step = (W2 - W1) / (nbp - 1);
 410       myU1 = U1; myU2 = U1;
 411       Standard_Real pmin = W1, pmax = W1, plim = W2+.1*Step;
 412       for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
 413         P = myCurve->Value(par);
 414         ElSLib::Parameters( Cylinder, P, U, V);
 415         U += Delta;
 416         d = U - U1;
 417         if(d > eps) {
 418           U -= Delta;
 419           Delta -= 2*M_PI;
 420           U += Delta;
 421           d = U - U1;
 422         }
 423         else if(d < -eps) {
 424           U -= Delta;
 425           Delta += 2*M_PI;
 426           U += Delta;
 427           d = U - U1;
 428         }
 429         dmax = Max(dmax, Abs(d));
 430         if(U < myU1) {myU1 = U; pmin = par;}
 431         if(U > myU2) {myU2 = U; pmax = par;}
 432         U1 = U;
 433       }
 434
 435       if(!(Abs(pmin - W1) <= Precision::PConfusion() ||
 436         Abs(pmin - W2) <= Precision::PConfusion())  ) myU1 -= dmax*.5;
 437       if(!(Abs(pmax - W1) <= Precision::PConfusion() ||
 438         Abs(pmax - W2) <= Precision::PConfusion())  ) myU2 += dmax*.5;
 439
 440       if((myU1 >=0. && myU1 <= 2*M_PI) &&
 441         (myU2 >=0. && myU2 <= 2*M_PI)    ) {
 442           UCouture = Standard_False;
 443       }
 444       else {
 445         U = ( myU1 + myU2 ) /2.;
 446         myU1 = U - M_PI;
 447         myU2 = U + M_PI;
 448         UCouture = Standard_True;
 449       }
 450     }
 451   }
 452   break;
 453   //
 454   case GeomAbs_Sphere:{
 455     VCouture = Standard_False;
 456     gp_Sphere SP = mySurface->Sphere();
 457     if ( myCurve->GetType() == GeomAbs_Circle) {
 458       UCouture = Standard_True;
 459
 460       // on cherche a savoir le nombre de fois que la couture est
 461       // traversee.
 462       // si 0 ou 2 fois : la PCurve est fermee et dans l`intervalle
 463       //                  [Uc-PI, Uc+PI] (Uc: U du centre du cercle)
 464       // si 1 fois      : la PCurve est ouverte et dans l`intervalle
 465       //                  [U1, U1 +/- 2*PI]
 466
 467       // pour determiner le nombre de solution, on resoud le systeme
 468       // x^2 + y^2 + z^2     = R^2  (1)
 469       // A x + B y + C z + D = 0    (2)
 470       // x > 0                      (3)
 471       // y = 0                      (4)
 472       // REM : (1) (2)     : equation du cercle
 473       //       (1) (3) (4) : equation de la couture.
 474       Standard_Integer NbSolutions = 0;
 475       Standard_Real A, B, C, D, R, Tol = 1.e-10;
 476       Standard_Real U1, U2, V1, V2;
 477       gp_Trsf Trsf;
 478       //
 479       gp_Circ Circle = myCurve->Circle();
 480       Trsf.SetTransformation(SP.Position());
 481       Circle.Transform(Trsf);
 482       //
 483       R = SP.Radius();
 484       gp_Pln Plane( gp_Ax3(Circle.Position()));
 485       Plane.Coefficients(A,B,C,D);
 486       //
 487       if ( Abs(C) < Tol) {
 488         if ( Abs(A) > Tol) {
 489           if ( (D/A) < 0.) {
 490             if      ( ( R - Abs(D/A))  > Tol)  NbSolutions = 2;
 491             else if ( Abs(R - Abs(D/A))< Tol)  NbSolutions = 1;
 492             else                               NbSolutions = 0;
 493           }
 494         }
 495       }
 496       else {
 497         Standard_Real delta = R*R*(A*A+C*C) - D*D;
 498         delta *= C*C;
 499         if ( Abs(delta) < Tol*Tol) {
 500           if ( A*D > 0.) NbSolutions = 1;
 501         }
 502         else if ( delta > 0) {
 503           Standard_Real xx;
 504           delta = Sqrt(delta);
 505           xx = -A*D+delta;
 506           //
 507           if ( xx > Tol) NbSolutions++;
 508           xx = -A*D-delta;
 509           //
 510           if ( xx > Tol) NbSolutions++;
 511         }
 512       }
 513       //
 514
 515       // box+sphere >>
 516       Standard_Real UU = 0.;
 517       ElSLib::Parameters(SP, P1, U1, V1);
 518       Standard_Real eps = 10.*Epsilon(1.);
 519       Standard_Real dt = Max(Precision::PConfusion(), 0.01*(W2-W1));
 520       if(Abs(U1) < eps)
 521       {
 522         //May be U1 must be equal 2*PI?
 523         gp_Pnt Pd  = myCurve->Value(W1+dt);
 524         Standard_Real ud, vd;
 525         ElSLib::Parameters(SP, Pd, ud, vd);
 526         if(Abs(U1 - ud) > M_PI)
 527         {
 528           U1 = 2.*M_PI;
 529         }
 530       }
 531       else if(Abs(2.*M_PI - U1) < eps)
 532       {
 533         //maybe U1 = 0.?
 534         gp_Pnt Pd  = myCurve->Value(W1+dt);
 535         Standard_Real ud, vd;
 536         ElSLib::Parameters(SP, Pd, ud, vd);
 537         if(Abs(U1 - ud) > M_PI)
 538         {
 539           U1 = 0.;
 540         }
 541       }
 542       //
 543       ElSLib::Parameters(SP, P2, U2, V1);
 544       if(Abs(U2) < eps)
 545       {
 546         //May be U2 must be equal 2*PI?
 547         gp_Pnt Pd  = myCurve->Value(W2-dt);
 548         Standard_Real ud, vd;
 549         ElSLib::Parameters(SP, Pd, ud, vd);
 550         if(Abs(U2 - ud) > M_PI)
 551         {
 552           U2 = 2.*M_PI;
 553         }
 554       }
 555       else if(Abs(2.*M_PI - U2) < eps)
 556       {
 557         //maybe U2 = 0.?
 558         gp_Pnt Pd  = myCurve->Value(W2-dt);
 559         Standard_Real ud, vd;
 560         ElSLib::Parameters(SP, Pd, ud, vd);
 561         if(Abs(U2 - ud) > M_PI)
 562         {
 563           U2 = 0.;
 564         }
 565       }
 566       //
 567       ElSLib::Parameters(SP, P, UU, V1);
 568       //+This fragment was the reason of bug # 26008.
 569       //+It has been deleted on April, 03 2015.
 570       //Standard_Real UUmi = Min(Min(U1,UU),Min(UU,U2));
 571       //Standard_Real UUma = Max(Max(U1,UU),Max(UU,U2));
 572       //Standard_Boolean reCalc = ((UUmi >= 0. && UUmi <= M_PI) && (UUma >= 0. && UUma <= M_PI));
 573       // box+sphere <<
 574       P2 = myCurve->Value(W1+M_PI/8);
 575       ElSLib::Parameters(SP,P2,U2,V2);
 576       //
 577       if ( NbSolutions == 1) {
 578         if ( Abs(U1-U2) > M_PI) { // on traverse la couture
 579           if ( U1 > M_PI) {
 580             myU1 = U1;
 581             myU2 = U1+2*M_PI;
 582           }
 583           else {
 584             myU2 = U1;
 585             myU1 = U1-2*M_PI;
 586           }
 587         }
 588         else { // on ne traverse pas la couture
 589           if ( U1 > U2) {
 590             myU2 = U1;
 591             myU1 = U1-2*M_PI;
 592           }
 593           else {
 594             myU1 = U1;
 595             myU2 = U1+2*M_PI;
 596           }
 597         }
 598       }
 599       else { // 0 ou 2 solutions
 600         gp_Pnt Center = Circle.Location();
 601         Standard_Real U,V;
 602         ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
 603         myU1 = U-M_PI;
 604         myU2 = U+M_PI;
 605       }
 606       //
 607       // eval the VCouture.
 608       if ( (C==0) || Abs(Abs(D/C)-R) > 1.e-10) {
 609         VCouture = Standard_False;
 610       }
 611       else {
 612         VCouture = Standard_True;
 613         UCouture = Standard_True;
 614
 615         if ( D/C < 0.) {
 616           myV1 =   - M_PI / 2.;
 617           myV2 = 3 * M_PI / 2.;
 618         }
 619         else {
 620           myV1 = -3 * M_PI / 2.;
 621           myV2 =      M_PI / 2.;
 622         }
 623
 624         // si P1.Z() vaut +/- R on est sur le sommet : pas significatif.
 625         gp_Pnt pp = P1.Transformed(Trsf);
 626
 627         if ( Abs( Abs(pp.Z()) - R) < Tol) {
 628           gp_Pnt Center = Circle.Location();
 629           Standard_Real U,V;
 630           ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
 631           myU1 = U-M_PI;
 632           myU2 = U+M_PI;
 633           VCouture = Standard_False;
 634         }
 635       }
 636
 637       // box+sphere >>
 638       myV1 = -1.e+100; myV2 = 1.e+100;
 639
 640       //+This fragment was the reason of bug # 26008.
 641       //+It has been deleted on April, 03 2015.
 642       //Standard_Real UU1 = myU1, UU2 = myU2;
 643       //if((Abs(UU1) <= (2.*M_PI) && Abs(UU2) <= (2.*M_PI)) && NbSolutions == 1 && reCalc) {
 644       //  gp_Pnt Center = Circle.Location();
 645       //  Standard_Real U,V;
 646       //  ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
 647       //  myU1 = U-M_PI;
 648       //  myU1 = Min(UU1,myU1);
 649       //  myU2 = myU1 + 2.*M_PI;
 650       //}
 651       // box+sphere <<
 652
 653     }//if ( myCurve->GetType() == GeomAbs_Circle)
 654
 655     else {
 656       Standard_Real U1, V1, U , V;
 657       ElSLib::Parameters( SP, P1, U1, V1);
 658       Standard_Real Step = .1, Delta = 0.;
 659       Standard_Real eps = M_PI, dmax = 0., d = 0.;
 660       Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
 661       nbp = Max(nbp, 3);
 662       Step = (W2 - W1) / (nbp - 1);
 663       myU1 = U1; myU2 = U1;
 664       Standard_Real pmin = W1, pmax = W1, plim = W2+.1*Step;
 665       for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
 666         P = myCurve->Value(par);
 667         ElSLib::Parameters( SP, P, U, V);
 668         U += Delta;
 669         d = U - U1;
 670         if(d > eps) {
 671           U -= Delta;
 672           Delta -= 2*M_PI;
 673           U += Delta;
 674           d = U - U1;
 675         }
 676         else if(d < -eps) {
 677           U -= Delta;
 678           Delta += 2*M_PI;
 679           U += Delta;
 680           d = U - U1;
 681         }
 682         dmax = Max(dmax, Abs(d));
 683         if(U < myU1) {myU1 = U; pmin = par;}
 684         if(U > myU2) {myU2 = U; pmax = par;}
 685         U1 = U;
 686       }
 687
 688       if(!(Abs(pmin - W1) <= Precision::PConfusion() ||
 689         Abs(pmin - W2) <= Precision::PConfusion())  ) myU1 -= dmax*.5;
 690       if(!(Abs(pmax - W1) <= Precision::PConfusion() ||
 691         Abs(pmax - W2) <= Precision::PConfusion())  ) myU2 += dmax*.5;
 692
 693       if((myU1 >=0. && myU1 <= 2*M_PI) &&
 694         (myU2 >=0. && myU2 <= 2*M_PI)    ) {
 695           myU1 = 0.;
 696           myU2 = 2.*M_PI;
 697           UCouture = Standard_False;
 698       }
 699       else {
 700         U = ( myU1 + myU2 ) /2.;
 701         myU1 = U - M_PI;
 702         myU2 = U + M_PI;
 703         UCouture = Standard_True;
 704       }
 705
 706       VCouture = Standard_False;
 707     }
 708   }
 709   break;
 710   //
 711   case GeomAbs_Torus:{
 712     gp_Torus TR = mySurface->Torus();
 713     Standard_Real U1, V1, U , V;
 714     ElSLib::Parameters( TR, P1, U1, V1);
 715     Standard_Real Step = .1, DeltaU = 0., DeltaV = 0.;
 716     Standard_Real eps = M_PI, dmaxU = 0., dU = 0., dmaxV = 0., dV = 0.;
 717     Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
 718     nbp = Max(nbp, 3);
 719     Step = (W2 - W1) / (nbp - 1);
 720     myU1 = U1; myU2 = U1;
 721     myV1 = V1; myV2 = V1;
 722     Standard_Real pminU = W1, pmaxU = W1, pminV = W1, pmaxV = W1,
 723       plim = W2+.1*Step;
 724     for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
 725       P = myCurve->Value(par);
 726       ElSLib::Parameters( TR, P, U, V);
 727       U += DeltaU;
 728       V += DeltaV;
 729       dU = U - U1;
 730       dV = V - V1;
 731       if(dU > eps) {
 732         U -= DeltaU;
 733         DeltaU -= 2*M_PI;
 734         U += DeltaU;
 735         dU = U - U1;
 736       }
 737       else if(dU < -eps) {
 738         U -= DeltaU;
 739         DeltaU += 2*M_PI;
 740         U += DeltaU;
 741         dU = U - U1;
 742       }
 743       if(dV > eps) {
 744         V -= DeltaV;
 745         DeltaV -= 2*M_PI;
 746         V += DeltaV;
 747         dV = V - V1;
 748       }
 749       else if(dV < -eps) {
 750         V -= DeltaV;
 751         DeltaV += 2*M_PI;
 752         V += DeltaV;
 753         dV = V - V1;
 754       }
 755       dmaxU = Max(dmaxU, Abs(dU));
 756       dmaxV = Max(dmaxV, Abs(dV));
 757       if(U < myU1) {myU1 = U; pminU = par;}
 758       if(U > myU2) {myU2 = U; pmaxU = par;}
 759       if(V < myV1) {myV1 = V; pminV = par;}
 760       if(V > myV2) {myV2 = V; pmaxV = par;}
 761       U1 = U;
 762       V1 = V;
 763     }
 764
 765     if(!(Abs(pminU - W1) <= Precision::PConfusion() ||
 766       Abs(pminU - W2) <= Precision::PConfusion())  ) myU1 -= dmaxU*.5;
 767     if(!(Abs(pmaxU - W1) <= Precision::PConfusion() ||
 768       Abs(pmaxU - W2) <= Precision::PConfusion())  ) myU2 += dmaxU*.5;
 769     if(!(Abs(pminV - W1) <= Precision::PConfusion() ||
 770       Abs(pminV - W2) <= Precision::PConfusion())  ) myV1 -= dmaxV*.5;
 771     if(!(Abs(pmaxV - W1) <= Precision::PConfusion() ||
 772       Abs(pmaxV - W2) <= Precision::PConfusion())  ) myV2 += dmaxV*.5;
 773
 774     if((myU1 >=0. && myU1 <= 2*M_PI) &&
 775       (myU2 >=0. && myU2 <= 2*M_PI)    ) {
 776         myU1 = 0.;
 777         myU2 = 2.*M_PI;
 778         UCouture = Standard_False;
 779     }
 780     else {
 781       U = ( myU1 + myU2 ) /2.;
 782       myU1 = U - M_PI;
 783       myU2 = U + M_PI;
 784       UCouture = Standard_True;
 785     }
 786     if((myV1 >=0. && myV1 <= 2*M_PI) &&
 787       (myV2 >=0. && myV2 <= 2*M_PI)    ) {
 788         VCouture = Standard_False;
 789     }
 790     else {
 791       V = ( myV1 + myV2 ) /2.;
 792       myV1 = V - M_PI;
 793       myV2 = V + M_PI;
 794       VCouture = Standard_True;
 795     }
 796
 797   }
 798   break;
 799
 800   default:
 801     {
 802       UCouture = Standard_False;
 803       VCouture = Standard_False;
 804     }
 805     break;
 806   }
 807 }
 808 //
 809 //
 810 //=======================================================================
 811 //classn : ProjLib_Function
 812 //purpose  :
 813 //=======================================================================
 814 class ProjLib_Function : public AppCont_Function
 815 {
 816   Handle(Adaptor3d_HCurve)   myCurve;
 817   Handle(Adaptor3d_HSurface) mySurface;
 818   Standard_Boolean myIsPeriodic[2];
 819   Standard_Real myPeriod[2];
 820   public :
 821
 822   Standard_Real    myU1,myU2,myV1,myV2;
 823   Standard_Boolean UCouture,VCouture;
 824
 825   ProjLib_Function(const Handle(Adaptor3d_HCurve)&   C,
 826                    const Handle(Adaptor3d_HSurface)& S)
 827 : myCurve(C),
 828   mySurface(S),
 829   myU1(0.0),
 830   myU2(0.0),
 831   myV1(0.0),
 832   myV2(0.0),
 833   UCouture(Standard_False),
 834   VCouture(Standard_False)
 835   {
 836     myNbPnt = 0;
 837     myNbPnt2d = 1;
 838     Function_SetUVBounds(myU1,myU2,myV1,myV2,UCouture,VCouture,myCurve,mySurface);
 839     myIsPeriodic[0] = mySurface->IsUPeriodic();
 840     myIsPeriodic[1] = mySurface->IsVPeriodic();
 841
 842     if (myIsPeriodic[0])
 843       myPeriod[0] = mySurface->UPeriod();
 844     else
 845       myPeriod[0] = 0.0;
 846
 847     if (myIsPeriodic[1])
 848       myPeriod[1] = mySurface->VPeriod();
 849     else
 850       myPeriod[1] = 0.0;
 851   }
 852
 853   void PeriodInformation(const Standard_Integer theDimIdx,
 854                          Standard_Boolean& IsPeriodic,
 855                          Standard_Real& thePeriod) const
 856   {
 857     IsPeriodic = myIsPeriodic[theDimIdx - 1];
 858     thePeriod = myPeriod[theDimIdx - 1];
 859   }
 860
 861   Standard_Real FirstParameter() const
 862   {
 863     return (myCurve->FirstParameter());
 864   }
 865
 866   Standard_Real LastParameter() const
 867   {
 868     return (myCurve->LastParameter());
 869   }
 870
 871   Standard_Boolean Value(const Standard_Real   theT,
 872                          NCollection_Array1<gp_Pnt2d>& thePnt2d,
 873                          NCollection_Array1<gp_Pnt>&   /*thePnt*/) const
 874   {
 875     thePnt2d(1) = Function_Value(theT, myCurve, mySurface, myU1, myU2, myV1, myV2, UCouture, VCouture);
 876     return Standard_True;
 877   }
 878
 879   gp_Pnt2d Value(const Standard_Real   theT) const
 880   {
 881     return Function_Value(theT, myCurve, mySurface, myU1, myU2, myV1, myV2, UCouture, VCouture);
 882   }
 883
 884   Standard_Boolean D1(const Standard_Real   theT,
 885                       NCollection_Array1<gp_Vec2d>& theVec2d,
 886                       NCollection_Array1<gp_Vec>&   /*theVec*/) const
 887   {
 888     gp_Pnt2d aPnt2d;
 889     gp_Vec2d aVec2d;
 890     Standard_Boolean isOk = Function_D1(theT, aPnt2d,aVec2d, myCurve, mySurface, myU1, myU2, myV1, myV2, UCouture, VCouture);
 891     theVec2d(1) = aVec2d;
 892     return isOk;
 893   }
 894 };
 895
 896 //=======================================================================
 897 //function : ProjLib_ComputeApprox
 898 //purpose  :
 899 //=======================================================================
 900
 901 ProjLib_ComputeApprox::ProjLib_ComputeApprox
 902   (const Handle(Adaptor3d_HCurve)   & C,
 903    const Handle(Adaptor3d_HSurface) & S,
 904    const Standard_Real              Tol )
 905 {
 906   // if the surface is a plane and the curve a BSpline or a BezierCurve,
 907   // don`t make an Approx but only the projection of the poles.
 908
 909   myTolerance = Max(Precision::PApproximation(),Tol);
 910   Standard_Integer NbKnots, NbPoles ;
 911   GeomAbs_CurveType   CType = C->GetType();
 912   GeomAbs_SurfaceType SType = S->GetType();
 913
 914   Standard_Boolean SurfIsAnal = (SType != GeomAbs_BSplineSurface) &&
 915                                 (SType != GeomAbs_BezierSurface)  &&
 916                                 (SType != GeomAbs_OtherSurface)     ;
 917
 918   Standard_Boolean CurvIsAnal = (CType != GeomAbs_BSplineCurve) &&
 919                                 (CType != GeomAbs_BezierCurve)  &&
 920                                 (CType != GeomAbs_OffsetCurve)  &&
 921                                 (CType != GeomAbs_OtherCurve)     ;
 922
 923   Standard_Boolean simplecase = SurfIsAnal && CurvIsAnal;
 924
 925   if (CType == GeomAbs_BSplineCurve &&
 926       SType == GeomAbs_Plane ) {
 927
 928     // get the poles and eventually the weights
 929     Handle(Geom_BSplineCurve) BS = C->BSpline();
 930     NbPoles = BS->NbPoles();
 931     TColgp_Array1OfPnt P3d( 1, NbPoles);
 932     TColgp_Array1OfPnt2d Poles( 1, NbPoles);
 933     TColStd_Array1OfReal Weights( 1, NbPoles);
 934     if ( BS->IsRational()) BS->Weights(Weights);
 935     BS->Poles( P3d);
 936     gp_Pln Plane = S->Plane();
 937     Standard_Real U,V;
 938     for ( Standard_Integer i = 1; i <= NbPoles; i++) {
 939       ElSLib::Parameters( Plane, P3d(i), U, V);
 940       Poles.SetValue(i,gp_Pnt2d(U,V));
 941     }
 942     NbKnots = BS->NbKnots();
 943     TColStd_Array1OfReal     Knots(1,NbKnots);
 944     TColStd_Array1OfInteger  Mults(1,NbKnots);
 945     BS->Knots(Knots) ;
 946     BS->Multiplicities(Mults) ;
 947     // get the knots and mults if BSplineCurve
 948     if ( BS->IsRational()) {
 949       myBSpline = new Geom2d_BSplineCurve(Poles,
 950                                           Weights,
 951                                           Knots,
 952                                           Mults,
 953                                           BS->Degree(),
 954                                           BS->IsPeriodic());
 955     }
 956     else {
 957       myBSpline = new Geom2d_BSplineCurve(Poles,
 958                                           Knots,
 959                                           Mults,
 960                                           BS->Degree(),
 961                                           BS->IsPeriodic());
 962     }
 963   }
 964   else if (CType == GeomAbs_BezierCurve &&
 965            SType == GeomAbs_Plane ) {
 966
 967     // get the poles and eventually the weights
 968     Handle(Geom_BezierCurve) BezierCurvePtr = C->Bezier() ;
 969     NbPoles = BezierCurvePtr->NbPoles();
 970     TColgp_Array1OfPnt P3d( 1, NbPoles);
 971     TColgp_Array1OfPnt2d  Poles( 1, NbPoles);
 972     TColStd_Array1OfReal   Weights( 1, NbPoles);
 973     if ( BezierCurvePtr->IsRational()) {
 974       BezierCurvePtr->Weights(Weights);
 975     }
 976     BezierCurvePtr->Poles( P3d);
 977
 978     // project the 3D-Poles on the plane
 979
 980     gp_Pln Plane = S->Plane();
 981     Standard_Real U,V;
 982     for ( Standard_Integer i = 1; i <= NbPoles; i++) {
 983       ElSLib::Parameters( Plane, P3d(i), U, V);
 984       Poles.SetValue(i,gp_Pnt2d(U,V));
 985     }
 986     if (  BezierCurvePtr->IsRational()) {
 987       myBezier =  new Geom2d_BezierCurve(Poles, Weights);
 988     }
 989     else {
 990       myBezier =  new Geom2d_BezierCurve(Poles);
 991     }
 992   }
 993   else {
 994     ProjLib_Function F( C, S);
 995
 996 #ifdef OCCT_DEBUG
 997     //if ( AffichValue) {
 998     //  Standard_Integer Nb = 20;
 999     //  Standard_Real U1, U2, dU, U;
1000     //  U1 = F.FirstParameter();
1001     //  U2 = F.LastParameter();
1002     //  dU = ( U2 - U1) / Nb;
1003     //  TColStd_Array1OfInteger Mults(1,Nb+1);
1004     //  TColStd_Array1OfReal    Knots(1,Nb+1);
1005     //  TColgp_Array1OfPnt2d    Poles(1,Nb+1);
1006     //  for ( Standard_Integer i = 1; i <= Nb+1; i++) {
1007     //   U = U1 + (i-1)*dU;
1008     //   Poles(i) = F.Value(U);
1009     //   cout << "i = " << i << ": U = " << U <<
1010     //      ", p(" << Poles(i).X() << ", " << Poles(i).Y() << ");" << endl;
1011     //   Knots(i) = i;
1012     //   Mults(i) = 1;
1013     //  }
1014     //  Mults(1)    = 2;
1015     //  Mults(Nb+1) = 2;
1016
1017     //2D-curve for showing in DRAW
1018     //  Handle(Geom2d_Curve) aCC = new Geom2d_BSplineCurve(Poles,Knots,Mults,1);
1019     //  AffichValue = Standard_False;
1020     //}
1021 #endif
1022
1023 //-----------
1024     Standard_Integer Deg1, Deg2;
1025     if(simplecase) {
1026       Deg1 = 8;
1027       Deg2 = 10;
1028     }
1029     else {
1030       Deg1 = 8;
1031       Deg2 = 12;
1032     }
1033 //-------------
1034     Approx_FitAndDivide2d Fit(F,Deg1,Deg2,myTolerance,myTolerance,
1035                               Standard_True);
1036     if(Fit.IsAllApproximated()) {
1037       Standard_Integer i;
1038       Standard_Integer NbCurves = Fit.NbMultiCurves();
1039
1040     // on essaie de rendre la courbe au moins C1
1041       Convert_CompBezierCurves2dToBSplineCurve2d Conv;
1042
1043       myTolerance = 0;
1044       Standard_Real Tol3d,Tol2d;
1045       for (i = 1; i <= NbCurves; i++) {
1046               Fit.Error(i,Tol3d, Tol2d);
1047               myTolerance = Max(myTolerance, Tol2d);
1048               AppParCurves_MultiCurve MC = Fit.Value( i);       //Charge la Ieme Curve
1049               TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Recupere les poles
1050               MC.Curve(1, Poles2d);
1051
1052               Conv.AddCurve(Poles2d);
1053       }
1054
1055       //mise a jour des fields de ProjLib_Approx
1056       Conv.Perform();
1057       NbPoles    = Conv.NbPoles();
1058       NbKnots    = Conv.NbKnots();
1059
1060       if(NbPoles <= 0 || NbPoles > 100000)
1061         return;
1062       if(NbKnots <= 0 || NbKnots > 100000)
1063         return;
1064
1065       TColgp_Array1OfPnt2d    NewPoles(1,NbPoles);
1066       TColStd_Array1OfReal    NewKnots(1,NbKnots);
1067       TColStd_Array1OfInteger NewMults(1,NbKnots);
1068
1069       Conv.KnotsAndMults(NewKnots,NewMults);
1070       Conv.Poles(NewPoles);
1071
1072       BSplCLib::Reparametrize(C->FirstParameter(),
1073                               C->LastParameter(),
1074                               NewKnots);
1075
1076       /*cout << endl;
1077       for (int i = 1; i <= NbPoles; i++)
1078       {
1079         cout << NewPoles.Value(i).X() << " " << NewPoles.Value(i).Y() << endl;
1080       }
1081       cout << endl; */
1082
1083       // il faut recadrer les poles de debut et de fin:
1084       // ( Car pour les problemes de couture, on a du ouvrir l`intervalle
1085       // de definition de la courbe.)
1086       // On choisit de calculer ces poles par prolongement de la courbe
1087       // approximee.
1088       myBSpline = new Geom2d_BSplineCurve (NewPoles,
1089                                            NewKnots,
1090                                            NewMults,
1091                                            Conv.Degree());
1092     }
1093     else {
1094       Standard_Integer NbCurves = Fit.NbMultiCurves();
1095       if(NbCurves != 0) {
1096               Standard_Real Tol3d,Tol2d;
1097               Fit.Error(NbCurves,Tol3d, Tol2d);
1098               myTolerance = Tol2d;
1099       }
1100     }
1101
1102     //Return curve home
1103     Standard_Real UFirst = F.FirstParameter();
1104     gp_Pnt P3d = C->Value( UFirst );
1105     Standard_Real u = 0., v = 0.;
1106     switch (SType)
1107     {
1108     case GeomAbs_Plane:
1109       {
1110         gp_Pln Plane = S->Plane();
1111         ElSLib::Parameters( Plane, P3d, u, v );
1112         break;
1113       }
1114     case GeomAbs_Cylinder:
1115       {
1116         gp_Cylinder Cylinder = S->Cylinder();
1117         ElSLib::Parameters( Cylinder, P3d, u, v );
1118         break;
1119       }
1120     case GeomAbs_Cone:
1121       {
1122         gp_Cone Cone = S->Cone();
1123         ElSLib::Parameters( Cone, P3d, u, v );
1124         break;
1125       }
1126     case GeomAbs_Sphere:
1127       {
1128         gp_Sphere Sphere = S->Sphere();
1129         ElSLib::Parameters( Sphere, P3d, u, v );
1130         break;
1131       }
1132     case GeomAbs_Torus:
1133       {
1134         gp_Torus Torus = S->Torus();
1135         ElSLib::Parameters( Torus, P3d, u, v );
1136         break;
1137       }
1138     default:
1139       Standard_NoSuchObject::Raise("ProjLib_ComputeApprox::Value");
1140     }
1141     Standard_Boolean ToMirror = Standard_False;
1142     Standard_Real du = 0., dv = 0.;
1143     Standard_Integer number;
1144     if (F.VCouture)
1145     {
1146       if (SType == GeomAbs_Sphere && Abs(u-F.myU1) > M_PI)
1147       {
1148         ToMirror = Standard_True;
1149         dv = -M_PI;
1150         v = M_PI - v;
1151       }
1152       Standard_Real newV = ElCLib::InPeriod( v, F.myV1, F.myV2 );
1153       number = (Standard_Integer) (Floor((newV-v)/(F.myV2-F.myV1)));
1154       dv -= number*(F.myV2-F.myV1);
1155     }
1156     if (F.UCouture || (F.VCouture && SType == GeomAbs_Sphere))
1157     {
1158       Standard_Real aNbPer;
1159       gp_Pnt2d P2d = F.Value(UFirst);
1160       du = u - P2d.X();
1161       du = (du < 0) ? (du - Precision::PConfusion()) :
1162         (du + Precision::PConfusion());
1163       modf(du/M_PI, &aNbPer);
1164       number = (Standard_Integer)aNbPer;
1165       du = number*M_PI;
1166     }
1167
1168     if (!myBSpline.IsNull())
1169     {
1170       if (du != 0. || dv != 0.)
1171         myBSpline->Translate( gp_Vec2d(du,dv) );
1172       if (ToMirror)
1173       {
1174         gp_Ax2d Axe( gp_Pnt2d(0.,0.), gp_Dir2d(1.,0.) );
1175         myBSpline->Mirror( Axe );
1176       }
1177     }
1178   }
1179 }
1180
1181 //=======================================================================
1182 //function : BSpline
1183 //purpose  :
1184 //=======================================================================
1185
1186 Handle(Geom2d_BSplineCurve)  ProjLib_ComputeApprox::BSpline() const
1187
1188 {
1189   return myBSpline ;
1190 }
1191
1192 //=======================================================================
1193 //function : Bezier
1194 //purpose  :
1195 //=======================================================================
1196
1197 Handle(Geom2d_BezierCurve)  ProjLib_ComputeApprox::Bezier()  const
1198
1199 {
1200   return myBezier ;
1201 }
1202
1203
1204 //=======================================================================
1205 //function : Tolerance
1206 //purpose  :
1207 //=======================================================================
1208
1209 Standard_Real ProjLib_ComputeApprox::Tolerance() const
1210 {
1211   return myTolerance;
1212 }
1213
1214