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