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_Function.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 OCCT_DEBUG
  78 #ifdef DRAW
  79 #include <DrawTrSurf.hxx>
  80 #include <Geom2d_Curve.hxx>
  81 #endif
  82 //static Standard_Integer compteur = 0;
  83 #endif
  84
  85 struct aFuncStruct
  86 {
  87   Handle(Adaptor3d_HSurface) mySurf; // Surface where to project.
  88   Handle(Adaptor3d_HCurve)   myCurve; // Curve to project.
  89   Handle(Adaptor2d_HCurve2d) myInitCurve2d; // Initial 2dcurve projection.
  90   Standard_Real mySqProjOrtTol; // Used to filter non-orthogonal projected point.
  91   Standard_Real myTolU;
  92   Standard_Real myTolV;
  93   Standard_Real myPeriod[2]; // U and V period correspondingly.
  94 };
  95
  96 //=======================================================================
  97 //function : computePeriodicity
  98 //purpose  : Compute period information on adaptor.
  99 //=======================================================================
 100 static void computePeriodicity(const Handle(Adaptor3d_HSurface)& theSurf,
 101                                Standard_Real &theUPeriod,
 102                                Standard_Real &theVPeriod)
 103 {
 104   theUPeriod = 0.0;
 105   theVPeriod = 0.0;
 106
 107   // Compute once information about periodicity.
 108   // Param space may be reduced in case of rectangular trimmed surface,
 109   // in this case really trimmed bounds should be set as unperiodic.
 110   Standard_Real aTrimF, aTrimL, aBaseF, aBaseL, aDummyF, aDummyL;
 111   Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theSurf->Surface(), Standard_False); // Not trim.
 112   // U param space.
 113   if (theSurf->IsUPeriodic())
 114   {
 115     theUPeriod = theSurf->UPeriod();
 116   }
 117   else if(theSurf->IsUClosed())
 118   {
 119     theUPeriod = theSurf->LastUParameter() - theSurf->FirstUParameter();
 120   }
 121   if (theUPeriod != 0.0)
 122   {
 123     aTrimF = theSurf->FirstUParameter(); // Trimmed first
 124     aTrimL = theSurf->LastUParameter(); // Trimmed last
 125     aS->Bounds(aBaseF, aBaseL, aDummyF, aDummyL); // Non-trimmed values.
 126     if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
 127     {
 128       // Param space reduced.
 129       theUPeriod = 0.0;
 130     }
 131   }
 132
 133   // V param space.
 134   if (theSurf->IsVPeriodic())
 135   {
 136     theVPeriod = theSurf->VPeriod();
 137   }
 138   else if(theSurf->IsVClosed())
 139   {
 140     theVPeriod = theSurf->LastVParameter() - theSurf->FirstVParameter();
 141   }
 142   if (theVPeriod != 0.0)
 143   {
 144     aTrimF = theSurf->FirstVParameter(); // Trimmed first
 145     aTrimL = theSurf->LastVParameter(); // Trimmed last
 146     aS->Bounds(aDummyF, aDummyL, aBaseF, aBaseL); // Non-trimmed values.
 147     if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
 148     {
 149       // Param space reduced.
 150       theVPeriod = 0.0;
 151     }
 152   }
 153 }
 154
 155 //=======================================================================
 156 //function : aFuncValue
 157 //purpose  : compute functional value in (theU,theV) point
 158 //=======================================================================
 159 static Standard_Real anOrthogSqValue(const gp_Pnt& aBasePnt,
 160                                      const Handle(Adaptor3d_HSurface)& Surf,
 161                                      const Standard_Real theU,
 162                                      const Standard_Real theV)
 163 {
 164   // Since find projection, formula is:
 165   // F1 = Dot(S_U, Vec(aBasePnt, aProjPnt))
 166   // F2 = Dot(S_V, Vec(aBasePnt, aProjPnt))
 167
 168   gp_Pnt aProjPnt;
 169   gp_Vec aSu, aSv;
 170
 171   Surf->D1(theU, theV, aProjPnt, aSu, aSv);
 172   gp_Vec aBaseVec(aBasePnt, aProjPnt);
 173
 174   if (aSu.SquareMagnitude() > Precision::SquareConfusion())
 175     aSu.Normalize();
 176
 177   if (aSv.SquareMagnitude() > Precision::SquareConfusion())
 178     aSv.Normalize();
 179
 180   Standard_Real aFirstPart = aSu.Dot(aBaseVec);
 181   Standard_Real aSecondPart = aSv.Dot(aBaseVec);
 182   return (aFirstPart * aFirstPart + aSecondPart * aSecondPart);
 183 }
 184
 185 //=======================================================================
 186 //function : Value
 187 //purpose  : (OCC217 - apo)- Compute Point2d that project on polar surface(<Surf>) 3D<Curve>
 188 //            <InitCurve2d> use for calculate start 2D point.
 189 //=======================================================================
 190 static gp_Pnt2d Function_Value(const Standard_Real theU,
 191                                const aFuncStruct& theData)
 192 {
 193   gp_Pnt2d p2d = theData.myInitCurve2d->Value(theU) ;
 194   gp_Pnt p = theData.myCurve->Value(theU);
 195   gp_Pnt aSurfPnt = theData.mySurf->Value(p2d.X(), p2d.Y());
 196   Standard_Real aSurfPntDist = aSurfPnt.SquareDistance(p);
 197
 198   Standard_Real Uinf, Usup, Vinf, Vsup;
 199   Uinf = theData.mySurf->Surface().FirstUParameter();
 200   Usup = theData.mySurf->Surface().LastUParameter();
 201   Vinf = theData.mySurf->Surface().FirstVParameter();
 202   Vsup = theData.mySurf->Surface().LastVParameter();
 203
 204   // Check case when curve is close to co-parametrized isoline on surf.
 205   if (Abs (p2d.X() - Uinf) < Precision::PConfusion() ||
 206       Abs (p2d.X() - Usup) < Precision::PConfusion() )
 207   {
 208     // V isoline.
 209     gp_Pnt aPnt;
 210     theData.mySurf->D0(p2d.X(), theU, aPnt);
 211     if (aPnt.SquareDistance(p) < aSurfPntDist)
 212       p2d.SetY(theU);
 213   }
 214
 215   if (Abs (p2d.Y() - Vinf) < Precision::PConfusion() ||
 216       Abs (p2d.Y() - Vsup) < Precision::PConfusion() )
 217   {
 218     // U isoline.
 219     gp_Pnt aPnt;
 220     theData.mySurf->D0(theU, p2d.Y(), aPnt);
 221     if (aPnt.SquareDistance(p) < aSurfPntDist)
 222       p2d.SetX(theU);
 223   }
 224
 225   Standard_Integer decalU = 0, decalV = 0;
 226   Standard_Real U0 = p2d.X(), V0 = p2d.Y();
 227
 228   GeomAbs_SurfaceType Type = theData.mySurf->GetType();
 229   if((Type != GeomAbs_BSplineSurface) &&
 230      (Type != GeomAbs_BezierSurface)  &&
 231      (Type != GeomAbs_OffsetSurface)    )
 232   {
 233     // Analytical cases.
 234     Standard_Real S = 0., T = 0.;
 235     switch (Type)
 236     {
 237     case GeomAbs_Cylinder:
 238       {
 239         gp_Cylinder Cylinder = theData.mySurf->Cylinder();
 240         ElSLib::Parameters( Cylinder, p, S, T);
 241         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 242         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 243         S += decalU*2*M_PI;
 244         break;
 245       }
 246     case GeomAbs_Cone:
 247       {
 248         gp_Cone Cone = theData.mySurf->Cone();
 249         ElSLib::Parameters( Cone, p, S, T);
 250         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 251         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 252         S += decalU*2*M_PI;
 253         break;
 254       }
 255     case GeomAbs_Sphere:
 256       {
 257         gp_Sphere Sphere = theData.mySurf->Sphere();
 258         ElSLib::Parameters( Sphere, p, S, T);
 259         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 260         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 261         S += decalU*2*M_PI;
 262         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
 263         if(V0 > (Vsup+(Vsup-Vinf))) decalV =  int((V0 - Vsup+(Vsup-Vinf))/(2*M_PI))+1;
 264         T += decalV*2*M_PI;
 265         if(0.4*M_PI < Abs(U0 - S) && Abs(U0 - S) < 1.6*M_PI)
 266         {
 267           T = M_PI - T;
 268           if(U0 < S)
 269             S -= M_PI;
 270           else
 271             S += M_PI;
 272         }
 273         break;
 274       }
 275     case GeomAbs_Torus:
 276       {
 277         gp_Torus Torus = theData.mySurf->Torus();
 278         ElSLib::Parameters( Torus, p, S, T);
 279         if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
 280         if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
 281         if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
 282         if(V0 > Vsup) decalV =  int((V0 - Vsup)/(2*M_PI))+1;
 283         S += decalU*2*M_PI; T += decalV*2*M_PI;
 284         break;
 285       }
 286     default:
 287       throw Standard_NoSuchObject("ProjLib_ComputeApproxOnPolarSurface::Value");
 288     }
 289     return gp_Pnt2d(S, T);
 290   }
 291
 292   // Non-analytical case.
 293   Standard_Real Dist2Min = RealLast();
 294   Standard_Real uperiod = theData.myPeriod[0],
 295                 vperiod = theData.myPeriod[1],
 296                 u, v;
 297
 298   // U0 and V0 are the points within the initialized period.
 299   if(U0 < Uinf)
 300   {
 301     if(!uperiod)
 302       U0 = Uinf;
 303     else
 304     {
 305       decalU = int((Uinf - U0)/uperiod)+1;
 306       U0 += decalU*uperiod;
 307     }
 308   }
 309   if(U0 > Usup)
 310   {
 311     if(!uperiod)
 312       U0 = Usup;
 313     else
 314     {
 315       decalU = -(int((U0 - Usup)/uperiod)+1);
 316       U0 += decalU*uperiod;
 317     }
 318   }
 319   if(V0 < Vinf)
 320   {
 321     if(!vperiod)
 322       V0 = Vinf;
 323     else
 324     {
 325       decalV = int((Vinf - V0)/vperiod)+1;
 326       V0 += decalV*vperiod;
 327     }
 328   }
 329   if(V0 > Vsup)
 330   {
 331     if(!vperiod)
 332       V0 = Vsup;
 333     else
 334     {
 335       decalV = -int((V0 - Vsup)/vperiod)-1;
 336       V0 += decalV*vperiod;
 337     }
 338   }
 339
 340   // The surface around (U0,V0) is reduced.
 341   Standard_Real uLittle = (Usup - Uinf)/10, vLittle = (Vsup - Vinf)/10;
 342   Standard_Real uInfLi = 0, vInfLi = 0,uSupLi = 0, vSupLi = 0;
 343   if((U0 - Uinf) > uLittle) uInfLi = U0 - uLittle; else uInfLi = Uinf;
 344   if((V0 - Vinf) > vLittle) vInfLi = V0 - vLittle; else vInfLi = Vinf;
 345   if((Usup - U0) > uLittle) uSupLi = U0 + uLittle; else uSupLi = Usup;
 346   if((Vsup - V0) > vLittle) vSupLi = V0 + vLittle; else vSupLi = Vsup;
 347
 348   GeomAdaptor_Surface SurfLittle;
 349   if (Type == GeomAbs_BSplineSurface)
 350   {
 351     Handle(Geom_Surface) GBSS(theData.mySurf->Surface().BSpline());
 352     SurfLittle.Load(GBSS, uInfLi, uSupLi, vInfLi, vSupLi);
 353   }
 354   else if (Type == GeomAbs_BezierSurface)
 355   {
 356     Handle(Geom_Surface) GS(theData.mySurf->Surface().Bezier());
 357     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
 358   }
 359   else if (Type == GeomAbs_OffsetSurface)
 360   {
 361     Handle(Geom_Surface) GS = GeomAdaptor::MakeSurface(theData.mySurf->Surface());
 362     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
 363   }
 364   else
 365   {
 366     throw Standard_NoSuchObject ("ProjLib_ComputeApproxOnPolarSurface::ProjectUsingInitialCurve2d() - unknown surface type");
 367   }
 368
 369   // Try to run simple search with initial point (U0, V0).
 370   Extrema_GenLocateExtPS  locext(SurfLittle, theData.myTolU, theData.myTolV);
 371   locext.Perform(p, U0, V0);
 372   if (locext.IsDone())
 373   {
 374     locext.Point().Parameter(u, v);
 375     Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
 376     if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
 377         locext.SquareDistance() < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
 378     {
 379       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
 380       return pnt;
 381     }
 382   }
 383
 384   // Perform whole param space search.
 385   Extrema_ExtPS  ext(p, SurfLittle, theData.myTolU, theData.myTolV);
 386   if (ext.IsDone() && ext.NbExt() >= 1)
 387   {
 388     Dist2Min = ext.SquareDistance(1);
 389     Standard_Integer GoodValue = 1;
 390     for (Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ )
 391     {
 392       if( Dist2Min > ext.SquareDistance(i))
 393       {
 394         Dist2Min = ext.SquareDistance(i);
 395         GoodValue = i;
 396       }
 397     }
 398     ext.Point(GoodValue).Parameter(u, v);
 399     Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
 400     if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
 401         ext.SquareDistance(GoodValue) < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
 402     {
 403       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
 404       return pnt;
 405     }
 406   }
 407
 408   // Both searches return bad values, use point from initial 2dcurve.
 409   return p2d;
 410 }
 411
 412
 413 //=======================================================================
 414 //function : ProjLib_PolarFunction
 415 //purpose  : (OCC217 - apo)- This class produce interface to call "gp_Pnt2d Function_Value(...)"
 416 //=======================================================================
 417
 418 class ProjLib_PolarFunction : public AppCont_Function
 419 {
 420   aFuncStruct myStruct;
 421
 422   public :
 423
 424   ProjLib_PolarFunction(const Handle(Adaptor3d_HCurve) & C,
 425                         const Handle(Adaptor3d_HSurface)& Surf,
 426                         const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 427                         const Standard_Real Tol3d)
 428   {
 429     myNbPnt = 0;
 430     myNbPnt2d = 1;
 431
 432     computePeriodicity(Surf, myStruct.myPeriod[0], myStruct.myPeriod[1]);
 433
 434     myStruct.myCurve = C;
 435     myStruct.myInitCurve2d = InitialCurve2d;
 436     myStruct.mySurf = Surf;
 437     myStruct.mySqProjOrtTol = 10000.0 * Tol3d * Tol3d;
 438     myStruct.myTolU = Surf->UResolution(Tol3d);
 439     myStruct.myTolV = Surf->VResolution(Tol3d);
 440   }
 441
 442   ~ProjLib_PolarFunction() {}
 443
 444   Standard_Real FirstParameter() const
 445   {
 446     return myStruct.myCurve->FirstParameter();
 447   }
 448
 449   Standard_Real LastParameter() const
 450   {
 451     return myStruct.myCurve->LastParameter();
 452   }
 453
 454   gp_Pnt2d Value(const Standard_Real t) const
 455   {
 456     return Function_Value(t, myStruct);
 457   }
 458
 459   Standard_Boolean Value(const Standard_Real   theT,
 460                          NCollection_Array1<gp_Pnt2d>& thePnt2d,
 461                          NCollection_Array1<gp_Pnt>&   /*thePnt*/) const
 462   {
 463     thePnt2d(1) = Function_Value(theT, myStruct);
 464     return Standard_True;
 465   }
 466
 467   Standard_Boolean D1(const Standard_Real   /*theT*/,
 468                       NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
 469                       NCollection_Array1<gp_Vec>&   /*theVec*/) const
 470     {return Standard_False;}
 471 };
 472
 473 //=======================================================================
 474 //function : ProjLib_ComputeApproxOnPolarSurface
 475 //purpose  :
 476 //=======================================================================
 477
 478 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface()
 479 : myProjIsDone(Standard_False),
 480   myTolerance(Precision::Approximation()),
 481   myTolReached(-1.0),
 482   myDegMin(-1), myDegMax(-1),
 483   myMaxSegments(-1),
 484   myMaxDist(-1.),
 485   myBndPnt(AppParCurves_TangencyPoint),
 486   myDist(0.)
 487 {
 488 }
 489
 490 //=======================================================================
 491 //function : ProjLib_ComputeApproxOnPolarSurface
 492 //purpose  :
 493 //=======================================================================
 494
 495 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 496                     (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
 497                      const Handle(Adaptor3d_HCurve)&   theCurve,
 498                      const Handle(Adaptor3d_HSurface)& theSurface,
 499                      const Standard_Real               theTolerance3D)
 500 : myProjIsDone(Standard_False),
 501   myTolerance(theTolerance3D),
 502   myTolReached(-1.0),
 503   myDegMin(-1), myDegMax(-1),
 504   myMaxSegments(-1),
 505   myMaxDist(-1.),
 506   myBndPnt(AppParCurves_TangencyPoint),
 507   myDist(0.)
 508 {
 509   myBSpline = Perform(theInitialCurve2d, theCurve, theSurface);
 510 }
 511
 512 //=======================================================================
 513 //function : ProjLib_ComputeApproxOnPolarSurface
 514 //purpose  : case without curve of initialization
 515 //=======================================================================
 516
 517 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 518                       (const Handle(Adaptor3d_HCurve)&   theCurve,
 519                        const Handle(Adaptor3d_HSurface)& theSurface,
 520                        const Standard_Real               theTolerance3D)
 521 : myProjIsDone(Standard_False),
 522   myTolerance(theTolerance3D),
 523   myTolReached(-1.0),
 524   myDegMin(-1), myDegMax(-1),
 525   myMaxSegments(-1),
 526   myMaxDist(-1.),
 527   myBndPnt(AppParCurves_TangencyPoint),
 528   myDist(0.)
 529 {
 530   const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
 531   myBSpline = Perform(anInitCurve2d, theCurve, theSurface);
 532 }
 533
 534 //=======================================================================
 535 //function : ProjLib_ComputeApproxOnPolarSurface
 536 //purpose  : Process the case of sewing
 537 //=======================================================================
 538
 539 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 540                 (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
 541                  const Handle(Adaptor2d_HCurve2d)& theInitialCurve2dBis,
 542                  const Handle(Adaptor3d_HCurve)&   theCurve,
 543                  const Handle(Adaptor3d_HSurface)& theSurface,
 544                  const Standard_Real               theTolerance3D)
 545 : myProjIsDone(Standard_False),
 546   myTolerance(theTolerance3D),
 547   myTolReached(-1.0),
 548   myDegMin(-1), myDegMax(-1),
 549   myMaxSegments(-1),
 550   myMaxDist(-1.),
 551   myBndPnt(AppParCurves_TangencyPoint),
 552   myDist(0.)
 553 {
 554   // InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing
 555   Handle(Geom2d_BSplineCurve) bsc =
 556     Perform(theInitialCurve2d, theCurve, theSurface);
 557
 558   if(myProjIsDone) {
 559     gp_Pnt2d P2dproj, P2d, P2dBis;
 560     P2dproj = bsc->StartPoint();
 561     P2d    = theInitialCurve2d->Value(theInitialCurve2d->FirstParameter());
 562     P2dBis = theInitialCurve2dBis->Value(theInitialCurve2dBis->FirstParameter());
 563
 564     Standard_Real Dist, DistBis;
 565     Dist    = P2dproj.Distance(P2d);
 566     DistBis = P2dproj.Distance(P2dBis);
 567     if( Dist < DistBis)  {
 568       // myBSpline2d is the pcurve that is found. It is translated to obtain myCurve2d
 569       myBSpline = bsc;
 570       Handle(Geom2d_Geometry) GG = myBSpline->Translated(P2d, P2dBis);
 571       my2ndCurve = Handle(Geom2d_Curve)::DownCast(GG);
 572     }
 573     else {
 574       my2ndCurve = bsc;
 575       Handle(Geom2d_Geometry) GG = my2ndCurve->Translated(P2dBis, P2d);
 576       myBSpline = Handle(Geom2d_BSplineCurve)::DownCast(GG);
 577     }
 578   }
 579 }
 580
 581
 582 //=======================================================================
 583 //function : Concat
 584 //purpose  :
 585 //=======================================================================
 586
 587 static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
 588                                           Handle(Geom2d_BSplineCurve) C2,
 589                                           Standard_Real theUJump,
 590                                           Standard_Real theVJump)
 591 {
 592   Standard_Integer deg, deg1, deg2;
 593   deg1 = C1->Degree();
 594   deg2 = C2->Degree();
 595
 596   if ( deg1 < deg2) {
 597     C1->IncreaseDegree(deg2);
 598     deg = deg2;
 599   }
 600   else if ( deg2 < deg1) {
 601     C2->IncreaseDegree(deg1);
 602     deg = deg1;
 603   }
 604   else deg = deg1;
 605
 606   Standard_Integer np1,np2,nk1,nk2,np,nk;
 607   np1 = C1->NbPoles();
 608   nk1 = C1->NbKnots();
 609   np2 = C2->NbPoles();
 610   nk2 = C2->NbKnots();
 611   nk = nk1 + nk2 -1;
 612   np = np1 + np2 -1;
 613
 614   TColStd_Array1OfReal    K1(1,nk1); C1->Knots(K1);
 615   TColStd_Array1OfInteger M1(1,nk1); C1->Multiplicities(M1);
 616   TColgp_Array1OfPnt2d    P1(1,np1); C1->Poles(P1);
 617   TColStd_Array1OfReal    K2(1,nk2); C2->Knots(K2);
 618   TColStd_Array1OfInteger M2(1,nk2); C2->Multiplicities(M2);
 619   TColgp_Array1OfPnt2d    P2(1,np2); C2->Poles(P2);
 620
 621   // Compute the new BSplineCurve
 622   TColStd_Array1OfReal    K(1,nk);
 623   TColStd_Array1OfInteger M(1,nk);
 624   TColgp_Array1OfPnt2d    P(1,np);
 625
 626   Standard_Integer i, count = 0;
 627   // Set Knots and Mults
 628   for ( i = 1; i <= nk1; i++) {
 629     count++;
 630     K(count) = K1(i);
 631     M(count) = M1(i);
 632   }
 633   M(count) = deg;
 634   for ( i = 2; i <= nk2; i++) {
 635     count++;
 636     K(count) = K2(i);
 637     M(count) = M2(i);
 638   }
 639   // Set the Poles
 640   count = 0;
 641   for (i = 1; i <= np1; i++) {
 642     count++;
 643     P(count) = P1(i);
 644   }
 645   for (i = 2; i <= np2; i++) {
 646     count++;
 647     P(count).SetX(P2(i).X() + theUJump);
 648     P(count).SetY(P2(i).Y() + theVJump);
 649   }
 650
 651   Handle(Geom2d_BSplineCurve) BS =
 652     new Geom2d_BSplineCurve(P,K,M,deg);
 653   return BS;
 654 }
 655
 656 //=======================================================================
 657 //function : Perform
 658 //purpose  :
 659 //=======================================================================
 660
 661 void ProjLib_ComputeApproxOnPolarSurface::Perform
 662 (const Handle(Adaptor3d_HCurve)& Curve, const Handle(Adaptor3d_HSurface)& S)
 663 {
 664   const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
 665   myBSpline = Perform(anInitCurve2d, Curve, S);
 666 }
 667
 668 //=======================================================================
 669 //function : Perform
 670 //purpose  :
 671 //=======================================================================
 672
 673 Handle(Geom2d_BSplineCurve) ProjLib_ComputeApproxOnPolarSurface::Perform
 674 (const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 675  const Handle(Adaptor3d_HCurve)& Curve,
 676  const Handle(Adaptor3d_HSurface)& S)
 677 {
 678   //OCC217
 679   Standard_Real Tol3d = myTolerance;
 680   Standard_Real ParamTol = Precision::PApproximation();
 681
 682   Handle(Adaptor2d_HCurve2d) AHC2d = InitialCurve2d;
 683   Handle(Adaptor3d_HCurve) AHC = Curve;
 684
 685 // if the curve 3d is a BSpline with degree C0, it is cut into sections with degree C1
 686 // -> bug cts18237
 687   GeomAbs_CurveType typeCurve = Curve->GetType();
 688   if(typeCurve == GeomAbs_BSplineCurve) {
 689     TColStd_ListOfTransient LOfBSpline2d;
 690     Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
 691     Standard_Integer nbInter = Curve->NbIntervals(GeomAbs_C1);
 692     if(nbInter > 1) {
 693       Standard_Integer i, j;
 694       Handle(Geom_TrimmedCurve) GTC;
 695       Handle(Geom2d_TrimmedCurve) G2dTC;
 696       TColStd_Array1OfReal Inter(1,nbInter+1);
 697       Curve->Intervals(Inter,GeomAbs_C1);
 698       Standard_Real firstinter = Inter.Value(1), secondinter = Inter.Value(2);
 699       // initialization 3d
 700       GTC = new Geom_TrimmedCurve(BSC, firstinter, secondinter);
 701       AHC = new GeomAdaptor_HCurve(GTC);
 702
 703       // if there is an initialization curve:
 704       // - either this is a BSpline C0, with discontinuity at the same parameters of nodes
 705       // and the sections C1 are taken
 706       // - or this is a curve C1 and the sections of intrest are taken otherwise the curve is created.
 707
 708       // initialization 2d
 709       Standard_Integer nbInter2d;
 710       Standard_Boolean C2dIsToCompute;
 711       C2dIsToCompute = InitialCurve2d.IsNull();
 712       Handle(Geom2d_BSplineCurve) BSC2d;
 713       Handle(Geom2d_Curve) G2dC;
 714
 715       if(!C2dIsToCompute) {
 716         nbInter2d = InitialCurve2d->NbIntervals(GeomAbs_C1);
 717         TColStd_Array1OfReal Inter2d(1,nbInter2d+1);
 718         InitialCurve2d->Intervals(Inter2d,GeomAbs_C1);
 719         j = 1;
 720         for(i = 1,j = 1;i <= nbInter;i++)
 721           if(Abs(Inter.Value(i) - Inter2d.Value(j)) < ParamTol) { //OCC217
 722           //if(Abs(Inter.Value(i) - Inter2d.Value(j)) < myTolerance) {
 723             if (j > nbInter2d) break;
 724             j++;
 725           }
 726         if(j != (nbInter2d+1)) {
 727           C2dIsToCompute = Standard_True;
 728         }
 729       }
 730
 731       if(C2dIsToCompute) {
 732         AHC2d = BuildInitialCurve2d(AHC, S);
 733       }
 734       else {
 735         typeCurve = InitialCurve2d->GetType();
 736         switch (typeCurve) {
 737         case GeomAbs_Line: {
 738           G2dC = new Geom2d_Line(InitialCurve2d->Line());
 739           break;
 740         }
 741         case GeomAbs_Circle: {
 742           G2dC = new Geom2d_Circle(InitialCurve2d->Circle());
 743           break;
 744         }
 745         case GeomAbs_Ellipse: {
 746           G2dC = new Geom2d_Ellipse(InitialCurve2d->Ellipse());
 747           break;
 748         }
 749         case GeomAbs_Hyperbola: {
 750           G2dC = new Geom2d_Hyperbola(InitialCurve2d->Hyperbola());
 751           break;
 752         }
 753         case GeomAbs_Parabola: {
 754           G2dC = new Geom2d_Parabola(InitialCurve2d->Parabola());
 755           break;
 756         }
 757         case GeomAbs_BezierCurve: {
 758           G2dC = InitialCurve2d->Bezier();
 759           break;
 760         }
 761         case GeomAbs_BSplineCurve: {
 762           G2dC = InitialCurve2d->BSpline();
 763           break;
 764         }
 765         case GeomAbs_OtherCurve:
 766         default:
 767           break;
 768         }
 769         gp_Pnt2d fp2d = G2dC->Value(firstinter), lp2d = G2dC->Value(secondinter);
 770         gp_Pnt fps, lps, fpc, lpc;
 771         S->D0(fp2d.X(), fp2d.Y(), fps);
 772         S->D0(lp2d.X(), lp2d.Y(), lps);
 773         Curve->D0(firstinter, fpc);
 774         Curve->D0(secondinter, lpc);
 775         //OCC217
 776         if((fps.IsEqual(fpc, Tol3d)) &&
 777            (lps.IsEqual(lpc, Tol3d))) {
 778         //if((fps.IsEqual(fpc, myTolerance)) &&
 779         //   (lps.IsEqual(lpc, myTolerance))) {
 780           G2dTC = new Geom2d_TrimmedCurve(G2dC, firstinter, secondinter);
 781           Geom2dAdaptor_Curve G2dAC(G2dTC);
 782           AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
 783           myProjIsDone = Standard_True;
 784         }
 785         else {
 786           AHC2d = BuildInitialCurve2d(AHC, S);
 787           C2dIsToCompute = Standard_True;
 788         }
 789       }
 790
 791       if(myProjIsDone) {
 792         BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 793         if(BSC2d.IsNull())
 794         {
 795             return Handle(Geom2d_BSplineCurve)();
 796         }
 797         LOfBSpline2d.Append(BSC2d);
 798       }
 799       else {
 800         return Handle(Geom2d_BSplineCurve)();
 801       }
 802
 803
 804
 805       Standard_Real iinter, ip1inter;
 806       Standard_Integer nbK2d, deg;
 807       nbK2d = BSC2d->NbKnots(); deg = BSC2d->Degree();
 808
 809       for(i = 2;i <= nbInter;i++) {
 810         iinter = Inter.Value(i);
 811         ip1inter = Inter.Value(i+1);
 812         // general case 3d
 813         GTC->SetTrim(iinter, ip1inter);
 814         AHC = new GeomAdaptor_HCurve(GTC);
 815
 816         // general case 2d
 817         if(C2dIsToCompute) {
 818           AHC2d = BuildInitialCurve2d(AHC, S);
 819         }
 820         else {
 821           gp_Pnt2d fp2d = G2dC->Value(iinter), lp2d = G2dC->Value(ip1inter);
 822           gp_Pnt fps, lps, fpc, lpc;
 823           S->D0(fp2d.X(), fp2d.Y(), fps);
 824           S->D0(lp2d.X(), lp2d.Y(), lps);
 825           Curve->D0(iinter, fpc);
 826           Curve->D0(ip1inter, lpc);
 827           //OCC217
 828           if((fps.IsEqual(fpc, Tol3d)) &&
 829              (lps.IsEqual(lpc, Tol3d))) {
 830           //if((fps.IsEqual(fpc, myTolerance)) &&
 831           //   (lps.IsEqual(lpc, myTolerance))) {
 832             G2dTC->SetTrim(iinter, ip1inter);
 833             Geom2dAdaptor_Curve G2dAC(G2dTC);
 834             AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
 835             myProjIsDone = Standard_True;
 836           }
 837           else {
 838             AHC2d = BuildInitialCurve2d(AHC, S);
 839           }
 840         }
 841         if(myProjIsDone) {
 842           BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 843           if(BSC2d.IsNull()) {
 844             return Handle(Geom2d_BSplineCurve)();
 845           }
 846           LOfBSpline2d.Append(BSC2d);
 847           nbK2d += BSC2d->NbKnots() - 1;
 848           deg = Max(deg, BSC2d->Degree());
 849         }
 850         else {
 851           return Handle(Geom2d_BSplineCurve)();
 852         }
 853       }
 854
 855       Standard_Real anUPeriod, anVPeriod;
 856       computePeriodicity(S, anUPeriod, anVPeriod);
 857       Standard_Integer NbC = LOfBSpline2d.Extent();
 858       Handle(Geom2d_BSplineCurve) CurBS;
 859       CurBS = Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
 860       LOfBSpline2d.RemoveFirst();
 861       for (Standard_Integer ii = 2; ii <= NbC; ii++)
 862       {
 863         Handle(Geom2d_BSplineCurve) BS =
 864           Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
 865
 866         //Check for period jump in point of contact.
 867         gp_Pnt2d aC1End = CurBS->Pole(CurBS->NbPoles()); // End of C1.
 868         gp_Pnt2d aC2Beg = BS->Pole(1); // Beginning of C2.
 869         Standard_Real anUJump = 0.0, anVJump = 0.0;
 870
 871         if (anUPeriod > 0.0 &&
 872             Abs (aC1End.X() - aC2Beg.X()) > (anUPeriod ) / 2.01)
 873         {
 874           Standard_Real aMultCoeff =  aC2Beg.X() < aC1End.X() ? 1.0 : -1.0;
 875           anUJump = (anUPeriod) * aMultCoeff;
 876         }
 877
 878         if (anVPeriod &&
 879             Abs (aC1End.Y() - aC2Beg.Y()) > (anVPeriod) / 2.01)
 880         {
 881           Standard_Real aMultCoeff =  aC2Beg.Y() < aC1End.Y() ? 1.0 : -1.0;
 882           anVJump = (anVPeriod) * aMultCoeff;
 883         }
 884
 885         CurBS = Concat(CurBS,BS, anUJump, anVJump);
 886         LOfBSpline2d.RemoveFirst();
 887       }
 888       return CurBS;
 889     }
 890   }
 891
 892   if(InitialCurve2d.IsNull()) {
 893     AHC2d = BuildInitialCurve2d(Curve, S);
 894     if(!myProjIsDone)
 895       return Handle(Geom2d_BSplineCurve)();
 896   }
 897   return ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 898 }
 899
 900
 901 //=======================================================================
 902 //function : ProjLib_BuildInitialCurve2d
 903 //purpose  :
 904 //=======================================================================
 905
 906 Handle(Adaptor2d_HCurve2d)
 907      ProjLib_ComputeApproxOnPolarSurface::
 908      BuildInitialCurve2d(const Handle(Adaptor3d_HCurve)&   Curve,
 909                          const Handle(Adaptor3d_HSurface)& Surf)
 910 {
 911   //  discretize the Curve with quasiuniform deflection
 912   //  density at least NbOfPnts points
 913   myProjIsDone = Standard_False;
 914
 915   //OCC217
 916   Standard_Real Tol3d = myTolerance;
 917   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
 918   Standard_Real DistTol3d = 100.0*Tol3d;
 919   if(myMaxDist > 0.)
 920   {
 921     DistTol3d = myMaxDist;
 922   }
 923   Standard_Real DistTol3d2 = DistTol3d * DistTol3d;
 924   Standard_Real uperiod = 0.0, vperiod = 0.0;
 925   computePeriodicity(Surf, uperiod, vperiod);
 926
 927   // NO myTol is Tol2d !!!!
 928   //Standard_Real TolU = myTolerance, TolV = myTolerance;
 929   //Standard_Real Tol3d = 100*myTolerance; // At random Balthazar.
 930
 931   Standard_Integer NbOfPnts = 61;
 932   GCPnts_QuasiUniformAbscissa QUA(Curve->GetCurve(),NbOfPnts);
 933   NbOfPnts = QUA.NbPoints();
 934   TColgp_Array1OfPnt Pts(1,NbOfPnts);
 935   TColStd_Array1OfReal Param(1,NbOfPnts);
 936   Standard_Integer i, j;
 937   for( i = 1; i <= NbOfPnts ; i++ ) {
 938     Param(i) = QUA.Parameter(i);
 939     Pts(i) = Curve->Value(Param(i));
 940   }
 941
 942   TColgp_Array1OfPnt2d Pts2d(1,NbOfPnts);
 943   TColStd_Array1OfInteger Mult(1,NbOfPnts);
 944   Mult.Init(1);
 945   Mult(1) = Mult(NbOfPnts) = 2;
 946
 947   Standard_Real Uinf, Usup, Vinf, Vsup;
 948   Uinf = Surf->Surface().FirstUParameter();
 949   Usup = Surf->Surface().LastUParameter();
 950   Vinf = Surf->Surface().FirstVParameter();
 951   Vsup = Surf->Surface().LastVParameter();
 952   GeomAbs_SurfaceType Type = Surf->GetType();
 953   if((Type != GeomAbs_BSplineSurface) && (Type != GeomAbs_BezierSurface) &&
 954      (Type != GeomAbs_OffsetSurface)) {
 955     Standard_Real S, T;
 956 //    Standard_Integer usens = 0, vsens = 0;
 957     // to know the position relatively to the period
 958     switch (Type) {
 959 //    case GeomAbs_Plane:
 960 //      {
 961 //      gp_Pln Plane = Surf->Plane();
 962 //      for ( i = 1 ; i <= NbOfPnts ; i++) {
 963 //        ElSLib::Parameters( Plane, Pts(i), S, T);
 964 //        Pts2d(i).SetCoord(S,T);
 965 //      }
 966 //      myProjIsDone = Standard_True;
 967 //      break;
 968 //      }
 969     case GeomAbs_Cylinder:
 970       {
 971 //      Standard_Real Sloc, Tloc;
 972         Standard_Real Sloc;
 973         Standard_Integer usens = 0;
 974         gp_Cylinder Cylinder = Surf->Cylinder();
 975         ElSLib::Parameters( Cylinder, Pts(1), S, T);
 976         Pts2d(1).SetCoord(S,T);
 977         for ( i = 2 ; i <= NbOfPnts ; i++) {
 978           Sloc = S;
 979           ElSLib::Parameters( Cylinder, Pts(i), S, T);
 980           if(Abs(Sloc - S) > M_PI) {
 981             if(Sloc > S)
 982               usens++;
 983             else
 984               usens--;
 985           }
 986           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
 987         }
 988         myProjIsDone = Standard_True;
 989         break;
 990       }
 991     case GeomAbs_Cone:
 992       {
 993 //      Standard_Real Sloc, Tloc;
 994         Standard_Real Sloc;
 995         Standard_Integer usens = 0;
 996         gp_Cone Cone = Surf->Cone();
 997         ElSLib::Parameters( Cone, Pts(1), S, T);
 998         Pts2d(1).SetCoord(S,T);
 999         for ( i = 2 ; i <= NbOfPnts ; i++) {
1000           Sloc = S;
1001           ElSLib::Parameters( Cone, Pts(i), S, T);
1002           if(Abs(Sloc - S) > M_PI) {
1003             if(Sloc > S)
1004               usens++;
1005             else
1006               usens--;
1007           }
1008           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
1009         }
1010         myProjIsDone = Standard_True;
1011         break;
1012       }
1013     case GeomAbs_Sphere:
1014       {
1015         Standard_Real Sloc, Tloc;
1016         Standard_Integer usens = 0, vsens = 0; //usens steps by half-period
1017         Standard_Boolean vparit = Standard_False;
1018         gp_Sphere Sphere = Surf->Sphere();
1019         ElSLib::Parameters( Sphere, Pts(1), S, T);
1020         Pts2d(1).SetCoord(S,T);
1021         for ( i = 2 ; i <= NbOfPnts ; i++) {
1022           Sloc = S;Tloc = T;
1023           ElSLib::Parameters( Sphere, Pts(i), S, T);
1024           if(1.6*M_PI < Abs(Sloc - S)) {
1025             if(Sloc > S)
1026               usens += 2;
1027             else
1028               usens -= 2;
1029     }
1030           if(1.6*M_PI > Abs(Sloc - S) && Abs(Sloc - S) > 0.4*M_PI) {
1031             vparit = !vparit;
1032             if(Sloc > S)
1033               usens++;
1034             else
1035               usens--;
1036             if(Abs(Tloc - Vsup) < (Vsup - Vinf)/5)
1037               vsens++;
1038             else
1039               vsens--;
1040           }
1041           if(vparit) {
1042             Pts2d(i).SetCoord(S+usens*M_PI,(M_PI - T)*(vsens-1));
1043           }
1044           else {
1045             Pts2d(i).SetCoord(S+usens*M_PI,T+vsens*M_PI);
1046
1047           }
1048         }
1049         myProjIsDone = Standard_True;
1050         break;
1051       }
1052     case GeomAbs_Torus:
1053       {
1054         Standard_Real Sloc, Tloc;
1055         Standard_Integer usens = 0, vsens = 0;
1056         gp_Torus Torus = Surf->Torus();
1057         ElSLib::Parameters( Torus, Pts(1), S, T);
1058         Pts2d(1).SetCoord(S,T);
1059         for ( i = 2 ; i <= NbOfPnts ; i++) {
1060           Sloc = S; Tloc = T;
1061           ElSLib::Parameters( Torus, Pts(i), S, T);
1062           if(Abs(Sloc - S) > M_PI) {
1063             if(Sloc > S)
1064               usens++;
1065             else
1066               usens--;
1067     }
1068           if(Abs(Tloc - T) > M_PI) {
1069             if(Tloc > T)
1070               vsens++;
1071             else
1072               vsens--;
1073     }
1074           Pts2d(i).SetCoord(S+usens*2*M_PI,T+vsens*2*M_PI);
1075         }
1076         myProjIsDone = Standard_True;
1077         break;
1078       }
1079     default:
1080       throw Standard_NoSuchObject("ProjLib_ComputeApproxOnPolarSurface::BuildInitialCurve2d");
1081     }
1082   }
1083   else {
1084     myProjIsDone = Standard_False;
1085     Standard_Real Dist2Min = 1.e+200, u = 0., v = 0.;
1086     myDist = 0.;
1087     gp_Pnt pntproj;
1088
1089     TColgp_SequenceOfPnt2d Sols;
1090     Standard_Boolean areManyZeros = Standard_False;
1091
1092     Curve->D0(Param.Value(1), pntproj) ;
1093     Extrema_ExtPS  aExtPS(pntproj, Surf->Surface(), TolU, TolV) ;
1094     Standard_Real aMinSqDist = RealLast();
1095     if (aExtPS.IsDone())
1096     {
1097       for (i = 1; i <= aExtPS.NbExt(); i++)
1098       {
1099         Standard_Real aSqDist = aExtPS.SquareDistance(i);
1100         if (aSqDist < aMinSqDist)
1101           aMinSqDist = aSqDist;
1102       }
1103     }
1104     if (aMinSqDist > DistTol3d2) //try to project with less tolerance
1105     {
1106       TolU = Min(TolU, Precision::PConfusion());
1107       TolV = Min(TolV, Precision::PConfusion());
1108       aExtPS.Initialize(Surf->Surface(),
1109                         Surf->Surface().FirstUParameter(), Surf->Surface().LastUParameter(),
1110                         Surf->Surface().FirstVParameter(), Surf->Surface().LastVParameter(),
1111                         TolU, TolV);
1112       aExtPS.Perform(pntproj);
1113     }
1114
1115     if( aExtPS.IsDone() && aExtPS.NbExt() >= 1 ) {
1116
1117       Standard_Integer GoodValue = 1;
1118
1119       for ( i = 1 ; i <= aExtPS.NbExt() ; i++ ) {
1120         if( aExtPS.SquareDistance(i) < DistTol3d2 ) {
1121           if( aExtPS.SquareDistance(i) <= 1.e-18 ) {
1122             aExtPS.Point(i).Parameter(u,v);
1123             gp_Pnt2d p2d(u,v);
1124             Standard_Boolean isSame = Standard_False;
1125             for( j = 1; j <= Sols.Length(); j++ ) {
1126               if( p2d.SquareDistance( Sols.Value(j) ) <= 1.e-18 ) {
1127                 isSame = Standard_True;
1128                 break;
1129               }
1130             }
1131             if( !isSame ) Sols.Append( p2d );
1132           }
1133           if( Dist2Min > aExtPS.SquareDistance(i) ) {
1134             Dist2Min = aExtPS.SquareDistance(i);
1135             GoodValue = i;
1136           }
1137         }
1138       }
1139
1140       if( Sols.Length() > 1 ) areManyZeros = Standard_True;
1141
1142       if( Dist2Min <= DistTol3d2) {
1143         if( !areManyZeros ) {
1144           aExtPS.Point(GoodValue).Parameter(u,v);
1145           Pts2d(1).SetCoord(u,v);
1146           myProjIsDone = Standard_True;
1147         }
1148         else {
1149           Standard_Integer nbSols = Sols.Length();
1150           Standard_Real Dist2Max = -1.e+200;
1151           for( i = 1; i <= nbSols; i++ ) {
1152             const gp_Pnt2d& aP1 = Sols.Value(i);
1153             for( j = i+1; j <= nbSols; j++ ) {
1154               const gp_Pnt2d& aP2 = Sols.Value(j);
1155               Standard_Real aDist2 = aP1.SquareDistance(aP2);
1156               if( aDist2 > Dist2Max ) Dist2Max = aDist2;
1157             }
1158           }
1159       Standard_Real aMaxT2 = Max(TolU,TolV);
1160       aMaxT2 *= aMaxT2;
1161           if( Dist2Max > aMaxT2 ) {
1162             Standard_Integer tPp = 0;
1163             for( i = 1; i <= 5; i++ ) {
1164               Standard_Integer nbExtOk = 0;
1165               Standard_Integer indExt = 0;
1166               Standard_Integer iT = 1 + (NbOfPnts - 1)/5*i;
1167               Curve->D0( Param.Value(iT), pntproj );
1168               Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1169               Dist2Min = 1.e+200;
1170               if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1171                 for( j = 1 ; j <= aTPS.NbExt() ; j++ ) {
1172                   if( aTPS.SquareDistance(j) < DistTol3d2 ) {
1173                     nbExtOk++;
1174                     if( aTPS.SquareDistance(j) < Dist2Min ) {
1175                       Dist2Min = aTPS.SquareDistance(j);
1176                       indExt = j;
1177                     }
1178                   }
1179                 }
1180               }
1181               if( nbExtOk == 1 ) {
1182                 tPp = iT;
1183                 aTPS.Point(indExt).Parameter(u,v);
1184                 break;
1185               }
1186             }
1187
1188             if (tPp != 0 && tPp != NbOfPnts) {
1189               gp_Pnt2d aPp = gp_Pnt2d(u,v);
1190               gp_Pnt2d aPn;
1191               Standard_Boolean isFound = Standard_False;
1192               for (j = tPp + 1; j <= NbOfPnts; ++j)
1193               {
1194                 Curve->D0( Param.Value(j), pntproj );
1195                 Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1196                 Dist2Min = RealLast();
1197                 if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1198                   Standard_Integer indExt = 0;
1199                   for (i = 1; i <= aTPS.NbExt(); i++) {
1200                     if (aTPS.SquareDistance(i) < DistTol3d2 && aTPS.SquareDistance(i) < Dist2Min) {
1201                       Dist2Min = aTPS.SquareDistance(i);
1202                       indExt = i;
1203                     }
1204                   }
1205                   if (indExt > 0) {
1206                     aTPS.Point(indExt).Parameter(u,v);
1207                     aPn = gp_Pnt2d(u,v);
1208                     isFound = Standard_True;
1209                     break;
1210                   }
1211                 }
1212               }
1213
1214               if( isFound ) {
1215                 gp_Vec2d atV(aPp,aPn);
1216                 Standard_Boolean isChosen = Standard_False;
1217                 for( i = 1; i <= nbSols; i++ ) {
1218                   const gp_Pnt2d& aP1 = Sols.Value(i);
1219                   gp_Vec2d asV(aP1,aPp);
1220                   if( asV.Dot(atV) > 0. ) {
1221                     isChosen = Standard_True;
1222                     Pts2d(1).SetCoord(aP1.X(),aP1.Y());
1223                     myProjIsDone = Standard_True;
1224                     break;
1225                   }
1226                 }
1227                 if( !isChosen ) {
1228                   aExtPS.Point(GoodValue).Parameter(u,v);
1229                   Pts2d(1).SetCoord(u,v);
1230                   myProjIsDone = Standard_True;
1231                 }
1232               }
1233               else {
1234                 aExtPS.Point(GoodValue).Parameter(u,v);
1235                 Pts2d(1).SetCoord(u,v);
1236                 myProjIsDone = Standard_True;
1237               }
1238             }
1239             else {
1240               aExtPS.Point(GoodValue).Parameter(u,v);
1241               Pts2d(1).SetCoord(u,v);
1242               myProjIsDone = Standard_True;
1243             }
1244           }
1245           else {
1246             aExtPS.Point(GoodValue).Parameter(u,v);
1247             Pts2d(1).SetCoord(u,v);
1248             myProjIsDone = Standard_True;
1249           }
1250         }
1251       }
1252
1253       //  calculate the following points with GenLocate_ExtPS
1254       // (and store the result and each parameter in a sequence)
1255       Standard_Integer usens = 0, vsens = 0;
1256       // to know the position relatively to the period
1257       Standard_Real U0 = u, V0 = v, U1 = u, V1 = v;
1258       // U0 and V0 are the points in the initialized period
1259       // (period with u and v),
1260       // U1 and V1 are the points for construction of poles
1261       myDist = Dist2Min;
1262       for ( i = 2 ; i <= NbOfPnts ; i++)
1263         if(myProjIsDone) {
1264           myProjIsDone = Standard_False;
1265           Dist2Min = RealLast();
1266           Curve->D0(Param.Value(i), pntproj);
1267           Extrema_GenLocateExtPS  aLocateExtPS(Surf->Surface(), TolU, TolV);
1268           aLocateExtPS.Perform(pntproj, U0, V0);
1269
1270           if (aLocateExtPS.IsDone())
1271           {
1272             if (aLocateExtPS.SquareDistance() < DistTol3d2)
1273             {  //OCC217
1274               //if (aLocateExtPS.SquareDistance() < Tol3d * Tol3d) {
1275               if (aLocateExtPS.SquareDistance() > myDist)
1276               {
1277                 myDist = aLocateExtPS.SquareDistance();
1278               }
1279               (aLocateExtPS.Point()).Parameter(U0,V0);
1280               U1 = U0 + usens*uperiod;
1281               V1 = V0 + vsens*vperiod;
1282               Pts2d(i).SetCoord(U1,V1);
1283               myProjIsDone = Standard_True;
1284             }
1285             else
1286             {
1287               Extrema_ExtPS aGlobalExtr(pntproj, Surf->Surface(), TolU, TolV);
1288               if (aGlobalExtr.IsDone())
1289               {
1290                 Standard_Real LocalMinSqDist = RealLast();
1291                 Standard_Integer imin = 0;
1292                 for (Standard_Integer isol = 1; isol <= aGlobalExtr.NbExt(); isol++)
1293                 {
1294                   Standard_Real aSqDist = aGlobalExtr.SquareDistance(isol);
1295                   if (aSqDist < LocalMinSqDist)
1296                   {
1297                     LocalMinSqDist = aSqDist;
1298                     imin = isol;
1299                   }
1300                 }
1301                 if (LocalMinSqDist < DistTol3d2)
1302                 {
1303                   if (LocalMinSqDist > myDist)
1304                   {
1305                     myDist = LocalMinSqDist;
1306                   }
1307                   Standard_Real LocalU, LocalV;
1308                   aGlobalExtr.Point(imin).Parameter(LocalU, LocalV);
1309                   if (uperiod > 0. && Abs(U0 - LocalU) >= uperiod/2.)
1310                   {
1311                     if (LocalU > U0)
1312                       usens = -1;
1313                     else
1314                       usens = 1;
1315                   }
1316                   if (vperiod > 0. && Abs(V0 - LocalV) >= vperiod/2.)
1317                   {
1318                     if (LocalV > V0)
1319                       vsens = -1;
1320                     else
1321                       vsens = 1;
1322                   }
1323                   U0 = LocalU; V0 = LocalV;
1324                   U1 = U0 + usens*uperiod;
1325                   V1 = V0 + vsens*vperiod;
1326                   Pts2d(i).SetCoord(U1,V1);
1327                   myProjIsDone = Standard_True;
1328
1329                   if((i == 2) && (!IsEqual(uperiod, 0.0) || !IsEqual(vperiod, 0.0)))
1330                   {//Make 1st point more precise for periodic surfaces
1331                     const Standard_Integer aSize = 3;
1332                     const gp_Pnt2d aP(Pts2d(2));
1333                     Standard_Real aUpar[aSize], aVpar[aSize];
1334                     Pts2d(1).Coord(aUpar[1], aVpar[1]);
1335                     aUpar[0] = aUpar[1] - uperiod;
1336                     aUpar[2] = aUpar[1] + uperiod;
1337                     aVpar[0] = aVpar[1] - vperiod;
1338                     aVpar[2] = aVpar[1] + vperiod;
1339
1340                     Standard_Real aSQdistMin = RealLast();
1341                     Standard_Integer aBestUInd = 1, aBestVInd = 1;
1342                     const Standard_Integer  aSizeU = IsEqual(uperiod, 0.0) ? 1 : aSize,
1343                                             aSizeV = IsEqual(vperiod, 0.0) ? 1 : aSize;
1344                     for(Standard_Integer uInd = 0; uInd < aSizeU; uInd++)
1345                     {
1346                       for(Standard_Integer vInd = 0; vInd < aSizeV; vInd++)
1347                       {
1348                         Standard_Real aSQdist = aP.SquareDistance(gp_Pnt2d(aUpar[uInd], aVpar[vInd]));
1349                         if(aSQdist < aSQdistMin)
1350                         {
1351                           aSQdistMin = aSQdist;
1352                           aBestUInd = uInd;
1353                           aBestVInd = vInd;
1354                         }
1355                       }
1356                     }
1357
1358                     Pts2d(1).SetCoord(aUpar[aBestUInd], aVpar[aBestVInd]);
1359                   }//if(i == 2) condition
1360                 }
1361               }
1362             }
1363           }
1364           if(!myProjIsDone && uperiod) {
1365             Standard_Real aUinf, aUsup, Uaux;
1366             aUinf = Surf->Surface().FirstUParameter();
1367             aUsup = Surf->Surface().LastUParameter();
1368             if((aUsup - U0) > (U0 - aUinf))
1369               Uaux = 2*aUinf - U0 + uperiod;
1370             else
1371               Uaux = 2*aUsup - U0 - uperiod;
1372
1373             Extrema_GenLocateExtPS  locext(Surf->Surface(), TolU, TolV);
1374             locext.Perform(pntproj, Uaux, V0);
1375
1376             if (locext.IsDone())
1377               if (locext.SquareDistance() < DistTol3d2) {  //OCC217
1378               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1379                 if (locext.SquareDistance() > myDist)
1380                 {
1381                   myDist = locext.SquareDistance();
1382                 }
1383                 (locext.Point()).Parameter(u,v);
1384                 if((aUsup - U0) > (U0 - aUinf))
1385                   usens--;
1386                 else
1387                   usens++;
1388                 U0 = u; V0 = v;
1389                 U1 = U0 + usens*uperiod;
1390                 V1 = V0 + vsens*vperiod;
1391                 Pts2d(i).SetCoord(U1,V1);
1392                 myProjIsDone = Standard_True;
1393               }
1394           }
1395           if(!myProjIsDone && vperiod) {
1396             Standard_Real aVinf, aVsup, Vaux;
1397             aVinf = Surf->Surface().FirstVParameter();
1398             aVsup = Surf->Surface().LastVParameter();
1399             if((aVsup - V0) > (V0 - aVinf))
1400               Vaux = 2*aVinf - V0 + vperiod;
1401             else
1402               Vaux = 2*aVsup - V0 - vperiod;
1403
1404             Extrema_GenLocateExtPS  locext(Surf->Surface(), TolU, TolV);
1405             locext.Perform(pntproj, U0, Vaux);
1406
1407             if (locext.IsDone())
1408               if (locext.SquareDistance() < DistTol3d2) {  //OCC217
1409               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1410                 if (locext.SquareDistance() > myDist)
1411                 {
1412                   myDist = locext.SquareDistance();
1413                 }
1414                 (locext.Point()).Parameter(u, v);
1415                 if((aVsup - V0) > (V0 - aVinf))
1416                   vsens--;
1417                 else
1418                   vsens++;
1419                 U0 = u; V0 = v;
1420                 U1 = U0 + usens*uperiod;
1421                 V1 = V0 + vsens*vperiod;
1422                 Pts2d(i).SetCoord(U1,V1);
1423                 myProjIsDone = Standard_True;
1424               }
1425           }
1426           if(!myProjIsDone && uperiod && vperiod) {
1427             Standard_Real Uaux, Vaux;
1428             if((Usup - U0) > (U0 - Uinf))
1429               Uaux = 2*Uinf - U0 + uperiod;
1430             else
1431               Uaux = 2*Usup - U0 - uperiod;
1432             if((Vsup - V0) > (V0 - Vinf))
1433               Vaux = 2*Vinf - V0 + vperiod;
1434             else
1435               Vaux = 2*Vsup - V0 - vperiod;
1436
1437             Extrema_GenLocateExtPS  locext(Surf->Surface(), TolU, TolV);
1438             locext.Perform(pntproj, Uaux, Vaux);
1439
1440             if (locext.IsDone())
1441               if (locext.SquareDistance() < DistTol3d2) {
1442               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1443                 if (locext.SquareDistance() > myDist)
1444                 {
1445                   myDist = locext.SquareDistance();
1446                 }
1447                 (locext.Point()).Parameter(u, v);
1448                 if((Usup - U0) > (U0 - Uinf))
1449                   usens--;
1450                 else
1451                   usens++;
1452                 if((Vsup - V0) > (V0 - Vinf))
1453                   vsens--;
1454                 else
1455                   vsens++;
1456                 U0 = u; V0 = v;
1457                 U1 = U0 + usens*uperiod;
1458                 V1 = V0 + vsens*vperiod;
1459                 Pts2d(i).SetCoord(U1,V1);
1460                 myProjIsDone = Standard_True;
1461               }
1462           }
1463           if(!myProjIsDone) {
1464             Extrema_ExtPS ext(pntproj, Surf->Surface(), TolU, TolV) ;
1465             if (ext.IsDone()) {
1466               Dist2Min = ext.SquareDistance(1);
1467               Standard_Integer aGoodValue = 1;
1468               for ( j = 2 ; j <= ext.NbExt() ; j++ )
1469                 if( Dist2Min > ext.SquareDistance(j)) {
1470                   Dist2Min = ext.SquareDistance(j);
1471                   aGoodValue = j;
1472                 }
1473               if (Dist2Min < DistTol3d2) {
1474               //if (Dist2Min < Tol3d * Tol3d) {
1475                 if (Dist2Min > myDist)
1476                 {
1477                   myDist = Dist2Min;
1478                 }
1479                 (ext.Point(aGoodValue)).Parameter(u, v);
1480                 if(uperiod) {
1481                   if((U0 - u) > (2*uperiod/3)) {
1482                     usens++;
1483                   }
1484                   else
1485                     if((u - U0) > (2*uperiod/3)) {
1486                       usens--;
1487                     }
1488     }
1489                 if(vperiod) {
1490                   if((V0 - v) > (vperiod/2)) {
1491                     vsens++;
1492                   }
1493                   else
1494                     if((v - V0) > (vperiod/2)) {
1495                       vsens--;
1496                     }
1497       }
1498                 U0 = u; V0 = v;
1499                 U1 = U0 + usens*uperiod;
1500                 V1 = V0 + vsens*vperiod;
1501                 Pts2d(i).SetCoord(U1,V1);
1502                 myProjIsDone = Standard_True;
1503               }
1504             }
1505           }
1506         }
1507         else break;
1508     }
1509   }
1510   // -- Pnts2d is transformed into Geom2d_BSplineCurve, with the help of Param and Mult
1511   if(myProjIsDone) {
1512     myBSpline = new Geom2d_BSplineCurve(Pts2d,Param,Mult,1);
1513     //jgv: put the curve into parametric range
1514     gp_Pnt2d MidPoint = myBSpline->Value(0.5*(myBSpline->FirstParameter() + myBSpline->LastParameter()));
1515     Standard_Real TestU = MidPoint.X(), TestV = MidPoint.Y();
1516     Standard_Real sense = 0.;
1517     if (uperiod)
1518     {
1519       if (TestU < Uinf - TolU)
1520         sense = 1.;
1521       else if (TestU > Usup + TolU)
1522         sense = -1;
1523       while (TestU < Uinf - TolU || TestU > Usup + TolU)
1524         TestU += sense * uperiod;
1525     }
1526     if (vperiod)
1527     {
1528       sense = 0.;
1529       if (TestV < Vinf - TolV)
1530         sense = 1.;
1531       else if (TestV > Vsup + TolV)
1532         sense = -1.;
1533       while (TestV < Vinf - TolV || TestV > Vsup + TolV)
1534         TestV += sense * vperiod;
1535     }
1536     gp_Vec2d Offset(TestU - MidPoint.X(), TestV - MidPoint.Y());
1537     if (Abs(Offset.X()) > gp::Resolution() ||
1538         Abs(Offset.Y()) > gp::Resolution())
1539       myBSpline->Translate(Offset);
1540     //////////////////////////////////////////
1541     Geom2dAdaptor_Curve GAC(myBSpline);
1542     Handle(Adaptor2d_HCurve2d) IC2d = new Geom2dAdaptor_HCurve(GAC);
1543 #ifdef OCCT_DEBUG
1544 //    char name [100];
1545 //    sprintf(name,"%s_%d","build",compteur++);
1546 //    DrawTrSurf::Set(name,myBSpline);
1547 #endif
1548     return IC2d;
1549   }
1550   else {
1551 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:50 2002 Begin
1552 //     Standard_NoSuchObject_Raise_if(1,"ProjLib_Compu: build echec");
1553 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:51 2002 End
1554     return Handle(Adaptor2d_HCurve2d)();
1555   }
1556 //  myProjIsDone = Standard_False;
1557 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:01 2002 Begin
1558 //   Standard_NoSuchObject_Raise_if(1,"ProjLib_ComputeOnPS: build echec");
1559 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:02 2002 End
1560 }
1561
1562 //=======================================================================
1563 //function : ProjLib_ProjectUsingInitialCurve2d
1564 //purpose  :
1565 //=======================================================================
1566 Handle(Geom2d_BSplineCurve)
1567      ProjLib_ComputeApproxOnPolarSurface::
1568      ProjectUsingInitialCurve2d(const Handle(Adaptor3d_HCurve)& Curve,
1569                                 const Handle(Adaptor3d_HSurface)& Surf,
1570                                 const Handle(Adaptor2d_HCurve2d)& InitCurve2d)
1571 {
1572   //OCC217
1573   Standard_Real Tol3d = myTolerance;
1574   Standard_Real DistTol3d = 100.0*Tol3d;
1575   if(myMaxDist > 0.)
1576   {
1577     DistTol3d = myMaxDist;
1578   }
1579   Standard_Real DistTol3d2 = DistTol3d * DistTol3d;
1580   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
1581   Standard_Real Tol2d = Max(Sqrt(TolU*TolU + TolV*TolV), Precision::PConfusion());
1582
1583   Standard_Integer i;
1584   GeomAbs_SurfaceType TheTypeS = Surf->GetType();
1585   GeomAbs_CurveType TheTypeC = Curve->GetType();
1586   if(TheTypeS == GeomAbs_Plane) {
1587     Standard_Real S, T;
1588     gp_Pln Plane = Surf->Plane();
1589     if(TheTypeC == GeomAbs_BSplineCurve) {
1590       myTolReached = Precision::Confusion();
1591       Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1592       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1593       for(i = 1;i <= Curve->NbPoles();i++) {
1594         ElSLib::Parameters( Plane, BSC->Pole(i), S, T);
1595         Poles2d(i).SetCoord(S,T);
1596       }
1597       TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1598       BSC->Knots(Knots);
1599       TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1600       BSC->Multiplicities(Mults);
1601       if(BSC->IsRational()) {
1602         TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1603         BSC->Weights(Weights);
1604         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1605                                        BSC->Degree(), BSC->IsPeriodic()) ;
1606       }
1607       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1608                                      BSC->Degree(), BSC->IsPeriodic()) ;
1609
1610     }
1611     if(TheTypeC == GeomAbs_BezierCurve) {
1612       myTolReached = Precision::Confusion();
1613       Handle(Geom_BezierCurve) BC = Curve->Bezier();
1614       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1615       for(i = 1;i <= Curve->NbPoles();i++) {
1616         ElSLib::Parameters( Plane, BC->Pole(i), S, T);
1617         Poles2d(i).SetCoord(S,T);
1618       }
1619       TColStd_Array1OfReal Knots(1, 2);
1620       Knots.SetValue(1,0.0);
1621       Knots.SetValue(2,1.0);
1622       TColStd_Array1OfInteger Mults(1, 2);
1623       Mults.Init(BC->NbPoles());
1624       if(BC->IsRational()) {
1625         TColStd_Array1OfReal Weights(1, BC->NbPoles());
1626         BC->Weights(Weights);
1627         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1628                                        BC->Degree(), BC->IsPeriodic()) ;
1629       }
1630       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1631                                      BC->Degree(), BC->IsPeriodic()) ;
1632     }
1633   }
1634   if(TheTypeS == GeomAbs_BSplineSurface) {
1635     Handle(Geom_BSplineSurface) BSS = Surf->BSpline();
1636     if((BSS->MaxDegree() == 1) &&
1637        (BSS->NbUPoles() == 2) &&
1638        (BSS->NbVPoles() == 2)) {
1639       gp_Pnt p11 = BSS->Pole(1,1);
1640       gp_Pnt p12 = BSS->Pole(1,2);
1641       gp_Pnt p21 = BSS->Pole(2,1);
1642       gp_Pnt p22 = BSS->Pole(2,2);
1643       gp_Vec V1(p11,p12);
1644       gp_Vec V2(p21,p22);
1645       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){
1646         Standard_Integer Dist2Min = IntegerLast();
1647         Standard_Real u,v;
1648         if(TheTypeC == GeomAbs_BSplineCurve) {
1649           myTolReached = Tol3d;
1650           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1651           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1652           for(i = 1;i <= Curve->NbPoles();i++) {
1653             myProjIsDone = Standard_False;
1654             Dist2Min = IntegerLast();
1655
1656             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1657             extrloc.Perform(BSC->Pole(i), (p11.X()+p22.X())/2, (p11.Y()+p22.Y())/2);
1658
1659             if (extrloc.IsDone()) {
1660               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1661               if (Dist2Min < DistTol3d2) {
1662                 (extrloc.Point()).Parameter(u,v);
1663                 Poles2d(i).SetCoord(u,v);
1664                 myProjIsDone = Standard_True;
1665               }
1666               else break;
1667             }
1668             else break;
1669             if(!myProjIsDone)
1670               break;
1671           }
1672           if(myProjIsDone) {
1673             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1674             BSC->Knots(Knots);
1675             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1676             BSC->Multiplicities(Mults);
1677             if(BSC->IsRational()) {
1678               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1679               BSC->Weights(Weights);
1680               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1681                                              BSC->Degree(), BSC->IsPeriodic()) ;
1682             }
1683             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1684                                            BSC->Degree(), BSC->IsPeriodic()) ;
1685
1686
1687           }
1688         }
1689         if(TheTypeC == GeomAbs_BezierCurve) {
1690           myTolReached = Tol3d;
1691           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1692           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1693           for(i = 1;i <= Curve->NbPoles();i++) {
1694             Dist2Min = IntegerLast();
1695
1696             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1697             extrloc.Perform(BC->Pole(i), 0.5, 0.5);
1698
1699             if (extrloc.IsDone()) {
1700               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1701               if (Dist2Min < DistTol3d2) {
1702                 (extrloc.Point()).Parameter(u,v);
1703                 Poles2d(i).SetCoord(u,v);
1704                 myProjIsDone = Standard_True;
1705               }
1706               else break;
1707             }
1708             else break;
1709             if(myProjIsDone)
1710               myProjIsDone = Standard_False;
1711             else break;
1712           }
1713           if(myProjIsDone) {
1714             TColStd_Array1OfReal Knots(1, 2);
1715             Knots.SetValue(1,0.0);
1716             Knots.SetValue(2,1.0);
1717             TColStd_Array1OfInteger Mults(1, 2);
1718             Mults.Init(BC->NbPoles());
1719             if(BC->IsRational()) {
1720               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1721               BC->Weights(Weights);
1722               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1723                                                     BC->Degree(), BC->IsPeriodic()) ;
1724             }
1725             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1726                                                   BC->Degree(), BC->IsPeriodic()) ;
1727           }
1728         }
1729       }
1730     }
1731   }
1732   else if(TheTypeS == GeomAbs_BezierSurface) {
1733     Handle(Geom_BezierSurface) BS = Surf->Bezier();
1734     if((BS->MaxDegree() == 1) &&
1735        (BS->NbUPoles() == 2) &&
1736        (BS->NbVPoles() == 2)) {
1737       gp_Pnt p11 = BS->Pole(1,1);
1738       gp_Pnt p12 = BS->Pole(1,2);
1739       gp_Pnt p21 = BS->Pole(2,1);
1740       gp_Pnt p22 = BS->Pole(2,2);
1741       gp_Vec V1(p11,p12);
1742       gp_Vec V2(p21,p22);
1743       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){
1744         Standard_Integer Dist2Min = IntegerLast();
1745         Standard_Real u,v;
1746
1747 //      gp_Pnt pntproj;
1748         if(TheTypeC == GeomAbs_BSplineCurve) {
1749           myTolReached = Tol3d;
1750           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1751           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1752           for(i = 1;i <= Curve->NbPoles();i++) {
1753             myProjIsDone = Standard_False;
1754             Dist2Min = IntegerLast();
1755
1756             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1757             extrloc.Perform(BSC->Pole(i), (p11.X()+p22.X())/2, (p11.Y()+p22.Y())/2);
1758
1759             if (extrloc.IsDone()) {
1760               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1761               if (Dist2Min < DistTol3d2) {
1762                 (extrloc.Point()).Parameter(u,v);
1763                 Poles2d(i).SetCoord(u,v);
1764                 myProjIsDone = Standard_True;
1765               }
1766               else break;
1767             }
1768             else break;
1769             if(!myProjIsDone)
1770               break;
1771           }
1772           if(myProjIsDone) {
1773             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1774             BSC->Knots(Knots);
1775             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1776             BSC->Multiplicities(Mults);
1777             if(BSC->IsRational()) {
1778               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1779               BSC->Weights(Weights);
1780               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1781                                                     BSC->Degree(), BSC->IsPeriodic()) ;
1782             }
1783             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1784                                                   BSC->Degree(), BSC->IsPeriodic()) ;
1785
1786
1787           }
1788         }
1789         if(TheTypeC == GeomAbs_BezierCurve) {
1790           myTolReached = Tol3d;
1791           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1792           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1793           for(i = 1;i <= Curve->NbPoles();i++) {
1794             Dist2Min = IntegerLast();
1795
1796             Extrema_GenLocateExtPS  extrloc(Surf->Surface(), TolU, TolV);
1797             extrloc.Perform(BC->Pole(i), 0.5, 0.5);
1798
1799             if (extrloc.IsDone()) {
1800               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1801               if (Dist2Min < DistTol3d2) {
1802                 (extrloc.Point()).Parameter(u,v);
1803                 Poles2d(i).SetCoord(u,v);
1804                 myProjIsDone = Standard_True;
1805               }
1806               else break;
1807             }
1808             else break;
1809             if(myProjIsDone)
1810               myProjIsDone = Standard_False;
1811             else break;
1812           }
1813           if(myProjIsDone) {
1814             TColStd_Array1OfReal Knots(1, 2);
1815             Knots.SetValue(1,0.0);
1816             Knots.SetValue(2,1.0);
1817             TColStd_Array1OfInteger Mults(1, 2);
1818             Mults.Init(BC->NbPoles());
1819             if(BC->IsRational()) {
1820               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1821               BC->Weights(Weights);
1822               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1823                                                     BC->Degree(), BC->IsPeriodic()) ;
1824             }
1825             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1826                                                   BC->Degree(), BC->IsPeriodic()) ;
1827           }
1828         }
1829       }
1830     }
1831   }
1832
1833   ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, Tol3d) ;
1834
1835 #ifdef OCCT_DEBUG
1836   Standard_Integer Nb = 50;
1837
1838   Standard_Real U, U1, U2;
1839   U1 = F.FirstParameter();
1840   U2 = F.LastParameter();
1841
1842   TColgp_Array1OfPnt2d    DummyPoles(1,Nb+1);
1843   TColStd_Array1OfReal    DummyKnots(1,Nb+1);
1844   TColStd_Array1OfInteger DummyMults(1,Nb+1);
1845   DummyMults.Init(1);
1846   DummyMults(1) = 2;
1847   DummyMults(Nb+1) = 2;
1848   for (Standard_Integer ij = 0; ij <= Nb; ij++) {
1849     U = (Nb-ij)*U1 + ij*U2;
1850     U /= Nb;
1851     DummyPoles(ij+1) = F.Value(U);
1852     DummyKnots(ij+1) = ij;
1853   }
1854   Handle(Geom2d_BSplineCurve) DummyC2d =
1855     new Geom2d_BSplineCurve(DummyPoles, DummyKnots, DummyMults, 1);
1856 #ifdef DRAW
1857   Standard_CString Temp = "bs2d";
1858   DrawTrSurf::Set(Temp,DummyC2d);
1859 #endif
1860 //  DrawTrSurf::Set((Standard_CString ) "bs2d",DummyC2d);
1861   Handle(Geom2dAdaptor_HCurve) DDD =
1862     Handle(Geom2dAdaptor_HCurve)::DownCast(InitCurve2d);
1863
1864 #ifdef DRAW
1865   Temp = "initc2d";
1866   DrawTrSurf::Set(Temp,DDD->ChangeCurve2d().Curve());
1867 #endif
1868 //  DrawTrSurf::Set((Standard_CString ) "initc2d",DDD->ChangeCurve2d().Curve());
1869 #endif
1870
1871   Standard_Integer Deg1,Deg2;
1872   Deg1 = 2;
1873   Deg2 = 8;
1874   if(myDegMin > 0)
1875   {
1876     Deg1 = myDegMin;
1877   }
1878   if(myDegMax > 0)
1879   {
1880     Deg2 = myDegMax;
1881   }
1882   Standard_Integer aMaxSegments = 1000;
1883   if(myMaxSegments > 0)
1884   {
1885     aMaxSegments = myMaxSegments;
1886   }
1887   AppParCurves_Constraint aFistC = AppParCurves_TangencyPoint, aLastC = AppParCurves_TangencyPoint;
1888   if(myBndPnt != AppParCurves_TangencyPoint)
1889   {
1890     aFistC = myBndPnt;
1891     aLastC = myBndPnt;
1892   }
1893
1894   if (myDist > 10.*Tol3d)
1895   {
1896     aFistC = AppParCurves_PassPoint;
1897     aLastC = AppParCurves_PassPoint;
1898   }
1899
1900   Approx_FitAndDivide2d Fit(Deg1, Deg2, Tol3d, Tol2d, Standard_True, aFistC, aLastC);
1901   Fit.SetMaxSegments(aMaxSegments);
1902   if (InitCurve2d->GetType() == GeomAbs_Line)
1903   {
1904     Fit.SetInvOrder(Standard_False);
1905   }
1906   Fit.Perform(F);
1907
1908   Standard_Real anOldTol2d = Tol2d;
1909   Standard_Real aNewTol2d = 0;
1910   if(Fit.IsAllApproximated()) {
1911     Standard_Integer j;
1912     Standard_Integer NbCurves = Fit.NbMultiCurves();
1913     Standard_Integer MaxDeg = 0;
1914     // To transform the MultiCurve into BSpline, it is required that all
1915     // Bezier constituing it have the same degree -> Calculation of MaxDeg
1916     Standard_Integer NbPoles  = 1;
1917     for (j = 1; j <= NbCurves; j++) {
1918       Standard_Integer Deg = Fit.Value(j).Degree();
1919       MaxDeg = Max ( MaxDeg, Deg);
1920       Fit.Error(j,Tol3d, Tol2d);
1921       aNewTol2d = Max(aNewTol2d, Tol2d);
1922     }
1923     //
1924     myTolReached = Max(myTolReached, myTolerance * (aNewTol2d / anOldTol2d));
1925     //
1926     NbPoles = MaxDeg * NbCurves + 1;               //Tops on the BSpline
1927     TColgp_Array1OfPnt2d  Poles( 1, NbPoles);
1928
1929     TColgp_Array1OfPnt2d TempPoles( 1, MaxDeg + 1);//to augment the degree
1930
1931     TColStd_Array1OfReal Knots( 1, NbCurves + 1);  //Nodes of the BSpline
1932
1933     Standard_Integer Compt = 1;
1934     for (i = 1; i <= NbCurves; i++) {
1935       Fit.Parameters(i, Knots(i), Knots(i+1));
1936       AppParCurves_MultiCurve MC = Fit.Value( i);       //Load the Ith Curve
1937       TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Retrieve the tops
1938       MC.Curve(1, Poles2d);
1939
1940       //Eventual augmentation of the degree
1941       Standard_Integer Inc = MaxDeg - MC.Degree();
1942       if ( Inc > 0) {
1943 //      BSplCLib::IncreaseDegree( Inc, Poles2d, PLib::NoWeights(),
1944         BSplCLib::IncreaseDegree( MaxDeg, Poles2d, BSplCLib::NoWeights(),
1945                          TempPoles, BSplCLib::NoWeights());
1946         //update of tops of the PCurve
1947         for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1948           Poles.SetValue( Compt, TempPoles( k));
1949           Compt++;
1950         }
1951       }
1952       else {
1953         //update of tops of the PCurve
1954         for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1955           Poles.SetValue( Compt, Poles2d( k));
1956           Compt++;
1957         }
1958       }
1959
1960       Compt--;
1961     }
1962
1963     //update of fields of  ProjLib_Approx
1964     Standard_Integer NbKnots = NbCurves + 1;
1965
1966     TColStd_Array1OfInteger   Mults( 1, NbKnots);
1967     Mults.Init(MaxDeg);
1968     Mults.SetValue( 1, MaxDeg + 1);
1969     Mults.SetValue(NbKnots, MaxDeg + 1);
1970     myProjIsDone = Standard_True;
1971     Handle(Geom2d_BSplineCurve) Dummy =
1972       new Geom2d_BSplineCurve(Poles,Knots,Mults,MaxDeg);
1973
1974     // try to smoother the Curve GeomAbs_C1.
1975
1976     Standard_Boolean OK = Standard_True;
1977     Standard_Real aSmoothTol = Max(Precision::Confusion(), aNewTol2d);
1978     if (myBndPnt == AppParCurves_PassPoint)
1979     {
1980       aSmoothTol *= 10.;
1981     }
1982     for (Standard_Integer ij = 2; ij < NbKnots; ij++) {
1983       OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1, aSmoothTol);
1984     }
1985 #ifdef OCCT_DEBUG
1986     if (!OK) {
1987       std::cout << "ProjLib_ComputeApproxOnPolarSurface : Smoothing echoue"<<std::endl;
1988     }
1989 #endif
1990     return Dummy;
1991   }
1992   return Handle(Geom2d_BSplineCurve)();
1993 }
1994
1995 //=======================================================================
1996 //function : BSpline
1997 //purpose  :
1998 //=======================================================================
1999
2000 Handle(Geom2d_BSplineCurve)
2001      ProjLib_ComputeApproxOnPolarSurface::BSpline() const
2002
2003 {
2004   return myBSpline ;
2005 }
2006
2007 //=======================================================================
2008 //function : Curve2d
2009 //purpose  :
2010 //=======================================================================
2011
2012 Handle(Geom2d_Curve)
2013      ProjLib_ComputeApproxOnPolarSurface::Curve2d() const
2014
2015 {
2016   Standard_NoSuchObject_Raise_if
2017     (!myProjIsDone,
2018      "ProjLib_ComputeApproxOnPolarSurface:2ndCurve2d");
2019   return my2ndCurve ;
2020 }
2021
2022
2023 //=======================================================================
2024 //function : IsDone
2025 //purpose  :
2026 //=======================================================================
2027
2028 Standard_Boolean ProjLib_ComputeApproxOnPolarSurface::IsDone() const
2029
2030 {
2031   return myProjIsDone;
2032 }
2033 //=======================================================================
2034 //function : SetTolerance
2035 //purpose  :
2036 //=======================================================================
2037
2038 void ProjLib_ComputeApproxOnPolarSurface::SetTolerance(const Standard_Real theTol)
2039
2040 {
2041   myTolerance = theTol;
2042 }
2043 //=======================================================================
2044 //function : SetDegree
2045 //purpose  :
2046 //=======================================================================
2047 void ProjLib_ComputeApproxOnPolarSurface::SetDegree(
2048                                        const Standard_Integer theDegMin,
2049                                        const Standard_Integer theDegMax)
2050 {
2051   myDegMin = theDegMin;
2052   myDegMax = theDegMax;
2053 }
2054 //=======================================================================
2055 //function : SetMaxSegments
2056 //purpose  :
2057 //=======================================================================
2058 void ProjLib_ComputeApproxOnPolarSurface::SetMaxSegments(
2059                                    const Standard_Integer theMaxSegments)
2060 {
2061   myMaxSegments = theMaxSegments;
2062 }
2063
2064 //=======================================================================
2065 //function : SetBndPnt
2066 //purpose  :
2067 //=======================================================================
2068 void ProjLib_ComputeApproxOnPolarSurface::SetBndPnt(
2069                                  const AppParCurves_Constraint theBndPnt)
2070 {
2071   myBndPnt = theBndPnt;
2072 }
2073
2074 //=======================================================================
2075 //function : SetMaxDist
2076 //purpose  :
2077 //=======================================================================
2078 void ProjLib_ComputeApproxOnPolarSurface::SetMaxDist(
2079                                           const Standard_Real theMaxDist)
2080 {
2081   myMaxDist = theMaxDist;
2082 }
2083
2084 //=======================================================================
2085 //function : Tolerance
2086 //purpose  :
2087 //=======================================================================
2088
2089 Standard_Real ProjLib_ComputeApproxOnPolarSurface::Tolerance() const
2090
2091 {
2092   return myTolReached;
2093 }
2094