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