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