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