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