src/ProjLib/ProjLib_ComputeApprox.cxx

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