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