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