src/ProjLib/ProjLib_ComputeApproxOnPolarSurface.cxx

   1 // Created by: Bruno DUMORTIER
   2 // Copyright (c) 1995-1999 Matra Datavision
   3 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #include <ProjLib_ComputeApproxOnPolarSurface.hxx>
  17 #include <AppCont_Function2d.hxx>
  18 #include <ElSLib.hxx>
  19 #include <ElCLib.hxx>
  20 #include <BSplCLib.hxx>
  21 #include <PLib.hxx>
  22 #include <Standard_NoSuchObject.hxx>
  23 #include <Geom_UndefinedDerivative.hxx>
  24 #include <gp_Trsf.hxx>
  25 #include <Precision.hxx>
  26 #include <Approx_FitAndDivide2d.hxx>
  27 #include <math.hxx>
  28 #include <AppParCurves_MultiCurve.hxx>
  29 #include <Geom_Surface.hxx>
  30 #include <Geom2d_BSplineCurve.hxx>
  31 #include <Geom2d_BezierCurve.hxx>
  32 #include <Geom2d_Line.hxx>
  33 #include <Geom2d_Circle.hxx>
  34 #include <Geom2d_Ellipse.hxx>
  35 #include <Geom2d_Hyperbola.hxx>
  36 #include <Geom2d_Parabola.hxx>
  37 #include <Geom2d_TrimmedCurve.hxx>
  38 #include <Geom_BSplineSurface.hxx>
  39 #include <Geom_BezierSurface.hxx>
  40 #include <Geom_BSplineCurve.hxx>
  41 #include <Geom_BezierCurve.hxx>
  42 #include <Geom_TrimmedCurve.hxx>
  43
  44 #include <TColgp_Array1OfPnt2d.hxx>
  45 #include <TColgp_Array2OfPnt2d.hxx>
  46 #include <TColgp_Array1OfPnt.hxx>
  47 #include <TColgp_SequenceOfPnt2d.hxx>
  48 #include <TColStd_Array1OfReal.hxx>
  49 #include <TColStd_Array1OfInteger.hxx>
  50 #include <TColStd_SequenceOfReal.hxx>
  51 #include <TColStd_ListOfTransient.hxx>
  52
  53 #include <GeomAbs_SurfaceType.hxx>
  54 #include <GeomAbs_CurveType.hxx>
  55 #include <Adaptor3d_Surface.hxx>
  56 #include <Adaptor3d_Curve.hxx>
  57 #include <Adaptor3d_HSurface.hxx>
  58 #include <Adaptor3d_HCurve.hxx>
  59 #include <Adaptor2d_HCurve2d.hxx>
  60 #include <Geom2dAdaptor_Curve.hxx>
  61 #include <Geom2dAdaptor_HCurve.hxx>
  62 #include <GeomAdaptor_HCurve.hxx>
  63 #include <GeomAdaptor.hxx>
  64 #include <GeomAdaptor_Surface.hxx>
  65 #include <TColgp_SequenceOfPnt.hxx>
  66
  67 #include <gp_Pnt.hxx>
  68 #include <gp_Pnt2d.hxx>
  69 #include <gp_Vec2d.hxx>
  70 #include <Extrema_GenLocateExtPS.hxx>
  71 #include <Extrema_ExtPS.hxx>
  72 #include <GCPnts_QuasiUniformAbscissa.hxx>
  73 #include <Standard_DomainError.hxx>
  74 //#include <GeomLib_IsIso.hxx>
  75 //#include <GeomLib_CheckSameParameter.hxx>
  76
  77 #ifdef DEB
  78 #ifdef DRAW
  79 #include <DrawTrSurf.hxx>
  80 #endif
  81 //static Standard_Integer compteur = 0;
  82 #endif
  83
  84 //=======================================================================
  85 //function : Value
  86 //purpose  : (OCC217 - apo)- Compute Point2d that project on polar surface(<Surf>) 3D<Curve>
  87 //            <InitCurve2d> use for calculate start 2D point.
  88 //=======================================================================
  89
  90 static gp_Pnt2d Function_Value(const Standard_Real U,
  91                                const Handle(Adaptor3d_HSurface)& Surf,
  92                                const Handle(Adaptor3d_HCurve)& Curve,
  93                                const Handle(Adaptor2d_HCurve2d)& InitCurve2d,
  94                                //OCC217
  95                                const Standard_Real DistTol3d, const Standard_Real tolU, const Standard_Real tolV)
  96                                //const Standard_Real Tolerance)
  97 {
  98   //OCC217
  99   //Standard_Real Tol3d = 100*Tolerance;
 100
 101   gp_Pnt2d p2d = InitCurve2d->Value(U) ;
 102   gp_Pnt p = Curve->Value(U);
 103 //  Curve->D0(U,p) ;
 104   Standard_Real Uinf, Usup, Vinf, Vsup;
 105   Uinf = Surf->Surface().FirstUParameter();
 106   Usup = Surf->Surface().LastUParameter();
 107   Vinf = Surf->Surface().FirstVParameter();
 108   Vsup = Surf->Surface().LastVParameter();
 109   Standard_Integer decalU = 0, decalV = 0;
 110   Standard_Real U0 = p2d.X(), V0 = p2d.Y();
 111
 112   GeomAbs_SurfaceType Type = Surf->GetType();
 113   if((Type != GeomAbs_BSplineSurface) &&
 114      (Type != GeomAbs_BezierSurface)  &&
 115      (Type != GeomAbs_OffsetSurface)    ) {
 116     Standard_Real S = 0., T = 0.;
 117     switch (Type) {
 118 //    case GeomAbs_Plane:
 119 //      {
 120 //      gp_Pln Plane = Surf->Plane();
 121 //      ElSLib::Parameters( Plane, p, S, T);
 122 //      break;
 123 //      }
 124     case GeomAbs_Cylinder:
 125       {
 126         gp_Cylinder Cylinder = Surf->Cylinder();
 127         ElSLib::Parameters( Cylinder, p, S, T);
 128         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 129         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 130         S += decalU*2*M_PI;
 131         break;
 132       }
 133     case GeomAbs_Cone:
 134       {
 135         gp_Cone Cone = Surf->Cone();
 136         ElSLib::Parameters( Cone, p, S, T);
 137         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 138         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 139         S += decalU*2*M_PI;
 140         break;
 141       }
 142     case GeomAbs_Sphere:
 143       {
 144         gp_Sphere Sphere = Surf->Sphere();
 145         ElSLib::Parameters( Sphere, p, S, T);
 146         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 147         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 148         S += decalU*2*M_PI;
 149         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
 150         if(V0 > (Vsup+(Vsup-Vinf))) decalV =  int((V0 - Vsup+(Vsup-Vinf))/(2*M_PI))+1;
 151         T += decalV*2*M_PI;
 152         if(0.4*M_PI < Abs(U0 - S) && Abs(U0 - S) < 1.6*M_PI) {
 153           T = M_PI - T;
 154           if(U0 < S)
 155             S -= M_PI;
 156           else
 157             S += M_PI;
 158         }
 159         break;
 160       }
 161     case GeomAbs_Torus:
 162       {
 163         gp_Torus Torus = Surf->Torus();
 164         ElSLib::Parameters( Torus, p, S, T);
 165         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 166         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 167         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
 168         if(V0 > Vsup) decalV =  int((V0 - Vsup)/(2*M_PI))+1;
 169         S += decalU*2*M_PI; T += decalV*2*M_PI;
 170         break;
 171       }
 172     default:
 173       Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::Value");
 174     }
 175     return gp_Pnt2d(S, T);
 176   }
 177
 178   //////////////////
 179   Standard_Real Dist2Min = RealLast();
 180   //OCC217
 181   //Standard_Real tolU,tolV ;
 182   //tolU = Tolerance;
 183   //tolV = Tolerance;
 184
 185   Standard_Real uperiod =0, vperiod = 0, u, v;
 186   // U0 and V0 are the points within the initialized period
 187   // (periode with u and v),
 188   // U1 and V1 are the points for construction of tops
 189
 190   if(Surf->IsUPeriodic() || Surf->IsUClosed()) {
 191     uperiod =  Surf->LastUParameter() - Surf->FirstUParameter();
 192   }
 193   if(Surf->IsVPeriodic() || Surf->IsVClosed()) {
 194     vperiod = Surf->LastVParameter() - Surf->FirstVParameter();
 195   }
 196   if(U0 < Uinf) {
 197     if(!uperiod)
 198       U0 = Uinf;
 199     else {
 200       decalU = int((Uinf - U0)/uperiod)+1;
 201       U0 += decalU*uperiod;
 202     }
 203   }
 204   if(U0 > Usup) {
 205     if(!uperiod)
 206       U0 = Usup;
 207     else {
 208       decalU = -(int((U0 - Usup)/uperiod)+1);
 209       U0 += decalU*uperiod;
 210     }
 211   }
 212   if(V0 < Vinf) {
 213     if(!vperiod)
 214       V0 = Vinf;
 215     else {
 216       decalV = int((Vinf - V0)/vperiod)+1;
 217       V0 += decalV*vperiod;
 218     }
 219   }
 220   if(V0 > Vsup) {
 221     if(!vperiod)
 222       V0 = Vsup;
 223     else {
 224       decalV = -int((V0 - Vsup)/vperiod)-1;
 225       V0 += decalV*vperiod;
 226     }
 227   }
 228
 229   // The surface around U0 is reduced
 230   Standard_Real uLittle = (Usup - Uinf)/10, vLittle = (Vsup - Vinf)/10;
 231   Standard_Real uInfLi = 0, vInfLi = 0,uSupLi = 0, vSupLi = 0;
 232   if((U0 - Uinf) > uLittle) uInfLi = U0 - uLittle; else uInfLi = Uinf;
 233   if((V0 - Vinf) > vLittle) vInfLi = V0 - vLittle; else vInfLi = Vinf;
 234   if((Usup - U0) > uLittle) uSupLi = U0 + uLittle; else uSupLi = Usup;
 235   if((Vsup - V0) > vLittle) vSupLi = V0 + vLittle; else vSupLi = Vsup;
 236
 237   //  const Adaptor3d_Surface GAS = Surf->Surface();
 238
 239   GeomAdaptor_Surface SurfLittle;
 240   if (Type == GeomAbs_BSplineSurface) {
 241     Handle(Geom_Surface) GBSS(Surf->Surface().BSpline());
 242     SurfLittle.Load(GBSS, uInfLi, uSupLi, vInfLi, vSupLi);
 243   }
 244   else if (Type == GeomAbs_BezierSurface) {
 245     Handle(Geom_Surface) GS(Surf->Surface().Bezier());
 246     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
 247   }
 248   else if (Type == GeomAbs_OffsetSurface) {
 249     Handle(Geom_Surface) GS = GeomAdaptor::MakeSurface(Surf->Surface());
 250     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
 251   }
 252   else {
 253     Standard_NoSuchObject::Raise("");
 254   }
 255
 256   Extrema_GenLocateExtPS  locext(p, SurfLittle, U0, V0, tolU, tolV);
 257   if (locext.IsDone()) {
 258     Dist2Min = locext.SquareDistance();
 259     if (Dist2Min < DistTol3d * DistTol3d) {
 260       (locext.Point()).Parameter(u, v);
 261       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
 262       return pnt;
 263     }
 264   }
 265
 266   Extrema_ExtPS  ext(p, SurfLittle, tolU, tolV) ;
 267   if (ext.IsDone() && ext.NbExt()>=1 ) {
 268     Dist2Min = ext.SquareDistance(1);
 269     Standard_Integer GoodValue = 1;
 270     for ( Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ )
 271       if( Dist2Min > ext.SquareDistance(i)) {
 272         Dist2Min = ext.SquareDistance(i);
 273         GoodValue = i;
 274       }
 275     if (Dist2Min < DistTol3d * DistTol3d) {
 276       (ext.Point(GoodValue)).Parameter(u,v);
 277       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
 278       return pnt;
 279     }
 280   }
 281
 282   return p2d;
 283 }
 284
 285
 286 //=======================================================================
 287 //function : ProjLib_PolarFunction
 288 //purpose  : (OCC217 - apo)- This class produce interface to call "gp_Pnt2d Function_Value(...)"
 289 //=======================================================================
 290
 291 class ProjLib_PolarFunction : public AppCont_Function2d
 292 {
 293   Handle(Adaptor3d_HCurve)           myCurve;
 294   Handle(Adaptor2d_HCurve2d)         myInitialCurve2d ;
 295   Handle(Adaptor3d_HSurface)         mySurface ;
 296   //OCC217
 297   Standard_Real                      myTolU, myTolV;
 298   Standard_Real                      myDistTol3d;
 299   //Standard_Real                    myTolerance ;
 300
 301   public :
 302
 303   ProjLib_PolarFunction(const Handle(Adaptor3d_HCurve) & C,
 304                         const Handle(Adaptor3d_HSurface)& Surf,
 305                         const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 306                         //OCC217
 307                         const Standard_Real Tol3d) :
 308                         //const Standard_Real Tolerance) :
 309   myCurve(C),
 310   myInitialCurve2d(InitialCurve2d),
 311   mySurface(Surf),
 312   //OCC217
 313   myTolU(Surf->UResolution(Tol3d)),
 314   myTolV(Surf->VResolution(Tol3d)),
 315   myDistTol3d(100.0*Tol3d){} ;
 316   //myTolerance(Tolerance){} ;
 317
 318   ~ProjLib_PolarFunction() {}
 319
 320   Standard_Real FirstParameter() const
 321   {return (myCurve->FirstParameter()/*+1.e-9*/);}
 322
 323   Standard_Real LastParameter() const
 324   {return (myCurve->LastParameter()/*-1.e-9*/);}
 325
 326   gp_Pnt2d Value(const Standard_Real t) const {
 327     return Function_Value
 328       (t,mySurface,myCurve,myInitialCurve2d,myDistTol3d,myTolU,myTolV) ;  //OCC217
 329       //(t,mySurface,myCurve,myInitialCurve2d,myTolerance) ;
 330   }
 331
 332 //  Standard_Boolean D1(const Standard_Real t, gp_Pnt2d& P, gp_Vec2d& V) const
 333   Standard_Boolean D1(const Standard_Real , gp_Pnt2d& , gp_Vec2d& ) const
 334     {return Standard_False;}
 335 };
 336
 337 //=======================================================================
 338 //function : ProjLib_ComputeApproxOnPolarSurface
 339 //purpose  :
 340 //=======================================================================
 341
 342 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface(){}
 343
 344
 345 //=======================================================================
 346 //function : ProjLib_ComputeApproxOnPolarSurface
 347 //purpose  :
 348 //=======================================================================
 349
 350 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 351 (const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 352  const Handle(Adaptor3d_HCurve)& Curve,
 353  const Handle(Adaptor3d_HSurface)& S,
 354  const Standard_Real tol3d):myProjIsDone(Standard_False)  //OCC217
 355  //const Standard_Real tolerance):myProjIsDone(Standard_False)
 356 {
 357   myTolerance = tol3d; //OCC217
 358   //myTolerance = Max(Tolerance,Precision::PApproximation());
 359   myBSpline = Perform(InitialCurve2d,Curve,S);
 360 }
 361 //=======================================================================
 362 //function : ProjLib_ComputeApproxOnPolarSurface
 363 //purpose  : Process the case of sewing
 364 //=======================================================================
 365
 366 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 367 (const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 368  const Handle(Adaptor2d_HCurve2d)& InitialCurve2dBis,
 369  const Handle(Adaptor3d_HCurve)& Curve,
 370  const Handle(Adaptor3d_HSurface)& S,
 371  const Standard_Real tol3d):myProjIsDone(Standard_False)  //OCC217
 372  //const Standard_Real tolerance):myProjIsDone(Standard_False)
 373 {// InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing
 374   myTolerance = tol3d; //OCC217
 375   //myTolerance = Max(tolerance,Precision::PApproximation());
 376   Handle(Geom2d_BSplineCurve) bsc = Perform(InitialCurve2d,Curve,S);
 377   if(myProjIsDone) {
 378     gp_Pnt2d P2dproj, P2d, P2dBis;
 379     P2dproj = bsc->StartPoint();
 380     P2d    = InitialCurve2d->Value(InitialCurve2d->FirstParameter());
 381     P2dBis = InitialCurve2dBis->Value(InitialCurve2dBis->FirstParameter());
 382
 383     Standard_Real Dist, DistBis;
 384     Dist    = P2dproj.Distance(P2d);
 385     DistBis = P2dproj.Distance(P2dBis);
 386     if( Dist < DistBis)  {
 387       // myBSpline2d is the pcurve that is found. It is translated to obtain myCurve2d
 388       myBSpline = bsc;
 389       Handle(Geom2d_Geometry) GG = myBSpline->Translated(P2d, P2dBis);
 390       my2ndCurve = Handle(Geom2d_Curve)::DownCast(GG);
 391     }
 392     else {
 393       my2ndCurve = bsc;
 394       Handle(Geom2d_Geometry) GG = my2ndCurve->Translated(P2dBis, P2d);
 395       myBSpline = Handle(Geom2d_BSplineCurve)::DownCast(GG);
 396     }
 397   }
 398 }
 399
 400 //=======================================================================
 401 //function : ProjLib_ComputeApproxOnPolarSurface
 402 //purpose  : case without curve of initialization
 403 //=======================================================================
 404
 405 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 406 (const Handle(Adaptor3d_HCurve)& Curve,
 407  const Handle(Adaptor3d_HSurface)& S,
 408  const Standard_Real tol3d):myProjIsDone(Standard_False)  //OCC217
 409  //const Standard_Real tolerance):myProjIsDone(Standard_False)
 410 {
 411   myTolerance = tol3d; //OCC217
 412   //myTolerance = Max(tolerance,Precision::PApproximation());
 413   const Handle(Adaptor2d_HCurve2d) InitCurve2d ;
 414   myBSpline = Perform(InitCurve2d,Curve,S);
 415 }
 416
 417 static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
 418                                           Handle(Geom2d_BSplineCurve) C2)
 419 {
 420   Standard_Integer deg, deg1, deg2;
 421   deg1 = C1->Degree();
 422   deg2 = C2->Degree();
 423
 424   if ( deg1 < deg2) {
 425     C1->IncreaseDegree(deg2);
 426     deg = deg2;
 427   }
 428   else if ( deg2 < deg1) {
 429     C2->IncreaseDegree(deg1);
 430     deg = deg1;
 431   }
 432   else deg = deg1;
 433
 434   Standard_Integer np1,np2,nk1,nk2,np,nk;
 435   np1 = C1->NbPoles();
 436   nk1 = C1->NbKnots();
 437   np2 = C2->NbPoles();
 438   nk2 = C2->NbKnots();
 439   nk = nk1 + nk2 -1;
 440   np = np1 + np2 -1;
 441
 442   TColStd_Array1OfReal    K1(1,nk1); C1->Knots(K1);
 443   TColStd_Array1OfInteger M1(1,nk1); C1->Multiplicities(M1);
 444   TColgp_Array1OfPnt2d    P1(1,np1); C1->Poles(P1);
 445   TColStd_Array1OfReal    K2(1,nk2); C2->Knots(K2);
 446   TColStd_Array1OfInteger M2(1,nk2); C2->Multiplicities(M2);
 447   TColgp_Array1OfPnt2d    P2(1,np2); C2->Poles(P2);
 448
 449   // Compute the new BSplineCurve
 450   TColStd_Array1OfReal    K(1,nk);
 451   TColStd_Array1OfInteger M(1,nk);
 452   TColgp_Array1OfPnt2d    P(1,np);
 453
 454   Standard_Integer i, count = 0;
 455   // Set Knots and Mults
 456   for ( i = 1; i <= nk1; i++) {
 457     count++;
 458     K(count) = K1(i);
 459     M(count) = M1(i);
 460   }
 461   M(count) = deg;
 462   for ( i = 2; i <= nk2; i++) {
 463     count++;
 464     K(count) = K2(i);
 465     M(count) = M2(i);
 466   }
 467   // Set the Poles
 468   count = 0;
 469   for (i = 1; i <= np1; i++) {
 470     count++;
 471     P(count) = P1(i);
 472   }
 473   for (i = 2; i <= np2; i++) {
 474     count++;
 475     P(count) = P2(i);
 476   }
 477
 478   Handle(Geom2d_BSplineCurve) BS =
 479     new Geom2d_BSplineCurve(P,K,M,deg);
 480   return BS;
 481 }
 482
 483
 484 //=======================================================================
 485 //function : Perform
 486 //purpose  :
 487 //=======================================================================
 488 Handle(Geom2d_BSplineCurve) ProjLib_ComputeApproxOnPolarSurface::Perform
 489 (const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 490  const Handle(Adaptor3d_HCurve)& Curve,
 491  const Handle(Adaptor3d_HSurface)& S)
 492 {
 493   //OCC217
 494   Standard_Real Tol3d = myTolerance;
 495   Standard_Real ParamTol = Precision::PApproximation();
 496
 497   Handle(Adaptor2d_HCurve2d) AHC2d = InitialCurve2d;
 498   Handle(Adaptor3d_HCurve) AHC = Curve;
 499
 500 // if the curve 3d is a BSpline with degree C0, it is cut into sections with degree C1
 501 // -> bug cts18237
 502   GeomAbs_CurveType typeCurve = Curve->GetType();
 503   if(typeCurve == GeomAbs_BSplineCurve) {
 504     TColStd_ListOfTransient LOfBSpline2d;
 505     Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
 506     Standard_Integer nbInter = Curve->NbIntervals(GeomAbs_C1);
 507     if(nbInter > 1) {
 508       Standard_Integer i, j;
 509       Handle(Geom_TrimmedCurve) GTC;
 510       Handle(Geom2d_TrimmedCurve) G2dTC;
 511       TColStd_Array1OfReal Inter(1,nbInter+1);
 512       Curve->Intervals(Inter,GeomAbs_C1);
 513       Standard_Real firstinter = Inter.Value(1), secondinter = Inter.Value(2);
 514       // initialization 3d
 515       GTC = new Geom_TrimmedCurve(BSC, firstinter, secondinter);
 516       AHC = new GeomAdaptor_HCurve(GTC);
 517
 518       // if there is an initialization curve:
 519       // - either this is a BSpline C0, with discontinuity at the same parameters of nodes
 520       // and the sections C1 are taken
 521       // - or this is a curve C1 and the sections of intrest are taken otherwise the curve is created.
 522
 523       // initialization 2d
 524       Standard_Integer nbInter2d;
 525       Standard_Boolean C2dIsToCompute;
 526       C2dIsToCompute = InitialCurve2d.IsNull();
 527       Handle(Geom2d_BSplineCurve) BSC2d;
 528       Handle(Geom2d_Curve) G2dC;
 529
 530       if(!C2dIsToCompute) {
 531         nbInter2d = InitialCurve2d->NbIntervals(GeomAbs_C1);
 532         TColStd_Array1OfReal Inter2d(1,nbInter2d+1);
 533         InitialCurve2d->Intervals(Inter2d,GeomAbs_C1);
 534         j = 1;
 535         for(i = 1,j = 1;i <= nbInter;i++)
 536           if(Abs(Inter.Value(i) - Inter2d.Value(j)) < ParamTol) { //OCC217
 537           //if(Abs(Inter.Value(i) - Inter2d.Value(j)) < myTolerance) {
 538             if (j > nbInter2d) break;
 539             j++;
 540           }
 541         if(j != (nbInter2d+1)) {
 542           C2dIsToCompute = Standard_True;
 543         }
 544       }
 545
 546       if(C2dIsToCompute) {
 547         AHC2d = BuildInitialCurve2d(AHC, S);
 548       }
 549       else {
 550         typeCurve = InitialCurve2d->GetType();
 551         switch (typeCurve) {
 552         case GeomAbs_Line: {
 553           G2dC = new Geom2d_Line(InitialCurve2d->Line());
 554           break;
 555         }
 556         case GeomAbs_Circle: {
 557           G2dC = new Geom2d_Circle(InitialCurve2d->Circle());
 558           break;
 559         }
 560         case GeomAbs_Ellipse: {
 561           G2dC = new Geom2d_Ellipse(InitialCurve2d->Ellipse());
 562           break;
 563         }
 564         case GeomAbs_Hyperbola: {
 565           G2dC = new Geom2d_Hyperbola(InitialCurve2d->Hyperbola());
 566           break;
 567         }
 568         case GeomAbs_Parabola: {
 569           G2dC = new Geom2d_Parabola(InitialCurve2d->Parabola());
 570           break;
 571         }
 572         case GeomAbs_BezierCurve: {
 573           G2dC = InitialCurve2d->Bezier();
 574           break;
 575         }
 576         case GeomAbs_BSplineCurve: {
 577           G2dC = InitialCurve2d->BSpline();
 578           break;
 579         }
 580         case GeomAbs_OtherCurve:
 581         default:
 582           break;
 583         }
 584         gp_Pnt2d fp2d = G2dC->Value(firstinter), lp2d = G2dC->Value(secondinter);
 585         gp_Pnt fps, lps, fpc, lpc;
 586         S->D0(fp2d.X(), fp2d.Y(), fps);
 587         S->D0(lp2d.X(), lp2d.Y(), lps);
 588         Curve->D0(firstinter, fpc);
 589         Curve->D0(secondinter, lpc);
 590         //OCC217
 591         if((fps.IsEqual(fpc, Tol3d)) &&
 592            (lps.IsEqual(lpc, Tol3d))) {
 593         //if((fps.IsEqual(fpc, myTolerance)) &&
 594         //   (lps.IsEqual(lpc, myTolerance))) {
 595           G2dTC = new Geom2d_TrimmedCurve(G2dC, firstinter, secondinter);
 596           Geom2dAdaptor_Curve G2dAC(G2dTC);
 597           AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
 598           myProjIsDone = Standard_True;
 599         }
 600         else {
 601           AHC2d = BuildInitialCurve2d(AHC, S);
 602           C2dIsToCompute = Standard_True;
 603         }
 604       }
 605
 606       if(myProjIsDone) {
 607         BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 608         if(BSC2d.IsNull()) return Handle(Geom2d_BSplineCurve)(); //IFV
 609         LOfBSpline2d.Append(BSC2d);
 610       }
 611       else {
 612         return Handle(Geom2d_BSplineCurve)();
 613       }
 614
 615
 616
 617       Standard_Real iinter, ip1inter;
 618       Standard_Integer nbK2d, deg;
 619       nbK2d = BSC2d->NbKnots(); deg = BSC2d->Degree();
 620
 621       for(i = 2;i <= nbInter;i++) {
 622         iinter = Inter.Value(i);
 623         ip1inter = Inter.Value(i+1);
 624         // general case 3d
 625         GTC->SetTrim(iinter, ip1inter);
 626         AHC = new GeomAdaptor_HCurve(GTC);
 627
 628         // general case 2d
 629         if(C2dIsToCompute) {
 630           AHC2d = BuildInitialCurve2d(AHC, S);
 631         }
 632         else {
 633           gp_Pnt2d fp2d = G2dC->Value(iinter), lp2d = G2dC->Value(ip1inter);
 634           gp_Pnt fps, lps, fpc, lpc;
 635           S->D0(fp2d.X(), fp2d.Y(), fps);
 636           S->D0(lp2d.X(), lp2d.Y(), lps);
 637           Curve->D0(iinter, fpc);
 638           Curve->D0(ip1inter, lpc);
 639           //OCC217
 640           if((fps.IsEqual(fpc, Tol3d)) &&
 641              (lps.IsEqual(lpc, Tol3d))) {
 642           //if((fps.IsEqual(fpc, myTolerance)) &&
 643           //   (lps.IsEqual(lpc, myTolerance))) {
 644             G2dTC->SetTrim(iinter, ip1inter);
 645             Geom2dAdaptor_Curve G2dAC(G2dTC);
 646             AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
 647             myProjIsDone = Standard_True;
 648           }
 649           else {
 650             AHC2d = BuildInitialCurve2d(AHC, S);
 651           }
 652         }
 653         if(myProjIsDone) {
 654           BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 655           if(BSC2d.IsNull()) {
 656             return Handle(Geom2d_BSplineCurve)();
 657           }
 658           LOfBSpline2d.Append(BSC2d);
 659           nbK2d += BSC2d->NbKnots() - 1;
 660           deg = Max(deg, BSC2d->Degree());
 661         }
 662         else {
 663           return Handle(Geom2d_BSplineCurve)();
 664         }
 665       }
 666
 667       Standard_Integer NbC = LOfBSpline2d.Extent();
 668       Handle(Geom2d_BSplineCurve) CurBS;
 669       CurBS = Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
 670       LOfBSpline2d.RemoveFirst();
 671       for (Standard_Integer ii = 2; ii <= NbC; ii++) {
 672         Handle(Geom2d_BSplineCurve) BS =
 673           Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
 674         CurBS = Concat(CurBS,BS);
 675         LOfBSpline2d.RemoveFirst();
 676       }
 677       return CurBS;
 678     }
 679   }
 680
 681   if(InitialCurve2d.IsNull()) {
 682     AHC2d = BuildInitialCurve2d(Curve, S);
 683     if(!myProjIsDone)
 684       return Handle(Geom2d_BSplineCurve)();
 685   }
 686   return ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 687 }
 688
 689 //=======================================================================
 690 //function : ProjLib_BuildInitialCurve2d
 691 //purpose  :
 692 //=======================================================================
 693
 694 Handle(Adaptor2d_HCurve2d)
 695      ProjLib_ComputeApproxOnPolarSurface::
 696      BuildInitialCurve2d(const Handle(Adaptor3d_HCurve)&   Curve,
 697                          const Handle(Adaptor3d_HSurface)& Surf)
 698 {
 699   //  discretize the Curve with quasiuniform deflection
 700   //  density at least NbOfPnts points
 701   myProjIsDone = Standard_False;
 702
 703   //OCC217
 704   Standard_Real Tol3d = myTolerance;
 705   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
 706   Standard_Real DistTol3d = 100.0*Tol3d;
 707
 708   Standard_Real uperiod = 0., vperiod = 0.;
 709   if(Surf->IsUPeriodic() || Surf->IsUClosed())
 710     uperiod = Surf->LastUParameter() - Surf->FirstUParameter();
 711
 712   if(Surf->IsVPeriodic() || Surf->IsVClosed())
 713     vperiod = Surf->LastVParameter() - Surf->FirstVParameter();
 714
 715
 716   // NO myTol is Tol2d !!!!
 717   //Standard_Real TolU = myTolerance, TolV = myTolerance;
 718   //Standard_Real Tol3d = 100*myTolerance; // At random Balthazar.
 719
 720   Standard_Integer NbOfPnts = 61;
 721   GCPnts_QuasiUniformAbscissa QUA(Curve->GetCurve(),NbOfPnts);
 722   TColgp_Array1OfPnt Pts(1,NbOfPnts);
 723   TColStd_Array1OfReal Param(1,NbOfPnts);
 724   Standard_Integer i, j;
 725   for( i = 1; i <= NbOfPnts ; i++ ) {
 726     Param(i) = QUA.Parameter(i);
 727     Pts(i) = Curve->Value(Param(i));
 728   }
 729
 730   TColgp_Array1OfPnt2d Pts2d(1,NbOfPnts);
 731   TColStd_Array1OfInteger Mult(1,NbOfPnts);
 732   Mult.Init(1);
 733   Mult(1) = Mult(NbOfPnts) = 2;
 734
 735   Standard_Real Uinf, Usup, Vinf, Vsup;
 736   Uinf = Surf->Surface().FirstUParameter();
 737   Usup = Surf->Surface().LastUParameter();
 738   Vinf = Surf->Surface().FirstVParameter();
 739   Vsup = Surf->Surface().LastVParameter();
 740   GeomAbs_SurfaceType Type = Surf->GetType();
 741   if((Type != GeomAbs_BSplineSurface) && (Type != GeomAbs_BezierSurface) &&
 742      (Type != GeomAbs_OffsetSurface)) {
 743     Standard_Real S, T;
 744 //    Standard_Integer usens = 0, vsens = 0;
 745     // to know the position relatively to the period
 746     switch (Type) {
 747 //    case GeomAbs_Plane:
 748 //      {
 749 //      gp_Pln Plane = Surf->Plane();
 750 //      for ( i = 1 ; i <= NbOfPnts ; i++) {
 751 //        ElSLib::Parameters( Plane, Pts(i), S, T);
 752 //        Pts2d(i).SetCoord(S,T);
 753 //      }
 754 //      myProjIsDone = Standard_True;
 755 //      break;
 756 //      }
 757     case GeomAbs_Cylinder:
 758       {
 759 //      Standard_Real Sloc, Tloc;
 760         Standard_Real Sloc;
 761         Standard_Integer usens = 0;
 762         gp_Cylinder Cylinder = Surf->Cylinder();
 763         ElSLib::Parameters( Cylinder, Pts(1), S, T);
 764         Pts2d(1).SetCoord(S,T);
 765         for ( i = 2 ; i <= NbOfPnts ; i++) {
 766           Sloc = S;
 767           ElSLib::Parameters( Cylinder, Pts(i), S, T);
 768           if(Abs(Sloc - S) > M_PI) {
 769             if(Sloc > S)
 770               usens++;
 771             else
 772               usens--;
 773           }
 774           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
 775         }
 776         myProjIsDone = Standard_True;
 777         break;
 778       }
 779     case GeomAbs_Cone:
 780       {
 781 //      Standard_Real Sloc, Tloc;
 782         Standard_Real Sloc;
 783         Standard_Integer usens = 0;
 784         gp_Cone Cone = Surf->Cone();
 785         ElSLib::Parameters( Cone, Pts(1), S, T);
 786         Pts2d(1).SetCoord(S,T);
 787         for ( i = 2 ; i <= NbOfPnts ; i++) {
 788           Sloc = S;
 789           ElSLib::Parameters( Cone, Pts(i), S, T);
 790           if(Abs(Sloc - S) > M_PI) {
 791             if(Sloc > S)
 792               usens++;
 793             else
 794               usens--;
 795           }
 796           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
 797         }
 798         myProjIsDone = Standard_True;
 799         break;
 800       }
 801     case GeomAbs_Sphere:
 802       {
 803         Standard_Real Sloc, Tloc;
 804         Standard_Integer usens = 0, vsens = 0; //usens steps by half-period
 805         Standard_Boolean vparit = Standard_False;
 806         gp_Sphere Sphere = Surf->Sphere();
 807         ElSLib::Parameters( Sphere, Pts(1), S, T);
 808         Pts2d(1).SetCoord(S,T);
 809         for ( i = 2 ; i <= NbOfPnts ; i++) {
 810           Sloc = S;Tloc = T;
 811           ElSLib::Parameters( Sphere, Pts(i), S, T);
 812           if(1.6*M_PI < Abs(Sloc - S)) {
 813             if(Sloc > S)
 814               usens += 2;
 815             else
 816               usens -= 2;
 817     }
 818           if(1.6*M_PI > Abs(Sloc - S) && Abs(Sloc - S) > 0.4*M_PI) {
 819             vparit = !vparit;
 820             if(Sloc > S)
 821               usens++;
 822             else
 823               usens--;
 824             if(Abs(Tloc - Vsup) < (Vsup - Vinf)/5)
 825               vsens++;
 826             else
 827               vsens--;
 828           }
 829           if(vparit) {
 830             Pts2d(i).SetCoord(S+usens*M_PI,(M_PI - T)*(vsens-1));
 831           }
 832           else {
 833             Pts2d(i).SetCoord(S+usens*M_PI,T+vsens*M_PI);
 834
 835           }
 836         }
 837         myProjIsDone = Standard_True;
 838         break;
 839       }
 840     case GeomAbs_Torus:
 841       {
 842         Standard_Real Sloc, Tloc;
 843         Standard_Integer usens = 0, vsens = 0;
 844         gp_Torus Torus = Surf->Torus();
 845         ElSLib::Parameters( Torus, Pts(1), S, T);
 846         Pts2d(1).SetCoord(S,T);
 847         for ( i = 2 ; i <= NbOfPnts ; i++) {
 848           Sloc = S; Tloc = T;
 849           ElSLib::Parameters( Torus, Pts(i), S, T);
 850           if(Abs(Sloc - S) > M_PI) {
 851             if(Sloc > S)
 852               usens++;
 853             else
 854               usens--;
 855     }
 856           if(Abs(Tloc - T) > M_PI) {
 857             if(Tloc > T)
 858               vsens++;
 859             else
 860               vsens--;
 861     }
 862           Pts2d(i).SetCoord(S+usens*2*M_PI,T+vsens*2*M_PI);
 863         }
 864         myProjIsDone = Standard_True;
 865         break;
 866       }
 867     default:
 868       Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::BuildInitialCurve2d");
 869     }
 870   }
 871   else {
 872     myProjIsDone = Standard_False;
 873     Standard_Real Dist2Min = 1.e+200, u = 0., v = 0.;
 874     gp_Pnt pntproj;
 875
 876     TColgp_SequenceOfPnt2d Sols;
 877     Standard_Boolean areManyZeros = Standard_False;
 878
 879     Curve->D0(Param.Value(1), pntproj) ;
 880     Extrema_ExtPS  aExtPS(pntproj, Surf->Surface(), TolU, TolV) ;
 881     Standard_Real aMinSqDist = RealLast();
 882     if (aExtPS.IsDone())
 883     {
 884       for (i = 1; i <= aExtPS.NbExt(); i++)
 885       {
 886         Standard_Real aSqDist = aExtPS.SquareDistance(i);
 887         if (aSqDist < aMinSqDist)
 888           aMinSqDist = aSqDist;
 889       }
 890     }
 891     if (aMinSqDist > DistTol3d * DistTol3d) //try to project with less tolerance
 892     {
 893       TolU = Min(TolU, Precision::PConfusion());
 894       TolV = Min(TolV, Precision::PConfusion());
 895       aExtPS.Initialize(Surf->Surface(),
 896                         Surf->Surface().FirstUParameter(), Surf->Surface().LastUParameter(),
 897                         Surf->Surface().FirstVParameter(), Surf->Surface().LastVParameter(),
 898                         TolU, TolV);
 899       aExtPS.Perform(pntproj);
 900     }
 901
 902     if( aExtPS.IsDone() && aExtPS.NbExt() >= 1 ) {
 903
 904       Standard_Integer GoodValue = 1;
 905
 906       for ( i = 1 ; i <= aExtPS.NbExt() ; i++ ) {
 907         if( aExtPS.SquareDistance(i) < DistTol3d * DistTol3d ) {
 908           if( aExtPS.SquareDistance(i) <= 1.e-18 ) {
 909             aExtPS.Point(i).Parameter(u,v);
 910             gp_Pnt2d p2d(u,v);
 911             Standard_Boolean isSame = Standard_False;
 912             for( j = 1; j <= Sols.Length(); j++ ) {
 913               if( p2d.SquareDistance( Sols.Value(j) ) <= 1.e-18 ) {
 914                 isSame = Standard_True;
 915                 break;
 916               }
 917             }
 918             if( !isSame ) Sols.Append( p2d );
 919           }
 920           if( Dist2Min > aExtPS.SquareDistance(i) ) {
 921             Dist2Min = aExtPS.SquareDistance(i);
 922             GoodValue = i;
 923           }
 924         }
 925       }
 926
 927       if( Sols.Length() > 1 ) areManyZeros = Standard_True;
 928
 929       if( Dist2Min <= DistTol3d * DistTol3d) {
 930         if( !areManyZeros ) {
 931           aExtPS.Point(GoodValue).Parameter(u,v);
 932           Pts2d(1).SetCoord(u,v);
 933           myProjIsDone = Standard_True;
 934         }
 935         else {
 936           Standard_Integer nbSols = Sols.Length();
 937           Standard_Real Dist2Max = -1.e+200;
 938           for( i = 1; i <= nbSols; i++ ) {
 939             const gp_Pnt2d& aP1 = Sols.Value(i);
 940             for( j = i+1; j <= nbSols; j++ ) {
 941               const gp_Pnt2d& aP2 = Sols.Value(j);
 942               Standard_Real aDist2 = aP1.SquareDistance(aP2);
 943               if( aDist2 > Dist2Max ) Dist2Max = aDist2;
 944             }
 945           }
 946       Standard_Real aMaxT2 = Max(TolU,TolV);
 947       aMaxT2 *= aMaxT2;
 948           if( Dist2Max > aMaxT2 ) {
 949             Standard_Integer tPp = 0;
 950             for( i = 1; i <= 5; i++ ) {
 951               Standard_Integer nbExtOk = 0;
 952               Standard_Integer indExt = 0;
 953               Standard_Integer iT = 1 + (NbOfPnts - 1)/5*i;
 954               Curve->D0( Param.Value(iT), pntproj );
 955               Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
 956               Dist2Min = 1.e+200;
 957               if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
 958                 for( j = 1 ; j <= aTPS.NbExt() ; j++ ) {
 959                   if( aTPS.SquareDistance(j) < DistTol3d * DistTol3d ) {
 960                     nbExtOk++;
 961                     if( aTPS.SquareDistance(j) < Dist2Min ) {
 962                       Dist2Min = aTPS.SquareDistance(j);
 963                       indExt = j;
 964                     }
 965                   }
 966                 }
 967               }
 968               if( nbExtOk == 1 ) {
 969                 tPp = iT;
 970                 aTPS.Point(indExt).Parameter(u,v);
 971                 break;
 972               }
 973             }
 974
 975             if( tPp != 0 ) {
 976               gp_Pnt2d aPp = gp_Pnt2d(u,v);
 977               gp_Pnt2d aPn;
 978               j = 1;
 979               Standard_Boolean isFound = Standard_False;
 980               while( !isFound ) {
 981                 Curve->D0( Param.Value(tPp+j), pntproj );
 982                 Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
 983                 Dist2Min = 1.e+200;
 984                 Standard_Integer indExt = 0;
 985                 if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
 986                   for( i = 1 ; i <= aTPS.NbExt() ; i++ ) {
 987                     if( aTPS.SquareDistance(i) < DistTol3d * DistTol3d && aTPS.SquareDistance(i) < Dist2Min ) {
 988                       Dist2Min = aTPS.SquareDistance(i);
 989                       indExt = i;
 990                       isFound = Standard_True;
 991                     }
 992                   }
 993                 }
 994                 if( isFound ) {
 995                   aTPS.Point(indExt).Parameter(u,v);
 996                   aPn = gp_Pnt2d(u,v);
 997                   break;
 998                 }
 999                 j++;
1000                 if( (tPp+j) > NbOfPnts ) break;
1001               }
1002
1003               if( isFound ) {
1004                 gp_Vec2d atV(aPp,aPn);
1005                 Standard_Boolean isChosen = Standard_False;
1006                 for( i = 1; i <= nbSols; i++ ) {
1007                   const gp_Pnt2d& aP1 = Sols.Value(i);
1008                   gp_Vec2d asV(aP1,aPp);
1009                   if( asV.Dot(atV) > 0. ) {
1010                     isChosen = Standard_True;
1011                     Pts2d(1).SetCoord(aP1.X(),aP1.Y());
1012                     myProjIsDone = Standard_True;
1013                     break;
1014                   }
1015                 }
1016                 if( !isChosen ) {
1017                   aExtPS.Point(GoodValue).Parameter(u,v);
1018                   Pts2d(1).SetCoord(u,v);
1019                   myProjIsDone = Standard_True;
1020                 }
1021               }
1022               else {
1023                 aExtPS.Point(GoodValue).Parameter(u,v);
1024                 Pts2d(1).SetCoord(u,v);
1025                 myProjIsDone = Standard_True;
1026               }
1027             }
1028             else {
1029               aExtPS.Point(GoodValue).Parameter(u,v);
1030               Pts2d(1).SetCoord(u,v);
1031               myProjIsDone = Standard_True;
1032             }
1033           }
1034           else {
1035             aExtPS.Point(GoodValue).Parameter(u,v);
1036             Pts2d(1).SetCoord(u,v);
1037             myProjIsDone = Standard_True;
1038           }
1039         }
1040       }
1041
1042       //  calculate the following points with GenLocate_ExtPS
1043       // (and store the result and each parameter in a sequence)
1044       Standard_Integer usens = 0, vsens = 0;
1045       // to know the position relatively to the period
1046       Standard_Real U0 = u, V0 = v, U1 = u, V1 = v;
1047       // U0 and V0 are the points in the initialized period
1048       // (period with u and v),
1049       // U1 and V1 are the points for construction of poles
1050
1051       for ( i = 2 ; i <= NbOfPnts ; i++)
1052         if(myProjIsDone) {
1053           myProjIsDone = Standard_False;
1054           Dist2Min = RealLast();
1055           Curve->D0(Param.Value(i), pntproj);
1056           Extrema_GenLocateExtPS  aLocateExtPS
1057             (pntproj, Surf->Surface(), U0, V0, TolU, TolV) ;
1058
1059           if (aLocateExtPS.IsDone())
1060           {
1061             if (aLocateExtPS.SquareDistance() < DistTol3d * DistTol3d)
1062             {  //OCC217
1063               //if (aLocateExtPS.SquareDistance() < Tol3d * Tol3d) {
1064               (aLocateExtPS.Point()).Parameter(U0,V0);
1065               U1 = U0 + usens*uperiod;
1066               V1 = V0 + vsens*vperiod;
1067               Pts2d(i).SetCoord(U1,V1);
1068               myProjIsDone = Standard_True;
1069             }
1070             else
1071             {
1072               Extrema_ExtPS aGlobalExtr(pntproj, Surf->Surface(), TolU, TolV);
1073               if (aGlobalExtr.IsDone())
1074               {
1075                 Standard_Real LocalMinSqDist = RealLast();
1076                 Standard_Integer imin = 0;
1077                 for (Standard_Integer isol = 1; isol <= aGlobalExtr.NbExt(); isol++)
1078                 {
1079                   Standard_Real aSqDist = aGlobalExtr.SquareDistance(isol);
1080                   if (aSqDist < LocalMinSqDist)
1081                   {
1082                     LocalMinSqDist = aSqDist;
1083                     imin = isol;
1084                   }
1085                 }
1086                 if (LocalMinSqDist < DistTol3d * DistTol3d)
1087                 {
1088                   Standard_Real LocalU, LocalV;
1089                   aGlobalExtr.Point(imin).Parameter(LocalU, LocalV);
1090                   if (uperiod > 0. && Abs(U0 - LocalU) >= uperiod/2.)
1091                   {
1092                     if (LocalU > U0)
1093                       usens = -1;
1094                     else
1095                       usens = 1;
1096                   }
1097                   if (vperiod > 0. && Abs(V0 - LocalV) >= vperiod/2.)
1098                   {
1099                     if (LocalV > V0)
1100                       vsens = -1;
1101                     else
1102                       vsens = 1;
1103                   }
1104                   U0 = LocalU; V0 = LocalV;
1105                   U1 = U0 + usens*uperiod;
1106                   V1 = V0 + vsens*vperiod;
1107                   Pts2d(i).SetCoord(U1,V1);
1108                   myProjIsDone = Standard_True;
1109                 }
1110               }
1111             }
1112           }
1113           if(!myProjIsDone && uperiod) {
1114             Standard_Real Uinf, Usup, Uaux;
1115             Uinf = Surf->Surface().FirstUParameter();
1116             Usup = Surf->Surface().LastUParameter();
1117             if((Usup - U0) > (U0 - Uinf))
1118               Uaux = 2*Uinf - U0 + uperiod;
1119             else
1120               Uaux = 2*Usup - U0 - uperiod;
1121             Extrema_GenLocateExtPS  locext(pntproj,
1122                                            Surf->Surface(),
1123                                            Uaux, V0, TolU, TolV);
1124             if (locext.IsDone())
1125               if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
1126               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1127                 (locext.Point()).Parameter(u,v);
1128                 if((Usup - U0) > (U0 - Uinf))
1129                   usens--;
1130                 else
1131                   usens++;
1132                 U0 = u; V0 = v;
1133                 U1 = U0 + usens*uperiod;
1134                 V1 = V0 + vsens*vperiod;
1135                 Pts2d(i).SetCoord(U1,V1);
1136                 myProjIsDone = Standard_True;
1137               }
1138           }
1139           if(!myProjIsDone && vperiod) {
1140             Standard_Real Vinf, Vsup, Vaux;
1141             Vinf = Surf->Surface().FirstVParameter();
1142             Vsup = Surf->Surface().LastVParameter();
1143             if((Vsup - V0) > (V0 - Vinf))
1144               Vaux = 2*Vinf - V0 + vperiod;
1145             else
1146               Vaux = 2*Vsup - V0 - vperiod;
1147             Extrema_GenLocateExtPS  locext(pntproj,
1148                                            Surf->Surface(),
1149                                            U0, Vaux, TolU, TolV) ;
1150             if (locext.IsDone())
1151               if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
1152               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1153                 (locext.Point()).Parameter(u,v);
1154                 if((Vsup - V0) > (V0 - Vinf))
1155                   vsens--;
1156                 else
1157                   vsens++;
1158                 U0 = u; V0 = v;
1159                 U1 = U0 + usens*uperiod;
1160                 V1 = V0 + vsens*vperiod;
1161                 Pts2d(i).SetCoord(U1,V1);
1162                 myProjIsDone = Standard_True;
1163               }
1164           }
1165           if(!myProjIsDone && uperiod && vperiod) {
1166             Standard_Real Uaux, Vaux;
1167             if((Usup - U0) > (U0 - Uinf))
1168               Uaux = 2*Uinf - U0 + uperiod;
1169             else
1170               Uaux = 2*Usup - U0 - uperiod;
1171             if((Vsup - V0) > (V0 - Vinf))
1172               Vaux = 2*Vinf - V0 + vperiod;
1173             else
1174               Vaux = 2*Vsup - V0 - vperiod;
1175             Extrema_GenLocateExtPS  locext(pntproj,
1176                                            Surf->Surface(),
1177                                            Uaux, Vaux, TolU, TolV);
1178             if (locext.IsDone())
1179               if (locext.SquareDistance() < DistTol3d * DistTol3d) {
1180               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1181                 (locext.Point()).Parameter(u,v);
1182                 if((Usup - U0) > (U0 - Uinf))
1183                   usens--;
1184                 else
1185                   usens++;
1186                 if((Vsup - V0) > (V0 - Vinf))
1187                   vsens--;
1188                 else
1189                   vsens++;
1190                 U0 = u; V0 = v;
1191                 U1 = U0 + usens*uperiod;
1192                 V1 = V0 + vsens*vperiod;
1193                 Pts2d(i).SetCoord(U1,V1);
1194                 myProjIsDone = Standard_True;
1195               }
1196           }
1197           if(!myProjIsDone) {
1198             Extrema_ExtPS ext(pntproj, Surf->Surface(), TolU, TolV) ;
1199             if (ext.IsDone()) {
1200               Dist2Min = ext.SquareDistance(1);
1201               Standard_Integer GoodValue = 1;
1202               for ( j = 2 ; j <= ext.NbExt() ; j++ )
1203                 if( Dist2Min > ext.SquareDistance(j)) {
1204                   Dist2Min = ext.SquareDistance(j);
1205                   GoodValue = j;
1206                 }
1207               if (Dist2Min < DistTol3d * DistTol3d) {
1208               //if (Dist2Min < Tol3d * Tol3d) {
1209                 (ext.Point(GoodValue)).Parameter(u,v);
1210                 if(uperiod) {
1211                   if((U0 - u) > (2*uperiod/3)) {
1212                     usens++;
1213                   }
1214                   else
1215                     if((u - U0) > (2*uperiod/3)) {
1216                       usens--;
1217                     }
1218     }
1219                 if(vperiod) {
1220                   if((V0 - v) > (vperiod/2)) {
1221                     vsens++;
1222                   }
1223                   else
1224                     if((v - V0) > (vperiod/2)) {
1225                       vsens--;
1226                     }
1227       }
1228                 U0 = u; V0 = v;
1229                 U1 = U0 + usens*uperiod;
1230                 V1 = V0 + vsens*vperiod;
1231                 Pts2d(i).SetCoord(U1,V1);
1232                 myProjIsDone = Standard_True;
1233               }
1234             }
1235           }
1236         }
1237         else break;
1238     }
1239   }
1240   // -- Pnts2d is transformed into Geom2d_BSplineCurve, with the help of Param and Mult
1241   if(myProjIsDone) {
1242     myBSpline = new Geom2d_BSplineCurve(Pts2d,Param,Mult,1);
1243     //jgv: put the curve into parametric range
1244     gp_Pnt2d MidPoint = myBSpline->Value(0.5*(myBSpline->FirstParameter() + myBSpline->LastParameter()));
1245     Standard_Real TestU = MidPoint.X(), TestV = MidPoint.Y();
1246     Standard_Real sense = 0.;
1247     if (uperiod)
1248     {
1249       if (TestU < Uinf - TolU)
1250         sense = 1.;
1251       else if (TestU > Usup + TolU)
1252         sense = -1;
1253       while (TestU < Uinf - TolU || TestU > Usup + TolU)
1254         TestU += sense * uperiod;
1255     }
1256     if (vperiod)
1257     {
1258       sense = 0.;
1259       if (TestV < Vinf - TolV)
1260         sense = 1.;
1261       else if (TestV > Vsup + TolV)
1262         sense = -1.;
1263       while (TestV < Vinf - TolV || TestV > Vsup + TolV)
1264         TestV += sense * vperiod;
1265     }
1266     gp_Vec2d Offset(TestU - MidPoint.X(), TestV - MidPoint.Y());
1267     if (Abs(Offset.X()) > gp::Resolution() ||
1268         Abs(Offset.Y()) > gp::Resolution())
1269       myBSpline->Translate(Offset);
1270     //////////////////////////////////////////
1271     Geom2dAdaptor_Curve GAC(myBSpline);
1272     Handle(Adaptor2d_HCurve2d) IC2d = new Geom2dAdaptor_HCurve(GAC);
1273 #ifdef DEB
1274 //    char name [100];
1275 //    sprintf(name,"%s_%d","build",compteur++);
1276 //    DrawTrSurf::Set(name,myBSpline);
1277 #endif
1278     return IC2d;
1279   }
1280   else {
1281 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:50 2002 Begin
1282 //     Standard_NoSuchObject_Raise_if(1,"ProjLib_Compu: build echec");
1283 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:51 2002 End
1284     return Handle(Adaptor2d_HCurve2d)();
1285   }
1286 //  myProjIsDone = Standard_False;
1287 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:01 2002 Begin
1288 //   Standard_NoSuchObject_Raise_if(1,"ProjLib_ComputeOnPS: build echec");
1289 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:02 2002 End
1290 }
1291
1292
1293
1294
1295 //=======================================================================
1296 //function : ProjLib_ProjectUsingInitialCurve2d
1297 //purpose  :
1298 //=======================================================================
1299 Handle(Geom2d_BSplineCurve)
1300      ProjLib_ComputeApproxOnPolarSurface::
1301      ProjectUsingInitialCurve2d(const Handle(Adaptor3d_HCurve)& Curve,
1302                                 const Handle(Adaptor3d_HSurface)& Surf,
1303                                 const Handle(Adaptor2d_HCurve2d)& InitCurve2d)
1304 {
1305   //OCC217
1306   Standard_Real Tol3d = myTolerance;
1307   Standard_Real DistTol3d = 1.0*Tol3d;
1308   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
1309   Standard_Real Tol2d = Sqrt(TolU*TolU + TolV*TolV);
1310
1311   Standard_Integer i;
1312   GeomAbs_SurfaceType TheTypeS = Surf->GetType();
1313   GeomAbs_CurveType TheTypeC = Curve->GetType();
1314 //  Handle(Standard_Type) TheTypeS = Surf->DynamicType();
1315 //  Handle(Standard_Type) TheTypeC = Curve->DynamicType();   // si on a :
1316 //  if(TheTypeS == STANDARD_TYPE(Geom_BSplineSurface)) {
1317   if(TheTypeS == GeomAbs_Plane) {
1318     Standard_Real S, T;
1319     gp_Pln Plane = Surf->Plane();
1320     if(TheTypeC == GeomAbs_BSplineCurve) {
1321       Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1322       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1323       for(i = 1;i <= Curve->NbPoles();i++) {
1324         ElSLib::Parameters( Plane, BSC->Pole(i), S, T);
1325         Poles2d(i).SetCoord(S,T);
1326       }
1327       TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1328       BSC->Knots(Knots);
1329       TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1330       BSC->Multiplicities(Mults);
1331       if(BSC->IsRational()) {
1332         TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1333         BSC->Weights(Weights);
1334         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1335                                        BSC->Degree(), BSC->IsPeriodic()) ;
1336       }
1337       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1338                                      BSC->Degree(), BSC->IsPeriodic()) ;
1339
1340     }
1341     if(TheTypeC == GeomAbs_BezierCurve) {
1342       Handle(Geom_BezierCurve) BC = Curve->Bezier();
1343       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1344       for(i = 1;i <= Curve->NbPoles();i++) {
1345         ElSLib::Parameters( Plane, BC->Pole(i), S, T);
1346         Poles2d(i).SetCoord(S,T);
1347       }
1348       TColStd_Array1OfReal Knots(1, 2);
1349       Knots.SetValue(1,0.0);
1350       Knots.SetValue(2,1.0);
1351       TColStd_Array1OfInteger Mults(1, 2);
1352       Mults.Init(BC->NbPoles());
1353       if(BC->IsRational()) {
1354         TColStd_Array1OfReal Weights(1, BC->NbPoles());
1355         BC->Weights(Weights);
1356         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1357                                        BC->Degree(), BC->IsPeriodic()) ;
1358       }
1359       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1360                                      BC->Degree(), BC->IsPeriodic()) ;
1361     }
1362   }
1363   if(TheTypeS == GeomAbs_BSplineSurface) {
1364     Handle(Geom_BSplineSurface) BSS = Surf->BSpline();
1365     if((BSS->MaxDegree() == 1) &&
1366        (BSS->NbUPoles() == 2) &&
1367        (BSS->NbVPoles() == 2)) {
1368       gp_Pnt p11 = BSS->Pole(1,1);
1369       gp_Pnt p12 = BSS->Pole(1,2);
1370       gp_Pnt p21 = BSS->Pole(2,1);
1371       gp_Pnt p22 = BSS->Pole(2,2);
1372       gp_Vec V1(p11,p12);
1373       gp_Vec V2(p21,p22);
1374       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){  //OCC217
1375       //if(V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){
1376         //  so the polar surface is plane
1377         //  and if it is enough to projet the  poles of Curve
1378         Standard_Integer Dist2Min = IntegerLast();
1379         Standard_Real u,v;
1380         //OCC217
1381         //Standard_Real TolU = Surf->UResolution(myTolerance)
1382         //  , TolV = Surf->VResolution(myTolerance);
1383 //      gp_Pnt pntproj;
1384         if(TheTypeC == GeomAbs_BSplineCurve) {
1385           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1386           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1387           for(i = 1;i <= Curve->NbPoles();i++) {
1388             myProjIsDone = Standard_False;
1389             Dist2Min = IntegerLast();
1390             Extrema_GenLocateExtPS  extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1391                                             (p11.Y()+p22.Y())/2,TolU,TolV) ;
1392             if (extrloc.IsDone()) {
1393               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1394               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1395               //if (Dist2Min < myTolerance * myTolerance) {
1396                 (extrloc.Point()).Parameter(u,v);
1397                 Poles2d(i).SetCoord(u,v);
1398                 myProjIsDone = Standard_True;
1399               }
1400               else break;
1401             }
1402             else break;
1403             if(!myProjIsDone)
1404               break;
1405           }
1406           if(myProjIsDone) {
1407             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1408             BSC->Knots(Knots);
1409             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1410             BSC->Multiplicities(Mults);
1411             if(BSC->IsRational()) {
1412               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1413               BSC->Weights(Weights);
1414               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1415                                              BSC->Degree(), BSC->IsPeriodic()) ;
1416             }
1417             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1418                                            BSC->Degree(), BSC->IsPeriodic()) ;
1419
1420
1421           }
1422         }
1423         if(TheTypeC == GeomAbs_BezierCurve) {
1424           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1425           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1426           for(i = 1;i <= Curve->NbPoles();i++) {
1427             Dist2Min = IntegerLast();
1428             Extrema_GenLocateExtPS  extrloc(BC->Pole(i),Surf->Surface(),0.5,
1429                                             0.5,TolU,TolV) ;
1430             if (extrloc.IsDone()) {
1431               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1432               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1433               //if (Dist2Min < myTolerance * myTolerance) {
1434                 (extrloc.Point()).Parameter(u,v);
1435                 Poles2d(i).SetCoord(u,v);
1436                 myProjIsDone = Standard_True;
1437               }
1438               else break;
1439             }
1440             else break;
1441             if(myProjIsDone)
1442               myProjIsDone = Standard_False;
1443             else break;
1444           }
1445           if(myProjIsDone) {
1446             TColStd_Array1OfReal Knots(1, 2);
1447             Knots.SetValue(1,0.0);
1448             Knots.SetValue(2,1.0);
1449             TColStd_Array1OfInteger Mults(1, 2);
1450             Mults.Init(BC->NbPoles());
1451             if(BC->IsRational()) {
1452               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1453               BC->Weights(Weights);
1454               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1455                                                     BC->Degree(), BC->IsPeriodic()) ;
1456             }
1457             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1458                                                   BC->Degree(), BC->IsPeriodic()) ;
1459           }
1460         }
1461       }
1462     }
1463   }
1464   else if(TheTypeS == GeomAbs_BezierSurface) {
1465     Handle(Geom_BezierSurface) BS = Surf->Bezier();
1466     if((BS->MaxDegree() == 1) &&
1467        (BS->NbUPoles() == 2) &&
1468        (BS->NbVPoles() == 2)) {
1469       gp_Pnt p11 = BS->Pole(1,1);
1470       gp_Pnt p12 = BS->Pole(1,2);
1471       gp_Pnt p21 = BS->Pole(2,1);
1472       gp_Pnt p22 = BS->Pole(2,2);
1473       gp_Vec V1(p11,p12);
1474       gp_Vec V2(p21,p22);
1475       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){ //OCC217
1476       //if (V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){
1477         //    and if it is enough to project the poles of Curve
1478         Standard_Integer Dist2Min = IntegerLast();
1479         Standard_Real u,v;
1480         //OCC217
1481         //Standard_Real TolU = Surf->UResolution(myTolerance)
1482         //  , TolV = Surf->VResolution(myTolerance);
1483
1484 //      gp_Pnt pntproj;
1485         if(TheTypeC == GeomAbs_BSplineCurve) {
1486           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1487           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1488           for(i = 1;i <= Curve->NbPoles();i++) {
1489             myProjIsDone = Standard_False;
1490             Dist2Min = IntegerLast();
1491             Extrema_GenLocateExtPS  extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1492                                             (p11.Y()+p22.Y())/2,TolU,TolV) ;
1493             if (extrloc.IsDone()) {
1494               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1495               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1496               //if (Dist2Min < myTolerance * myTolerance) {
1497                 (extrloc.Point()).Parameter(u,v);
1498                 Poles2d(i).SetCoord(u,v);
1499                 myProjIsDone = Standard_True;
1500               }
1501               else break;
1502             }
1503             else break;
1504             if(!myProjIsDone)
1505               break;
1506           }
1507           if(myProjIsDone) {
1508             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1509             BSC->Knots(Knots);
1510             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1511             BSC->Multiplicities(Mults);
1512             if(BSC->IsRational()) {
1513               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1514               BSC->Weights(Weights);
1515               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1516                                                     BSC->Degree(), BSC->IsPeriodic()) ;
1517             }
1518             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1519                                                   BSC->Degree(), BSC->IsPeriodic()) ;
1520
1521
1522           }
1523         }
1524         if(TheTypeC == GeomAbs_BezierCurve) {
1525           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1526           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1527           for(i = 1;i <= Curve->NbPoles();i++) {
1528             Dist2Min = IntegerLast();
1529             Extrema_GenLocateExtPS  extrloc(BC->Pole(i),Surf->Surface(),0.5,
1530                                             0.5,TolU,TolV) ;
1531             if (extrloc.IsDone()) {
1532               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1533               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1534               //if (Dist2Min < myTolerance * myTolerance) {
1535                 (extrloc.Point()).Parameter(u,v);
1536                 Poles2d(i).SetCoord(u,v);
1537                 myProjIsDone = Standard_True;
1538               }
1539               else break;
1540             }
1541             else break;
1542             if(myProjIsDone)
1543               myProjIsDone = Standard_False;
1544             else break;
1545           }
1546           if(myProjIsDone) {
1547             TColStd_Array1OfReal Knots(1, 2);
1548             Knots.SetValue(1,0.0);
1549             Knots.SetValue(2,1.0);
1550             TColStd_Array1OfInteger Mults(1, 2);
1551             Mults.Init(BC->NbPoles());
1552             if(BC->IsRational()) {
1553               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1554               BC->Weights(Weights);
1555               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1556                                                     BC->Degree(), BC->IsPeriodic()) ;
1557             }
1558             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1559                                                   BC->Degree(), BC->IsPeriodic()) ;
1560           }
1561         }
1562       }
1563     }
1564   }
1565
1566   ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, Tol3d) ;  //OCC217
1567   //ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, myTolerance) ;
1568
1569 #ifdef DEB
1570   Standard_Integer Nb = 50;
1571
1572   Standard_Real U, U1, U2;
1573   U1 = F.FirstParameter();
1574   U2 = F.LastParameter();
1575
1576   TColgp_Array1OfPnt2d    DummyPoles(1,Nb+1);
1577   TColStd_Array1OfReal    DummyKnots(1,Nb+1);
1578   TColStd_Array1OfInteger DummyMults(1,Nb+1);
1579   DummyMults.Init(1);
1580   DummyMults(1) = 2;
1581   DummyMults(Nb+1) = 2;
1582   for (Standard_Integer ij = 0; ij <= Nb; ij++) {
1583     U = (Nb-ij)*U1 + ij*U2;
1584     U /= Nb;
1585     DummyPoles(ij+1) = F.Value(U);
1586     DummyKnots(ij+1) = ij;
1587   }
1588   Handle(Geom2d_BSplineCurve) DummyC2d =
1589     new Geom2d_BSplineCurve(DummyPoles, DummyKnots, DummyMults, 1);
1590 #ifdef DRAW
1591   Standard_CString Temp = "bs2d";
1592   DrawTrSurf::Set(Temp,DummyC2d);
1593 #endif
1594 //  DrawTrSurf::Set((Standard_CString ) "bs2d",DummyC2d);
1595   Handle(Geom2dAdaptor_HCurve) DDD =
1596     Handle(Geom2dAdaptor_HCurve)::DownCast(InitCurve2d);
1597
1598 #ifdef DRAW
1599   Temp = "initc2d";
1600   DrawTrSurf::Set(Temp,DDD->ChangeCurve2d().Curve());
1601 #endif
1602 //  DrawTrSurf::Set((Standard_CString ) "initc2d",DDD->ChangeCurve2d().Curve());
1603 #endif
1604
1605   Standard_Integer Deg1,Deg2;
1606 //  Deg1 = 8;
1607 //  Deg2 = 8;
1608   Deg1 = 2; //IFV
1609   Deg2 = 8; //IFV
1610
1611   Approx_FitAndDivide2d Fit(F,Deg1,Deg2,Tol3d,Tol2d, //OCC217
1612   //Approx_FitAndDivide2d Fit(F,Deg1,Deg2,myTolerance,myTolerance,
1613                             Standard_True);
1614
1615   if(Fit.IsAllApproximated()) {
1616     Standard_Integer i;
1617     Standard_Integer NbCurves = Fit.NbMultiCurves();
1618     Standard_Integer MaxDeg = 0;
1619     // To transform the MultiCurve into BSpline, it is required that all
1620     // Bezier constituing it have the same degree -> Calculation of MaxDeg
1621     Standard_Integer NbPoles  = 1;
1622     for (i = 1; i <= NbCurves; i++) {
1623       Standard_Integer Deg = Fit.Value(i).Degree();
1624       MaxDeg = Max ( MaxDeg, Deg);
1625     }
1626
1627     NbPoles = MaxDeg * NbCurves + 1;               //Tops on the BSpline
1628     TColgp_Array1OfPnt2d  Poles( 1, NbPoles);
1629
1630     TColgp_Array1OfPnt2d TempPoles( 1, MaxDeg + 1);//to augment the degree
1631
1632     TColStd_Array1OfReal Knots( 1, NbCurves + 1);  //Nodes of the BSpline
1633
1634     Standard_Integer Compt = 1;
1635     for (i = 1; i <= NbCurves; i++) {
1636       Fit.Parameters(i, Knots(i), Knots(i+1));
1637       AppParCurves_MultiCurve MC = Fit.Value( i);       //Load the Ith Curve
1638       TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Retrieve the tops
1639       MC.Curve(1, Poles2d);
1640
1641       //Eventual augmentation of the degree
1642       Standard_Integer Inc = MaxDeg - MC.Degree();
1643       if ( Inc > 0) {
1644 //      BSplCLib::IncreaseDegree( Inc, Poles2d, PLib::NoWeights(),
1645         BSplCLib::IncreaseDegree( MaxDeg, Poles2d, PLib::NoWeights(),
1646                          TempPoles, PLib::NoWeights());
1647         //update of tops of the PCurve
1648         for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
1649           Poles.SetValue( Compt, TempPoles( j));
1650           Compt++;
1651         }
1652       }
1653       else {
1654         //update of tops of the PCurve
1655         for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
1656           Poles.SetValue( Compt, Poles2d( j));
1657           Compt++;
1658         }
1659       }
1660
1661       Compt--;
1662     }
1663
1664     //update of fields of  ProjLib_Approx
1665     Standard_Integer NbKnots = NbCurves + 1;
1666
1667     // The start and end nodes are not correct : Cf: opening of the interval
1668     //Knots( 1) -= 1.e-9;
1669     //Knots(NbKnots) += 1.e-9;
1670
1671
1672     TColStd_Array1OfInteger   Mults( 1, NbKnots);
1673     Mults.Init(MaxDeg);
1674     Mults.SetValue( 1, MaxDeg + 1);
1675     Mults.SetValue(NbKnots, MaxDeg + 1);
1676     myProjIsDone = Standard_True;
1677     Handle(Geom2d_BSplineCurve) Dummy =
1678       new Geom2d_BSplineCurve(Poles,Knots,Mults,MaxDeg);
1679
1680     // try to smoother the Curve GeomAbs_C1.
1681
1682     Standard_Boolean OK = Standard_True;
1683
1684     for (Standard_Integer ij = 2; ij < NbKnots; ij++) {
1685       OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,Tol3d);  //OCC217
1686       //OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,myTolerance);
1687     }
1688 #ifdef DEB
1689     if (!OK) {
1690       cout << "ProjLib_ComputeApproxOnPolarSurface : Smoothing echoue"<<endl;
1691     }
1692 #endif
1693     return Dummy;
1694   }
1695   return Handle(Geom2d_BSplineCurve)();
1696 }
1697
1698 //=======================================================================
1699 //function : BSpline
1700 //purpose  :
1701 //=======================================================================
1702
1703 Handle(Geom2d_BSplineCurve)
1704      ProjLib_ComputeApproxOnPolarSurface::BSpline() const
1705
1706 {
1707 //  Modified by Sergey KHROMOV - Thu Apr 18 11:16:46 2002 End
1708 //   Standard_NoSuchObject_Raise_if
1709 //     (!myProjIsDone,
1710 //      "ProjLib_ComputeApproxOnPolarSurface:BSpline");
1711 //  Modified by Sergey KHROMOV - Thu Apr 18 11:16:47 2002 End
1712   return myBSpline ;
1713 }
1714
1715 //=======================================================================
1716 //function : Curve2d
1717 //purpose  :
1718 //=======================================================================
1719
1720 Handle(Geom2d_Curve)
1721      ProjLib_ComputeApproxOnPolarSurface::Curve2d() const
1722
1723 {
1724   Standard_NoSuchObject_Raise_if
1725     (!myProjIsDone,
1726      "ProjLib_ComputeApproxOnPolarSurface:2ndCurve2d");
1727   return my2ndCurve ;
1728 }
1729
1730
1731 //=======================================================================
1732 //function : IsDone
1733 //purpose  :
1734 //=======================================================================
1735
1736 Standard_Boolean ProjLib_ComputeApproxOnPolarSurface::IsDone() const
1737
1738 {
1739   return myProjIsDone;
1740 }