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       Standard_NoSuchObject::Raise("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     Standard_NoSuchObject::Raise("");
 367   }
 368
 369   // Try to run simple search with initial point (U0, V0).
 370   Extrema_GenLocateExtPS  locext(p, SurfLittle, U0, V0, theData.myTolU, theData.myTolV);
 371   if (locext.IsDone())
 372   {
 373     locext.Point().Parameter(u, v);
 374     Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
 375     if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
 376         locext.SquareDistance() < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
 377     {
 378       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
 379       return pnt;
 380     }
 381   }
 382
 383   // Perform whole param space search.
 384   Extrema_ExtPS  ext(p, SurfLittle, theData.myTolU, theData.myTolV);
 385   if (ext.IsDone() && ext.NbExt() >= 1)
 386   {
 387     Dist2Min = ext.SquareDistance(1);
 388     Standard_Integer GoodValue = 1;
 389     for (Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ )
 390     {
 391       if( Dist2Min > ext.SquareDistance(i))
 392       {
 393         Dist2Min = ext.SquareDistance(i);
 394         GoodValue = i;
 395       }
 396     }
 397     ext.Point(GoodValue).Parameter(u, v);
 398     Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
 399     if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
 400         ext.SquareDistance(GoodValue) < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
 401     {
 402       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
 403       return pnt;
 404     }
 405   }
 406
 407   // Both searches return bad values, use point from initial 2dcurve.
 408   return p2d;
 409 }
 410
 411
 412 //=======================================================================
 413 //function : ProjLib_PolarFunction
 414 //purpose  : (OCC217 - apo)- This class produce interface to call "gp_Pnt2d Function_Value(...)"
 415 //=======================================================================
 416
 417 class ProjLib_PolarFunction : public AppCont_Function
 418 {
 419   aFuncStruct myStruct;
 420
 421   public :
 422
 423   ProjLib_PolarFunction(const Handle(Adaptor3d_HCurve) & C,
 424                         const Handle(Adaptor3d_HSurface)& Surf,
 425                         const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 426                         const Standard_Real Tol3d)
 427   {
 428     myNbPnt = 0;
 429     myNbPnt2d = 1;
 430
 431     computePeriodicity(Surf, myStruct.myPeriod[0], myStruct.myPeriod[1]);
 432
 433     myStruct.myCurve = C;
 434     myStruct.myInitCurve2d = InitialCurve2d;
 435     myStruct.mySurf = Surf;
 436     myStruct.mySqProjOrtTol = 10000.0 * Tol3d * Tol3d;
 437     myStruct.myTolU = Surf->UResolution(Tol3d);
 438     myStruct.myTolV = Surf->VResolution(Tol3d);
 439   }
 440
 441   ~ProjLib_PolarFunction() {}
 442
 443   Standard_Real FirstParameter() const
 444   {
 445     return myStruct.myCurve->FirstParameter();
 446   }
 447
 448   Standard_Real LastParameter() const
 449   {
 450     return myStruct.myCurve->LastParameter();
 451   }
 452
 453   gp_Pnt2d Value(const Standard_Real t) const
 454   {
 455     return Function_Value(t, myStruct);
 456   }
 457
 458   Standard_Boolean Value(const Standard_Real   theT,
 459                          NCollection_Array1<gp_Pnt2d>& thePnt2d,
 460                          NCollection_Array1<gp_Pnt>&   /*thePnt*/) const
 461   {
 462     thePnt2d(1) = Function_Value(theT, myStruct);
 463     return Standard_True;
 464   }
 465
 466   Standard_Boolean D1(const Standard_Real   /*theT*/,
 467                       NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
 468                       NCollection_Array1<gp_Vec>&   /*theVec*/) const
 469     {return Standard_False;}
 470 };
 471
 472 //=======================================================================
 473 //function : ProjLib_ComputeApproxOnPolarSurface
 474 //purpose  :
 475 //=======================================================================
 476
 477 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface()
 478 : myProjIsDone(Standard_False),
 479   myTolerance (-1.0)
 480 {
 481 }
 482
 483
 484 //=======================================================================
 485 //function : ProjLib_ComputeApproxOnPolarSurface
 486 //purpose  :
 487 //=======================================================================
 488
 489 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 490                     (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
 491                      const Handle(Adaptor3d_HCurve)&   theCurve,
 492                      const Handle(Adaptor3d_HSurface)& theSurface,
 493                      const Standard_Real               theTolerance3D)
 494 : myProjIsDone(Standard_False),
 495   myTolerance (theTolerance3D)
 496 {
 497   myBSpline = Perform(theInitialCurve2d, theCurve, theSurface);
 498 }
 499 //=======================================================================
 500 //function : ProjLib_ComputeApproxOnPolarSurface
 501 //purpose  : Process the case of sewing
 502 //=======================================================================
 503
 504 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 505                 (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
 506                  const Handle(Adaptor2d_HCurve2d)& theInitialCurve2dBis,
 507                  const Handle(Adaptor3d_HCurve)&   theCurve,
 508                  const Handle(Adaptor3d_HSurface)& theSurface,
 509                  const Standard_Real               theTolerance3D)
 510 : myProjIsDone(Standard_False),
 511   myTolerance (theTolerance3D)
 512 {
 513   // InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing
 514   Handle(Geom2d_BSplineCurve) bsc =
 515     Perform(theInitialCurve2d, theCurve, theSurface);
 516
 517   if(myProjIsDone) {
 518     gp_Pnt2d P2dproj, P2d, P2dBis;
 519     P2dproj = bsc->StartPoint();
 520     P2d    = theInitialCurve2d->Value(theInitialCurve2d->FirstParameter());
 521     P2dBis = theInitialCurve2dBis->Value(theInitialCurve2dBis->FirstParameter());
 522
 523     Standard_Real Dist, DistBis;
 524     Dist    = P2dproj.Distance(P2d);
 525     DistBis = P2dproj.Distance(P2dBis);
 526     if( Dist < DistBis)  {
 527       // myBSpline2d is the pcurve that is found. It is translated to obtain myCurve2d
 528       myBSpline = bsc;
 529       Handle(Geom2d_Geometry) GG = myBSpline->Translated(P2d, P2dBis);
 530       my2ndCurve = Handle(Geom2d_Curve)::DownCast(GG);
 531     }
 532     else {
 533       my2ndCurve = bsc;
 534       Handle(Geom2d_Geometry) GG = my2ndCurve->Translated(P2dBis, P2d);
 535       myBSpline = Handle(Geom2d_BSplineCurve)::DownCast(GG);
 536     }
 537   }
 538 }
 539
 540 //=======================================================================
 541 //function : ProjLib_ComputeApproxOnPolarSurface
 542 //purpose  : case without curve of initialization
 543 //=======================================================================
 544
 545 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 546                       (const Handle(Adaptor3d_HCurve)&   theCurve,
 547                        const Handle(Adaptor3d_HSurface)& theSurface,
 548                        const Standard_Real               theTolerance3D)
 549 : myProjIsDone(Standard_False),
 550   myTolerance (theTolerance3D)
 551 {
 552   const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
 553   myBSpline = Perform(anInitCurve2d, theCurve, theSurface);
 554 }
 555
 556 //=======================================================================
 557 //function : Concat
 558 //purpose  :
 559 //=======================================================================
 560
 561 static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
 562                                           Handle(Geom2d_BSplineCurve) C2,
 563                                           Standard_Real theUJump,
 564                                           Standard_Real theVJump)
 565 {
 566   Standard_Integer deg, deg1, deg2;
 567   deg1 = C1->Degree();
 568   deg2 = C2->Degree();
 569
 570   if ( deg1 < deg2) {
 571     C1->IncreaseDegree(deg2);
 572     deg = deg2;
 573   }
 574   else if ( deg2 < deg1) {
 575     C2->IncreaseDegree(deg1);
 576     deg = deg1;
 577   }
 578   else deg = deg1;
 579
 580   Standard_Integer np1,np2,nk1,nk2,np,nk;
 581   np1 = C1->NbPoles();
 582   nk1 = C1->NbKnots();
 583   np2 = C2->NbPoles();
 584   nk2 = C2->NbKnots();
 585   nk = nk1 + nk2 -1;
 586   np = np1 + np2 -1;
 587
 588   TColStd_Array1OfReal    K1(1,nk1); C1->Knots(K1);
 589   TColStd_Array1OfInteger M1(1,nk1); C1->Multiplicities(M1);
 590   TColgp_Array1OfPnt2d    P1(1,np1); C1->Poles(P1);
 591   TColStd_Array1OfReal    K2(1,nk2); C2->Knots(K2);
 592   TColStd_Array1OfInteger M2(1,nk2); C2->Multiplicities(M2);
 593   TColgp_Array1OfPnt2d    P2(1,np2); C2->Poles(P2);
 594
 595   // Compute the new BSplineCurve
 596   TColStd_Array1OfReal    K(1,nk);
 597   TColStd_Array1OfInteger M(1,nk);
 598   TColgp_Array1OfPnt2d    P(1,np);
 599
 600   Standard_Integer i, count = 0;
 601   // Set Knots and Mults
 602   for ( i = 1; i <= nk1; i++) {
 603     count++;
 604     K(count) = K1(i);
 605     M(count) = M1(i);
 606   }
 607   M(count) = deg;
 608   for ( i = 2; i <= nk2; i++) {
 609     count++;
 610     K(count) = K2(i);
 611     M(count) = M2(i);
 612   }
 613   // Set the Poles
 614   count = 0;
 615   for (i = 1; i <= np1; i++) {
 616     count++;
 617     P(count) = P1(i);
 618   }
 619   for (i = 2; i <= np2; i++) {
 620     count++;
 621     P(count).SetX(P2(i).X() + theUJump);
 622     P(count).SetY(P2(i).Y() + theVJump);
 623   }
 624
 625   Handle(Geom2d_BSplineCurve) BS =
 626     new Geom2d_BSplineCurve(P,K,M,deg);
 627   return BS;
 628 }
 629
 630
 631 //=======================================================================
 632 //function : Perform
 633 //purpose  :
 634 //=======================================================================
 635 Handle(Geom2d_BSplineCurve) ProjLib_ComputeApproxOnPolarSurface::Perform
 636 (const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
 637  const Handle(Adaptor3d_HCurve)& Curve,
 638  const Handle(Adaptor3d_HSurface)& S)
 639 {
 640   //OCC217
 641   Standard_Real Tol3d = myTolerance;
 642   Standard_Real ParamTol = Precision::PApproximation();
 643
 644   Handle(Adaptor2d_HCurve2d) AHC2d = InitialCurve2d;
 645   Handle(Adaptor3d_HCurve) AHC = Curve;
 646
 647 // if the curve 3d is a BSpline with degree C0, it is cut into sections with degree C1
 648 // -> bug cts18237
 649   GeomAbs_CurveType typeCurve = Curve->GetType();
 650   if(typeCurve == GeomAbs_BSplineCurve) {
 651     TColStd_ListOfTransient LOfBSpline2d;
 652     Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
 653     Standard_Integer nbInter = Curve->NbIntervals(GeomAbs_C1);
 654     if(nbInter > 1) {
 655       Standard_Integer i, j;
 656       Handle(Geom_TrimmedCurve) GTC;
 657       Handle(Geom2d_TrimmedCurve) G2dTC;
 658       TColStd_Array1OfReal Inter(1,nbInter+1);
 659       Curve->Intervals(Inter,GeomAbs_C1);
 660       Standard_Real firstinter = Inter.Value(1), secondinter = Inter.Value(2);
 661       // initialization 3d
 662       GTC = new Geom_TrimmedCurve(BSC, firstinter, secondinter);
 663       AHC = new GeomAdaptor_HCurve(GTC);
 664
 665       // if there is an initialization curve:
 666       // - either this is a BSpline C0, with discontinuity at the same parameters of nodes
 667       // and the sections C1 are taken
 668       // - or this is a curve C1 and the sections of intrest are taken otherwise the curve is created.
 669
 670       // initialization 2d
 671       Standard_Integer nbInter2d;
 672       Standard_Boolean C2dIsToCompute;
 673       C2dIsToCompute = InitialCurve2d.IsNull();
 674       Handle(Geom2d_BSplineCurve) BSC2d;
 675       Handle(Geom2d_Curve) G2dC;
 676
 677       if(!C2dIsToCompute) {
 678         nbInter2d = InitialCurve2d->NbIntervals(GeomAbs_C1);
 679         TColStd_Array1OfReal Inter2d(1,nbInter2d+1);
 680         InitialCurve2d->Intervals(Inter2d,GeomAbs_C1);
 681         j = 1;
 682         for(i = 1,j = 1;i <= nbInter;i++)
 683           if(Abs(Inter.Value(i) - Inter2d.Value(j)) < ParamTol) { //OCC217
 684           //if(Abs(Inter.Value(i) - Inter2d.Value(j)) < myTolerance) {
 685             if (j > nbInter2d) break;
 686             j++;
 687           }
 688         if(j != (nbInter2d+1)) {
 689           C2dIsToCompute = Standard_True;
 690         }
 691       }
 692
 693       if(C2dIsToCompute) {
 694         AHC2d = BuildInitialCurve2d(AHC, S);
 695       }
 696       else {
 697         typeCurve = InitialCurve2d->GetType();
 698         switch (typeCurve) {
 699         case GeomAbs_Line: {
 700           G2dC = new Geom2d_Line(InitialCurve2d->Line());
 701           break;
 702         }
 703         case GeomAbs_Circle: {
 704           G2dC = new Geom2d_Circle(InitialCurve2d->Circle());
 705           break;
 706         }
 707         case GeomAbs_Ellipse: {
 708           G2dC = new Geom2d_Ellipse(InitialCurve2d->Ellipse());
 709           break;
 710         }
 711         case GeomAbs_Hyperbola: {
 712           G2dC = new Geom2d_Hyperbola(InitialCurve2d->Hyperbola());
 713           break;
 714         }
 715         case GeomAbs_Parabola: {
 716           G2dC = new Geom2d_Parabola(InitialCurve2d->Parabola());
 717           break;
 718         }
 719         case GeomAbs_BezierCurve: {
 720           G2dC = InitialCurve2d->Bezier();
 721           break;
 722         }
 723         case GeomAbs_BSplineCurve: {
 724           G2dC = InitialCurve2d->BSpline();
 725           break;
 726         }
 727         case GeomAbs_OtherCurve:
 728         default:
 729           break;
 730         }
 731         gp_Pnt2d fp2d = G2dC->Value(firstinter), lp2d = G2dC->Value(secondinter);
 732         gp_Pnt fps, lps, fpc, lpc;
 733         S->D0(fp2d.X(), fp2d.Y(), fps);
 734         S->D0(lp2d.X(), lp2d.Y(), lps);
 735         Curve->D0(firstinter, fpc);
 736         Curve->D0(secondinter, lpc);
 737         //OCC217
 738         if((fps.IsEqual(fpc, Tol3d)) &&
 739            (lps.IsEqual(lpc, Tol3d))) {
 740         //if((fps.IsEqual(fpc, myTolerance)) &&
 741         //   (lps.IsEqual(lpc, myTolerance))) {
 742           G2dTC = new Geom2d_TrimmedCurve(G2dC, firstinter, secondinter);
 743           Geom2dAdaptor_Curve G2dAC(G2dTC);
 744           AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
 745           myProjIsDone = Standard_True;
 746         }
 747         else {
 748           AHC2d = BuildInitialCurve2d(AHC, S);
 749           C2dIsToCompute = Standard_True;
 750         }
 751       }
 752
 753       if(myProjIsDone) {
 754         BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 755         if(BSC2d.IsNull()) return Handle(Geom2d_BSplineCurve)(); //IFV
 756         LOfBSpline2d.Append(BSC2d);
 757       }
 758       else {
 759         return Handle(Geom2d_BSplineCurve)();
 760       }
 761
 762
 763
 764       Standard_Real iinter, ip1inter;
 765       Standard_Integer nbK2d, deg;
 766       nbK2d = BSC2d->NbKnots(); deg = BSC2d->Degree();
 767
 768       for(i = 2;i <= nbInter;i++) {
 769         iinter = Inter.Value(i);
 770         ip1inter = Inter.Value(i+1);
 771         // general case 3d
 772         GTC->SetTrim(iinter, ip1inter);
 773         AHC = new GeomAdaptor_HCurve(GTC);
 774
 775         // general case 2d
 776         if(C2dIsToCompute) {
 777           AHC2d = BuildInitialCurve2d(AHC, S);
 778         }
 779         else {
 780           gp_Pnt2d fp2d = G2dC->Value(iinter), lp2d = G2dC->Value(ip1inter);
 781           gp_Pnt fps, lps, fpc, lpc;
 782           S->D0(fp2d.X(), fp2d.Y(), fps);
 783           S->D0(lp2d.X(), lp2d.Y(), lps);
 784           Curve->D0(iinter, fpc);
 785           Curve->D0(ip1inter, lpc);
 786           //OCC217
 787           if((fps.IsEqual(fpc, Tol3d)) &&
 788              (lps.IsEqual(lpc, Tol3d))) {
 789           //if((fps.IsEqual(fpc, myTolerance)) &&
 790           //   (lps.IsEqual(lpc, myTolerance))) {
 791             G2dTC->SetTrim(iinter, ip1inter);
 792             Geom2dAdaptor_Curve G2dAC(G2dTC);
 793             AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
 794             myProjIsDone = Standard_True;
 795           }
 796           else {
 797             AHC2d = BuildInitialCurve2d(AHC, S);
 798           }
 799         }
 800         if(myProjIsDone) {
 801           BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 802           if(BSC2d.IsNull()) {
 803             return Handle(Geom2d_BSplineCurve)();
 804           }
 805           LOfBSpline2d.Append(BSC2d);
 806           nbK2d += BSC2d->NbKnots() - 1;
 807           deg = Max(deg, BSC2d->Degree());
 808         }
 809         else {
 810           return Handle(Geom2d_BSplineCurve)();
 811         }
 812       }
 813
 814       Standard_Real anUPeriod, anVPeriod;
 815       computePeriodicity(S, anUPeriod, anVPeriod);
 816       Standard_Integer NbC = LOfBSpline2d.Extent();
 817       Handle(Geom2d_BSplineCurve) CurBS;
 818       CurBS = Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
 819       LOfBSpline2d.RemoveFirst();
 820       for (Standard_Integer ii = 2; ii <= NbC; ii++)
 821       {
 822         Handle(Geom2d_BSplineCurve) BS =
 823           Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
 824
 825         //Check for period jump in point of contact.
 826         gp_Pnt2d aC1End = CurBS->Pole(CurBS->NbPoles()); // End of C1.
 827         gp_Pnt2d aC2Beg = BS->Pole(1); // Beginning of C2.
 828         Standard_Real anUJump = 0.0, anVJump = 0.0;
 829
 830         if (anUPeriod > 0.0 &&
 831             Abs (aC1End.X() - aC2Beg.X()) > (anUPeriod ) / 2.01)
 832         {
 833           Standard_Real aMultCoeff =  aC2Beg.X() < aC1End.X() ? 1.0 : -1.0;
 834           anUJump = (anUPeriod) * aMultCoeff;
 835         }
 836
 837         if (anVPeriod &&
 838             Abs (aC1End.Y() - aC2Beg.Y()) > (anVPeriod) / 2.01)
 839         {
 840           Standard_Real aMultCoeff =  aC2Beg.Y() < aC1End.Y() ? 1.0 : -1.0;
 841           anVJump = (anVPeriod) * aMultCoeff;
 842         }
 843
 844         CurBS = Concat(CurBS,BS, anUJump, anVJump);
 845         LOfBSpline2d.RemoveFirst();
 846       }
 847       return CurBS;
 848     }
 849   }
 850
 851   if(InitialCurve2d.IsNull()) {
 852     AHC2d = BuildInitialCurve2d(Curve, S);
 853     if(!myProjIsDone)
 854       return Handle(Geom2d_BSplineCurve)();
 855   }
 856   return ProjectUsingInitialCurve2d(AHC, S, AHC2d);
 857 }
 858
 859 //=======================================================================
 860 //function : ProjLib_BuildInitialCurve2d
 861 //purpose  :
 862 //=======================================================================
 863
 864 Handle(Adaptor2d_HCurve2d)
 865      ProjLib_ComputeApproxOnPolarSurface::
 866      BuildInitialCurve2d(const Handle(Adaptor3d_HCurve)&   Curve,
 867                          const Handle(Adaptor3d_HSurface)& Surf)
 868 {
 869   //  discretize the Curve with quasiuniform deflection
 870   //  density at least NbOfPnts points
 871   myProjIsDone = Standard_False;
 872
 873   //OCC217
 874   Standard_Real Tol3d = myTolerance;
 875   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
 876   Standard_Real DistTol3d = 100.0*Tol3d;
 877
 878   Standard_Real uperiod = 0.0, vperiod = 0.0;
 879   computePeriodicity(Surf, uperiod, vperiod);
 880
 881   // NO myTol is Tol2d !!!!
 882   //Standard_Real TolU = myTolerance, TolV = myTolerance;
 883   //Standard_Real Tol3d = 100*myTolerance; // At random Balthazar.
 884
 885   Standard_Integer NbOfPnts = 61;
 886   GCPnts_QuasiUniformAbscissa QUA(Curve->GetCurve(),NbOfPnts);
 887   TColgp_Array1OfPnt Pts(1,NbOfPnts);
 888   TColStd_Array1OfReal Param(1,NbOfPnts);
 889   Standard_Integer i, j;
 890   for( i = 1; i <= NbOfPnts ; i++ ) {
 891     Param(i) = QUA.Parameter(i);
 892     Pts(i) = Curve->Value(Param(i));
 893   }
 894
 895   TColgp_Array1OfPnt2d Pts2d(1,NbOfPnts);
 896   TColStd_Array1OfInteger Mult(1,NbOfPnts);
 897   Mult.Init(1);
 898   Mult(1) = Mult(NbOfPnts) = 2;
 899
 900   Standard_Real Uinf, Usup, Vinf, Vsup;
 901   Uinf = Surf->Surface().FirstUParameter();
 902   Usup = Surf->Surface().LastUParameter();
 903   Vinf = Surf->Surface().FirstVParameter();
 904   Vsup = Surf->Surface().LastVParameter();
 905   GeomAbs_SurfaceType Type = Surf->GetType();
 906   if((Type != GeomAbs_BSplineSurface) && (Type != GeomAbs_BezierSurface) &&
 907      (Type != GeomAbs_OffsetSurface)) {
 908     Standard_Real S, T;
 909 //    Standard_Integer usens = 0, vsens = 0;
 910     // to know the position relatively to the period
 911     switch (Type) {
 912 //    case GeomAbs_Plane:
 913 //      {
 914 //      gp_Pln Plane = Surf->Plane();
 915 //      for ( i = 1 ; i <= NbOfPnts ; i++) {
 916 //        ElSLib::Parameters( Plane, Pts(i), S, T);
 917 //        Pts2d(i).SetCoord(S,T);
 918 //      }
 919 //      myProjIsDone = Standard_True;
 920 //      break;
 921 //      }
 922     case GeomAbs_Cylinder:
 923       {
 924 //      Standard_Real Sloc, Tloc;
 925         Standard_Real Sloc;
 926         Standard_Integer usens = 0;
 927         gp_Cylinder Cylinder = Surf->Cylinder();
 928         ElSLib::Parameters( Cylinder, Pts(1), S, T);
 929         Pts2d(1).SetCoord(S,T);
 930         for ( i = 2 ; i <= NbOfPnts ; i++) {
 931           Sloc = S;
 932           ElSLib::Parameters( Cylinder, Pts(i), S, T);
 933           if(Abs(Sloc - S) > M_PI) {
 934             if(Sloc > S)
 935               usens++;
 936             else
 937               usens--;
 938           }
 939           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
 940         }
 941         myProjIsDone = Standard_True;
 942         break;
 943       }
 944     case GeomAbs_Cone:
 945       {
 946 //      Standard_Real Sloc, Tloc;
 947         Standard_Real Sloc;
 948         Standard_Integer usens = 0;
 949         gp_Cone Cone = Surf->Cone();
 950         ElSLib::Parameters( Cone, Pts(1), S, T);
 951         Pts2d(1).SetCoord(S,T);
 952         for ( i = 2 ; i <= NbOfPnts ; i++) {
 953           Sloc = S;
 954           ElSLib::Parameters( Cone, Pts(i), S, T);
 955           if(Abs(Sloc - S) > M_PI) {
 956             if(Sloc > S)
 957               usens++;
 958             else
 959               usens--;
 960           }
 961           Pts2d(i).SetCoord(S+usens*2*M_PI,T);
 962         }
 963         myProjIsDone = Standard_True;
 964         break;
 965       }
 966     case GeomAbs_Sphere:
 967       {
 968         Standard_Real Sloc, Tloc;
 969         Standard_Integer usens = 0, vsens = 0; //usens steps by half-period
 970         Standard_Boolean vparit = Standard_False;
 971         gp_Sphere Sphere = Surf->Sphere();
 972         ElSLib::Parameters( Sphere, Pts(1), S, T);
 973         Pts2d(1).SetCoord(S,T);
 974         for ( i = 2 ; i <= NbOfPnts ; i++) {
 975           Sloc = S;Tloc = T;
 976           ElSLib::Parameters( Sphere, Pts(i), S, T);
 977           if(1.6*M_PI < Abs(Sloc - S)) {
 978             if(Sloc > S)
 979               usens += 2;
 980             else
 981               usens -= 2;
 982     }
 983           if(1.6*M_PI > Abs(Sloc - S) && Abs(Sloc - S) > 0.4*M_PI) {
 984             vparit = !vparit;
 985             if(Sloc > S)
 986               usens++;
 987             else
 988               usens--;
 989             if(Abs(Tloc - Vsup) < (Vsup - Vinf)/5)
 990               vsens++;
 991             else
 992               vsens--;
 993           }
 994           if(vparit) {
 995             Pts2d(i).SetCoord(S+usens*M_PI,(M_PI - T)*(vsens-1));
 996           }
 997           else {
 998             Pts2d(i).SetCoord(S+usens*M_PI,T+vsens*M_PI);
 999
1000           }
1001         }
1002         myProjIsDone = Standard_True;
1003         break;
1004       }
1005     case GeomAbs_Torus:
1006       {
1007         Standard_Real Sloc, Tloc;
1008         Standard_Integer usens = 0, vsens = 0;
1009         gp_Torus Torus = Surf->Torus();
1010         ElSLib::Parameters( Torus, Pts(1), S, T);
1011         Pts2d(1).SetCoord(S,T);
1012         for ( i = 2 ; i <= NbOfPnts ; i++) {
1013           Sloc = S; Tloc = T;
1014           ElSLib::Parameters( Torus, Pts(i), S, T);
1015           if(Abs(Sloc - S) > M_PI) {
1016             if(Sloc > S)
1017               usens++;
1018             else
1019               usens--;
1020     }
1021           if(Abs(Tloc - T) > M_PI) {
1022             if(Tloc > T)
1023               vsens++;
1024             else
1025               vsens--;
1026     }
1027           Pts2d(i).SetCoord(S+usens*2*M_PI,T+vsens*2*M_PI);
1028         }
1029         myProjIsDone = Standard_True;
1030         break;
1031       }
1032     default:
1033       Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::BuildInitialCurve2d");
1034     }
1035   }
1036   else {
1037     myProjIsDone = Standard_False;
1038     Standard_Real Dist2Min = 1.e+200, u = 0., v = 0.;
1039     gp_Pnt pntproj;
1040
1041     TColgp_SequenceOfPnt2d Sols;
1042     Standard_Boolean areManyZeros = Standard_False;
1043
1044     Curve->D0(Param.Value(1), pntproj) ;
1045     Extrema_ExtPS  aExtPS(pntproj, Surf->Surface(), TolU, TolV) ;
1046     Standard_Real aMinSqDist = RealLast();
1047     if (aExtPS.IsDone())
1048     {
1049       for (i = 1; i <= aExtPS.NbExt(); i++)
1050       {
1051         Standard_Real aSqDist = aExtPS.SquareDistance(i);
1052         if (aSqDist < aMinSqDist)
1053           aMinSqDist = aSqDist;
1054       }
1055     }
1056     if (aMinSqDist > DistTol3d * DistTol3d) //try to project with less tolerance
1057     {
1058       TolU = Min(TolU, Precision::PConfusion());
1059       TolV = Min(TolV, Precision::PConfusion());
1060       aExtPS.Initialize(Surf->Surface(),
1061                         Surf->Surface().FirstUParameter(), Surf->Surface().LastUParameter(),
1062                         Surf->Surface().FirstVParameter(), Surf->Surface().LastVParameter(),
1063                         TolU, TolV);
1064       aExtPS.Perform(pntproj);
1065     }
1066
1067     if( aExtPS.IsDone() && aExtPS.NbExt() >= 1 ) {
1068
1069       Standard_Integer GoodValue = 1;
1070
1071       for ( i = 1 ; i <= aExtPS.NbExt() ; i++ ) {
1072         if( aExtPS.SquareDistance(i) < DistTol3d * DistTol3d ) {
1073           if( aExtPS.SquareDistance(i) <= 1.e-18 ) {
1074             aExtPS.Point(i).Parameter(u,v);
1075             gp_Pnt2d p2d(u,v);
1076             Standard_Boolean isSame = Standard_False;
1077             for( j = 1; j <= Sols.Length(); j++ ) {
1078               if( p2d.SquareDistance( Sols.Value(j) ) <= 1.e-18 ) {
1079                 isSame = Standard_True;
1080                 break;
1081               }
1082             }
1083             if( !isSame ) Sols.Append( p2d );
1084           }
1085           if( Dist2Min > aExtPS.SquareDistance(i) ) {
1086             Dist2Min = aExtPS.SquareDistance(i);
1087             GoodValue = i;
1088           }
1089         }
1090       }
1091
1092       if( Sols.Length() > 1 ) areManyZeros = Standard_True;
1093
1094       if( Dist2Min <= DistTol3d * DistTol3d) {
1095         if( !areManyZeros ) {
1096           aExtPS.Point(GoodValue).Parameter(u,v);
1097           Pts2d(1).SetCoord(u,v);
1098           myProjIsDone = Standard_True;
1099         }
1100         else {
1101           Standard_Integer nbSols = Sols.Length();
1102           Standard_Real Dist2Max = -1.e+200;
1103           for( i = 1; i <= nbSols; i++ ) {
1104             const gp_Pnt2d& aP1 = Sols.Value(i);
1105             for( j = i+1; j <= nbSols; j++ ) {
1106               const gp_Pnt2d& aP2 = Sols.Value(j);
1107               Standard_Real aDist2 = aP1.SquareDistance(aP2);
1108               if( aDist2 > Dist2Max ) Dist2Max = aDist2;
1109             }
1110           }
1111       Standard_Real aMaxT2 = Max(TolU,TolV);
1112       aMaxT2 *= aMaxT2;
1113           if( Dist2Max > aMaxT2 ) {
1114             Standard_Integer tPp = 0;
1115             for( i = 1; i <= 5; i++ ) {
1116               Standard_Integer nbExtOk = 0;
1117               Standard_Integer indExt = 0;
1118               Standard_Integer iT = 1 + (NbOfPnts - 1)/5*i;
1119               Curve->D0( Param.Value(iT), pntproj );
1120               Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1121               Dist2Min = 1.e+200;
1122               if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1123                 for( j = 1 ; j <= aTPS.NbExt() ; j++ ) {
1124                   if( aTPS.SquareDistance(j) < DistTol3d * DistTol3d ) {
1125                     nbExtOk++;
1126                     if( aTPS.SquareDistance(j) < Dist2Min ) {
1127                       Dist2Min = aTPS.SquareDistance(j);
1128                       indExt = j;
1129                     }
1130                   }
1131                 }
1132               }
1133               if( nbExtOk == 1 ) {
1134                 tPp = iT;
1135                 aTPS.Point(indExt).Parameter(u,v);
1136                 break;
1137               }
1138             }
1139
1140             if( tPp != 0 ) {
1141               gp_Pnt2d aPp = gp_Pnt2d(u,v);
1142               gp_Pnt2d aPn;
1143               j = 1;
1144               Standard_Boolean isFound = Standard_False;
1145               while( !isFound ) {
1146                 Curve->D0( Param.Value(tPp+j), pntproj );
1147                 Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1148                 Dist2Min = 1.e+200;
1149                 Standard_Integer indExt = 0;
1150                 if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1151                   for( i = 1 ; i <= aTPS.NbExt() ; i++ ) {
1152                     if( aTPS.SquareDistance(i) < DistTol3d * DistTol3d && aTPS.SquareDistance(i) < Dist2Min ) {
1153                       Dist2Min = aTPS.SquareDistance(i);
1154                       indExt = i;
1155                       isFound = Standard_True;
1156                     }
1157                   }
1158                 }
1159                 if( isFound ) {
1160                   aTPS.Point(indExt).Parameter(u,v);
1161                   aPn = gp_Pnt2d(u,v);
1162                   break;
1163                 }
1164                 j++;
1165                 if( (tPp+j) > NbOfPnts ) break;
1166               }
1167
1168               if( isFound ) {
1169                 gp_Vec2d atV(aPp,aPn);
1170                 Standard_Boolean isChosen = Standard_False;
1171                 for( i = 1; i <= nbSols; i++ ) {
1172                   const gp_Pnt2d& aP1 = Sols.Value(i);
1173                   gp_Vec2d asV(aP1,aPp);
1174                   if( asV.Dot(atV) > 0. ) {
1175                     isChosen = Standard_True;
1176                     Pts2d(1).SetCoord(aP1.X(),aP1.Y());
1177                     myProjIsDone = Standard_True;
1178                     break;
1179                   }
1180                 }
1181                 if( !isChosen ) {
1182                   aExtPS.Point(GoodValue).Parameter(u,v);
1183                   Pts2d(1).SetCoord(u,v);
1184                   myProjIsDone = Standard_True;
1185                 }
1186               }
1187               else {
1188                 aExtPS.Point(GoodValue).Parameter(u,v);
1189                 Pts2d(1).SetCoord(u,v);
1190                 myProjIsDone = Standard_True;
1191               }
1192             }
1193             else {
1194               aExtPS.Point(GoodValue).Parameter(u,v);
1195               Pts2d(1).SetCoord(u,v);
1196               myProjIsDone = Standard_True;
1197             }
1198           }
1199           else {
1200             aExtPS.Point(GoodValue).Parameter(u,v);
1201             Pts2d(1).SetCoord(u,v);
1202             myProjIsDone = Standard_True;
1203           }
1204         }
1205       }
1206
1207       //  calculate the following points with GenLocate_ExtPS
1208       // (and store the result and each parameter in a sequence)
1209       Standard_Integer usens = 0, vsens = 0;
1210       // to know the position relatively to the period
1211       Standard_Real U0 = u, V0 = v, U1 = u, V1 = v;
1212       // U0 and V0 are the points in the initialized period
1213       // (period with u and v),
1214       // U1 and V1 are the points for construction of poles
1215
1216       for ( i = 2 ; i <= NbOfPnts ; i++)
1217         if(myProjIsDone) {
1218           myProjIsDone = Standard_False;
1219           Dist2Min = RealLast();
1220           Curve->D0(Param.Value(i), pntproj);
1221           Extrema_GenLocateExtPS  aLocateExtPS
1222             (pntproj, Surf->Surface(), U0, V0, TolU, TolV) ;
1223
1224           if (aLocateExtPS.IsDone())
1225           {
1226             if (aLocateExtPS.SquareDistance() < DistTol3d * DistTol3d)
1227             {  //OCC217
1228               //if (aLocateExtPS.SquareDistance() < Tol3d * Tol3d) {
1229               (aLocateExtPS.Point()).Parameter(U0,V0);
1230               U1 = U0 + usens*uperiod;
1231               V1 = V0 + vsens*vperiod;
1232               Pts2d(i).SetCoord(U1,V1);
1233               myProjIsDone = Standard_True;
1234             }
1235             else
1236             {
1237               Extrema_ExtPS aGlobalExtr(pntproj, Surf->Surface(), TolU, TolV);
1238               if (aGlobalExtr.IsDone())
1239               {
1240                 Standard_Real LocalMinSqDist = RealLast();
1241                 Standard_Integer imin = 0;
1242                 for (Standard_Integer isol = 1; isol <= aGlobalExtr.NbExt(); isol++)
1243                 {
1244                   Standard_Real aSqDist = aGlobalExtr.SquareDistance(isol);
1245                   if (aSqDist < LocalMinSqDist)
1246                   {
1247                     LocalMinSqDist = aSqDist;
1248                     imin = isol;
1249                   }
1250                 }
1251                 if (LocalMinSqDist < DistTol3d * DistTol3d)
1252                 {
1253                   Standard_Real LocalU, LocalV;
1254                   aGlobalExtr.Point(imin).Parameter(LocalU, LocalV);
1255                   if (uperiod > 0. && Abs(U0 - LocalU) >= uperiod/2.)
1256                   {
1257                     if (LocalU > U0)
1258                       usens = -1;
1259                     else
1260                       usens = 1;
1261                   }
1262                   if (vperiod > 0. && Abs(V0 - LocalV) >= vperiod/2.)
1263                   {
1264                     if (LocalV > V0)
1265                       vsens = -1;
1266                     else
1267                       vsens = 1;
1268                   }
1269                   U0 = LocalU; V0 = LocalV;
1270                   U1 = U0 + usens*uperiod;
1271                   V1 = V0 + vsens*vperiod;
1272                   Pts2d(i).SetCoord(U1,V1);
1273                   myProjIsDone = Standard_True;
1274
1275                   if((i == 2) && (!IsEqual(uperiod, 0.0) || !IsEqual(vperiod, 0.0)))
1276                   {//Make 1st point more precise for periodic surfaces
1277                     const Standard_Integer aSize = 3;
1278                     const gp_Pnt2d aP(Pts2d(2));
1279                     Standard_Real aUpar[aSize], aVpar[aSize];
1280                     Pts2d(1).Coord(aUpar[1], aVpar[1]);
1281                     aUpar[0] = aUpar[1] - uperiod;
1282                     aUpar[2] = aUpar[1] + uperiod;
1283                     aVpar[0] = aVpar[1] - vperiod;
1284                     aVpar[2] = aVpar[1] + vperiod;
1285
1286                     Standard_Real aSQdistMin = RealLast();
1287                     Standard_Integer aBestUInd = 1, aBestVInd = 1;
1288                     const Standard_Integer  aSizeU = IsEqual(uperiod, 0.0) ? 1 : aSize,
1289                                             aSizeV = IsEqual(vperiod, 0.0) ? 1 : aSize;
1290                     for(Standard_Integer uInd = 0; uInd < aSizeU; uInd++)
1291                     {
1292                       for(Standard_Integer vInd = 0; vInd < aSizeV; vInd++)
1293                       {
1294                         Standard_Real aSQdist = aP.SquareDistance(gp_Pnt2d(aUpar[uInd], aVpar[vInd]));
1295                         if(aSQdist < aSQdistMin)
1296                         {
1297                           aSQdistMin = aSQdist;
1298                           aBestUInd = uInd;
1299                           aBestVInd = vInd;
1300                         }
1301                       }
1302                     }
1303
1304                     Pts2d(1).SetCoord(aUpar[aBestUInd], aVpar[aBestVInd]);
1305                   }//if(i == 2) condition
1306                 }
1307               }
1308             }
1309           }
1310           if(!myProjIsDone && uperiod) {
1311             Standard_Real aUinf, aUsup, Uaux;
1312             aUinf = Surf->Surface().FirstUParameter();
1313             aUsup = Surf->Surface().LastUParameter();
1314             if((aUsup - U0) > (U0 - aUinf))
1315               Uaux = 2*aUinf - U0 + uperiod;
1316             else
1317               Uaux = 2*aUsup - U0 - uperiod;
1318             Extrema_GenLocateExtPS  locext(pntproj,
1319                                            Surf->Surface(),
1320                                            Uaux, V0, TolU, TolV);
1321             if (locext.IsDone())
1322               if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
1323               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1324                 (locext.Point()).Parameter(u,v);
1325                 if((aUsup - U0) > (U0 - aUinf))
1326                   usens--;
1327                 else
1328                   usens++;
1329                 U0 = u; V0 = v;
1330                 U1 = U0 + usens*uperiod;
1331                 V1 = V0 + vsens*vperiod;
1332                 Pts2d(i).SetCoord(U1,V1);
1333                 myProjIsDone = Standard_True;
1334               }
1335           }
1336           if(!myProjIsDone && vperiod) {
1337             Standard_Real aVinf, aVsup, Vaux;
1338             aVinf = Surf->Surface().FirstVParameter();
1339             aVsup = Surf->Surface().LastVParameter();
1340             if((aVsup - V0) > (V0 - aVinf))
1341               Vaux = 2*aVinf - V0 + vperiod;
1342             else
1343               Vaux = 2*aVsup - V0 - vperiod;
1344             Extrema_GenLocateExtPS  locext(pntproj,
1345                                            Surf->Surface(),
1346                                            U0, Vaux, TolU, TolV) ;
1347             if (locext.IsDone())
1348               if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
1349               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1350                 (locext.Point()).Parameter(u,v);
1351                 if((aVsup - V0) > (V0 - aVinf))
1352                   vsens--;
1353                 else
1354                   vsens++;
1355                 U0 = u; V0 = v;
1356                 U1 = U0 + usens*uperiod;
1357                 V1 = V0 + vsens*vperiod;
1358                 Pts2d(i).SetCoord(U1,V1);
1359                 myProjIsDone = Standard_True;
1360               }
1361           }
1362           if(!myProjIsDone && uperiod && vperiod) {
1363             Standard_Real Uaux, Vaux;
1364             if((Usup - U0) > (U0 - Uinf))
1365               Uaux = 2*Uinf - U0 + uperiod;
1366             else
1367               Uaux = 2*Usup - U0 - uperiod;
1368             if((Vsup - V0) > (V0 - Vinf))
1369               Vaux = 2*Vinf - V0 + vperiod;
1370             else
1371               Vaux = 2*Vsup - V0 - vperiod;
1372             Extrema_GenLocateExtPS  locext(pntproj,
1373                                            Surf->Surface(),
1374                                            Uaux, Vaux, TolU, TolV);
1375             if (locext.IsDone())
1376               if (locext.SquareDistance() < DistTol3d * DistTol3d) {
1377               //if (locext.SquareDistance() < Tol3d * Tol3d) {
1378                 (locext.Point()).Parameter(u,v);
1379                 if((Usup - U0) > (U0 - Uinf))
1380                   usens--;
1381                 else
1382                   usens++;
1383                 if((Vsup - V0) > (V0 - Vinf))
1384                   vsens--;
1385                 else
1386                   vsens++;
1387                 U0 = u; V0 = v;
1388                 U1 = U0 + usens*uperiod;
1389                 V1 = V0 + vsens*vperiod;
1390                 Pts2d(i).SetCoord(U1,V1);
1391                 myProjIsDone = Standard_True;
1392               }
1393           }
1394           if(!myProjIsDone) {
1395             Extrema_ExtPS ext(pntproj, Surf->Surface(), TolU, TolV) ;
1396             if (ext.IsDone()) {
1397               Dist2Min = ext.SquareDistance(1);
1398               Standard_Integer aGoodValue = 1;
1399               for ( j = 2 ; j <= ext.NbExt() ; j++ )
1400                 if( Dist2Min > ext.SquareDistance(j)) {
1401                   Dist2Min = ext.SquareDistance(j);
1402                   aGoodValue = j;
1403                 }
1404               if (Dist2Min < DistTol3d * DistTol3d) {
1405               //if (Dist2Min < Tol3d * Tol3d) {
1406                 (ext.Point(aGoodValue)).Parameter(u,v);
1407                 if(uperiod) {
1408                   if((U0 - u) > (2*uperiod/3)) {
1409                     usens++;
1410                   }
1411                   else
1412                     if((u - U0) > (2*uperiod/3)) {
1413                       usens--;
1414                     }
1415     }
1416                 if(vperiod) {
1417                   if((V0 - v) > (vperiod/2)) {
1418                     vsens++;
1419                   }
1420                   else
1421                     if((v - V0) > (vperiod/2)) {
1422                       vsens--;
1423                     }
1424       }
1425                 U0 = u; V0 = v;
1426                 U1 = U0 + usens*uperiod;
1427                 V1 = V0 + vsens*vperiod;
1428                 Pts2d(i).SetCoord(U1,V1);
1429                 myProjIsDone = Standard_True;
1430               }
1431             }
1432           }
1433         }
1434         else break;
1435     }
1436   }
1437   // -- Pnts2d is transformed into Geom2d_BSplineCurve, with the help of Param and Mult
1438   if(myProjIsDone) {
1439     myBSpline = new Geom2d_BSplineCurve(Pts2d,Param,Mult,1);
1440     //jgv: put the curve into parametric range
1441     gp_Pnt2d MidPoint = myBSpline->Value(0.5*(myBSpline->FirstParameter() + myBSpline->LastParameter()));
1442     Standard_Real TestU = MidPoint.X(), TestV = MidPoint.Y();
1443     Standard_Real sense = 0.;
1444     if (uperiod)
1445     {
1446       if (TestU < Uinf - TolU)
1447         sense = 1.;
1448       else if (TestU > Usup + TolU)
1449         sense = -1;
1450       while (TestU < Uinf - TolU || TestU > Usup + TolU)
1451         TestU += sense * uperiod;
1452     }
1453     if (vperiod)
1454     {
1455       sense = 0.;
1456       if (TestV < Vinf - TolV)
1457         sense = 1.;
1458       else if (TestV > Vsup + TolV)
1459         sense = -1.;
1460       while (TestV < Vinf - TolV || TestV > Vsup + TolV)
1461         TestV += sense * vperiod;
1462     }
1463     gp_Vec2d Offset(TestU - MidPoint.X(), TestV - MidPoint.Y());
1464     if (Abs(Offset.X()) > gp::Resolution() ||
1465         Abs(Offset.Y()) > gp::Resolution())
1466       myBSpline->Translate(Offset);
1467     //////////////////////////////////////////
1468     Geom2dAdaptor_Curve GAC(myBSpline);
1469     Handle(Adaptor2d_HCurve2d) IC2d = new Geom2dAdaptor_HCurve(GAC);
1470 #ifdef OCCT_DEBUG
1471 //    char name [100];
1472 //    sprintf(name,"%s_%d","build",compteur++);
1473 //    DrawTrSurf::Set(name,myBSpline);
1474 #endif
1475     return IC2d;
1476   }
1477   else {
1478 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:50 2002 Begin
1479 //     Standard_NoSuchObject_Raise_if(1,"ProjLib_Compu: build echec");
1480 //  Modified by Sergey KHROMOV - Thu Apr 18 10:57:51 2002 End
1481     return Handle(Adaptor2d_HCurve2d)();
1482   }
1483 //  myProjIsDone = Standard_False;
1484 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:01 2002 Begin
1485 //   Standard_NoSuchObject_Raise_if(1,"ProjLib_ComputeOnPS: build echec");
1486 //  Modified by Sergey KHROMOV - Thu Apr 18 10:58:02 2002 End
1487 }
1488
1489
1490
1491
1492 //=======================================================================
1493 //function : ProjLib_ProjectUsingInitialCurve2d
1494 //purpose  :
1495 //=======================================================================
1496 Handle(Geom2d_BSplineCurve)
1497      ProjLib_ComputeApproxOnPolarSurface::
1498      ProjectUsingInitialCurve2d(const Handle(Adaptor3d_HCurve)& Curve,
1499                                 const Handle(Adaptor3d_HSurface)& Surf,
1500                                 const Handle(Adaptor2d_HCurve2d)& InitCurve2d)
1501 {
1502   //OCC217
1503   Standard_Real Tol3d = myTolerance;
1504   Standard_Real DistTol3d = 1.0*Tol3d;
1505   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
1506   Standard_Real Tol2d = Max(Sqrt(TolU*TolU + TolV*TolV), Precision::PConfusion());
1507
1508   Standard_Integer i;
1509   GeomAbs_SurfaceType TheTypeS = Surf->GetType();
1510   GeomAbs_CurveType TheTypeC = Curve->GetType();
1511 //  Handle(Standard_Type) TheTypeS = Surf->DynamicType();
1512 //  Handle(Standard_Type) TheTypeC = Curve->DynamicType();   // si on a :
1513 //  if(TheTypeS == STANDARD_TYPE(Geom_BSplineSurface)) {
1514   if(TheTypeS == GeomAbs_Plane) {
1515     Standard_Real S, T;
1516     gp_Pln Plane = Surf->Plane();
1517     if(TheTypeC == GeomAbs_BSplineCurve) {
1518       Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1519       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1520       for(i = 1;i <= Curve->NbPoles();i++) {
1521         ElSLib::Parameters( Plane, BSC->Pole(i), S, T);
1522         Poles2d(i).SetCoord(S,T);
1523       }
1524       TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1525       BSC->Knots(Knots);
1526       TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1527       BSC->Multiplicities(Mults);
1528       if(BSC->IsRational()) {
1529         TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1530         BSC->Weights(Weights);
1531         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1532                                        BSC->Degree(), BSC->IsPeriodic()) ;
1533       }
1534       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1535                                      BSC->Degree(), BSC->IsPeriodic()) ;
1536
1537     }
1538     if(TheTypeC == GeomAbs_BezierCurve) {
1539       Handle(Geom_BezierCurve) BC = Curve->Bezier();
1540       TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1541       for(i = 1;i <= Curve->NbPoles();i++) {
1542         ElSLib::Parameters( Plane, BC->Pole(i), S, T);
1543         Poles2d(i).SetCoord(S,T);
1544       }
1545       TColStd_Array1OfReal Knots(1, 2);
1546       Knots.SetValue(1,0.0);
1547       Knots.SetValue(2,1.0);
1548       TColStd_Array1OfInteger Mults(1, 2);
1549       Mults.Init(BC->NbPoles());
1550       if(BC->IsRational()) {
1551         TColStd_Array1OfReal Weights(1, BC->NbPoles());
1552         BC->Weights(Weights);
1553         return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1554                                        BC->Degree(), BC->IsPeriodic()) ;
1555       }
1556       return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1557                                      BC->Degree(), BC->IsPeriodic()) ;
1558     }
1559   }
1560   if(TheTypeS == GeomAbs_BSplineSurface) {
1561     Handle(Geom_BSplineSurface) BSS = Surf->BSpline();
1562     if((BSS->MaxDegree() == 1) &&
1563        (BSS->NbUPoles() == 2) &&
1564        (BSS->NbVPoles() == 2)) {
1565       gp_Pnt p11 = BSS->Pole(1,1);
1566       gp_Pnt p12 = BSS->Pole(1,2);
1567       gp_Pnt p21 = BSS->Pole(2,1);
1568       gp_Pnt p22 = BSS->Pole(2,2);
1569       gp_Vec V1(p11,p12);
1570       gp_Vec V2(p21,p22);
1571       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){  //OCC217
1572       //if(V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){
1573         //  so the polar surface is plane
1574         //  and if it is enough to projet the  poles of Curve
1575         Standard_Integer Dist2Min = IntegerLast();
1576         Standard_Real u,v;
1577         //OCC217
1578         //Standard_Real TolU = Surf->UResolution(myTolerance)
1579         //  , TolV = Surf->VResolution(myTolerance);
1580 //      gp_Pnt pntproj;
1581         if(TheTypeC == GeomAbs_BSplineCurve) {
1582           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1583           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1584           for(i = 1;i <= Curve->NbPoles();i++) {
1585             myProjIsDone = Standard_False;
1586             Dist2Min = IntegerLast();
1587             Extrema_GenLocateExtPS  extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1588                                             (p11.Y()+p22.Y())/2,TolU,TolV) ;
1589             if (extrloc.IsDone()) {
1590               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1591               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1592               //if (Dist2Min < myTolerance * myTolerance) {
1593                 (extrloc.Point()).Parameter(u,v);
1594                 Poles2d(i).SetCoord(u,v);
1595                 myProjIsDone = Standard_True;
1596               }
1597               else break;
1598             }
1599             else break;
1600             if(!myProjIsDone)
1601               break;
1602           }
1603           if(myProjIsDone) {
1604             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1605             BSC->Knots(Knots);
1606             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1607             BSC->Multiplicities(Mults);
1608             if(BSC->IsRational()) {
1609               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1610               BSC->Weights(Weights);
1611               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1612                                              BSC->Degree(), BSC->IsPeriodic()) ;
1613             }
1614             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1615                                            BSC->Degree(), BSC->IsPeriodic()) ;
1616
1617
1618           }
1619         }
1620         if(TheTypeC == GeomAbs_BezierCurve) {
1621           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1622           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1623           for(i = 1;i <= Curve->NbPoles();i++) {
1624             Dist2Min = IntegerLast();
1625             Extrema_GenLocateExtPS  extrloc(BC->Pole(i),Surf->Surface(),0.5,
1626                                             0.5,TolU,TolV) ;
1627             if (extrloc.IsDone()) {
1628               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1629               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1630               //if (Dist2Min < myTolerance * myTolerance) {
1631                 (extrloc.Point()).Parameter(u,v);
1632                 Poles2d(i).SetCoord(u,v);
1633                 myProjIsDone = Standard_True;
1634               }
1635               else break;
1636             }
1637             else break;
1638             if(myProjIsDone)
1639               myProjIsDone = Standard_False;
1640             else break;
1641           }
1642           if(myProjIsDone) {
1643             TColStd_Array1OfReal Knots(1, 2);
1644             Knots.SetValue(1,0.0);
1645             Knots.SetValue(2,1.0);
1646             TColStd_Array1OfInteger Mults(1, 2);
1647             Mults.Init(BC->NbPoles());
1648             if(BC->IsRational()) {
1649               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1650               BC->Weights(Weights);
1651               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1652                                                     BC->Degree(), BC->IsPeriodic()) ;
1653             }
1654             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1655                                                   BC->Degree(), BC->IsPeriodic()) ;
1656           }
1657         }
1658       }
1659     }
1660   }
1661   else if(TheTypeS == GeomAbs_BezierSurface) {
1662     Handle(Geom_BezierSurface) BS = Surf->Bezier();
1663     if((BS->MaxDegree() == 1) &&
1664        (BS->NbUPoles() == 2) &&
1665        (BS->NbVPoles() == 2)) {
1666       gp_Pnt p11 = BS->Pole(1,1);
1667       gp_Pnt p12 = BS->Pole(1,2);
1668       gp_Pnt p21 = BS->Pole(2,1);
1669       gp_Pnt p22 = BS->Pole(2,2);
1670       gp_Vec V1(p11,p12);
1671       gp_Vec V2(p21,p22);
1672       if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){ //OCC217
1673       //if (V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){
1674         //    and if it is enough to project the poles of Curve
1675         Standard_Integer Dist2Min = IntegerLast();
1676         Standard_Real u,v;
1677         //OCC217
1678         //Standard_Real TolU = Surf->UResolution(myTolerance)
1679         //  , TolV = Surf->VResolution(myTolerance);
1680
1681 //      gp_Pnt pntproj;
1682         if(TheTypeC == GeomAbs_BSplineCurve) {
1683           Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1684           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1685           for(i = 1;i <= Curve->NbPoles();i++) {
1686             myProjIsDone = Standard_False;
1687             Dist2Min = IntegerLast();
1688             Extrema_GenLocateExtPS  extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1689                                             (p11.Y()+p22.Y())/2,TolU,TolV) ;
1690             if (extrloc.IsDone()) {
1691               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1692               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1693               //if (Dist2Min < myTolerance * myTolerance) {
1694                 (extrloc.Point()).Parameter(u,v);
1695                 Poles2d(i).SetCoord(u,v);
1696                 myProjIsDone = Standard_True;
1697               }
1698               else break;
1699             }
1700             else break;
1701             if(!myProjIsDone)
1702               break;
1703           }
1704           if(myProjIsDone) {
1705             TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1706             BSC->Knots(Knots);
1707             TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1708             BSC->Multiplicities(Mults);
1709             if(BSC->IsRational()) {
1710               TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1711               BSC->Weights(Weights);
1712               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1713                                                     BSC->Degree(), BSC->IsPeriodic()) ;
1714             }
1715             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1716                                                   BSC->Degree(), BSC->IsPeriodic()) ;
1717
1718
1719           }
1720         }
1721         if(TheTypeC == GeomAbs_BezierCurve) {
1722           Handle(Geom_BezierCurve) BC = Curve->Bezier();
1723           TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1724           for(i = 1;i <= Curve->NbPoles();i++) {
1725             Dist2Min = IntegerLast();
1726             Extrema_GenLocateExtPS  extrloc(BC->Pole(i),Surf->Surface(),0.5,
1727                                             0.5,TolU,TolV) ;
1728             if (extrloc.IsDone()) {
1729               Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1730               if (Dist2Min < DistTol3d * DistTol3d) {  //OCC217
1731               //if (Dist2Min < myTolerance * myTolerance) {
1732                 (extrloc.Point()).Parameter(u,v);
1733                 Poles2d(i).SetCoord(u,v);
1734                 myProjIsDone = Standard_True;
1735               }
1736               else break;
1737             }
1738             else break;
1739             if(myProjIsDone)
1740               myProjIsDone = Standard_False;
1741             else break;
1742           }
1743           if(myProjIsDone) {
1744             TColStd_Array1OfReal Knots(1, 2);
1745             Knots.SetValue(1,0.0);
1746             Knots.SetValue(2,1.0);
1747             TColStd_Array1OfInteger Mults(1, 2);
1748             Mults.Init(BC->NbPoles());
1749             if(BC->IsRational()) {
1750               TColStd_Array1OfReal Weights(1, BC->NbPoles());
1751               BC->Weights(Weights);
1752               return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1753                                                     BC->Degree(), BC->IsPeriodic()) ;
1754             }
1755             return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1756                                                   BC->Degree(), BC->IsPeriodic()) ;
1757           }
1758         }
1759       }
1760     }
1761   }
1762
1763   ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, Tol3d) ;  //OCC217
1764   //ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, myTolerance) ;
1765
1766 #ifdef OCCT_DEBUG
1767   Standard_Integer Nb = 50;
1768
1769   Standard_Real U, U1, U2;
1770   U1 = F.FirstParameter();
1771   U2 = F.LastParameter();
1772
1773   TColgp_Array1OfPnt2d    DummyPoles(1,Nb+1);
1774   TColStd_Array1OfReal    DummyKnots(1,Nb+1);
1775   TColStd_Array1OfInteger DummyMults(1,Nb+1);
1776   DummyMults.Init(1);
1777   DummyMults(1) = 2;
1778   DummyMults(Nb+1) = 2;
1779   for (Standard_Integer ij = 0; ij <= Nb; ij++) {
1780     U = (Nb-ij)*U1 + ij*U2;
1781     U /= Nb;
1782     DummyPoles(ij+1) = F.Value(U);
1783     DummyKnots(ij+1) = ij;
1784   }
1785   Handle(Geom2d_BSplineCurve) DummyC2d =
1786     new Geom2d_BSplineCurve(DummyPoles, DummyKnots, DummyMults, 1);
1787 #ifdef DRAW
1788   Standard_CString Temp = "bs2d";
1789   DrawTrSurf::Set(Temp,DummyC2d);
1790 #endif
1791 //  DrawTrSurf::Set((Standard_CString ) "bs2d",DummyC2d);
1792   Handle(Geom2dAdaptor_HCurve) DDD =
1793     Handle(Geom2dAdaptor_HCurve)::DownCast(InitCurve2d);
1794
1795 #ifdef DRAW
1796   Temp = "initc2d";
1797   DrawTrSurf::Set(Temp,DDD->ChangeCurve2d().Curve());
1798 #endif
1799 //  DrawTrSurf::Set((Standard_CString ) "initc2d",DDD->ChangeCurve2d().Curve());
1800 #endif
1801
1802   Standard_Integer Deg1,Deg2;
1803 //  Deg1 = 8;
1804 //  Deg2 = 8;
1805   Deg1 = 2; //IFV
1806   Deg2 = 8; //IFV
1807
1808   Approx_FitAndDivide2d Fit(F,Deg1,Deg2,Tol3d,Tol2d, //OCC217
1809   //Approx_FitAndDivide2d Fit(F,Deg1,Deg2,myTolerance,myTolerance,
1810                             Standard_True);
1811
1812   if(Fit.IsAllApproximated()) {
1813     Standard_Integer j;
1814     Standard_Integer NbCurves = Fit.NbMultiCurves();
1815     Standard_Integer MaxDeg = 0;
1816     // To transform the MultiCurve into BSpline, it is required that all
1817     // Bezier constituing it have the same degree -> Calculation of MaxDeg
1818     Standard_Integer NbPoles  = 1;
1819     for (j = 1; j <= NbCurves; j++) {
1820       Standard_Integer Deg = Fit.Value(j).Degree();
1821       MaxDeg = Max ( MaxDeg, Deg);
1822     }
1823
1824     NbPoles = MaxDeg * NbCurves + 1;               //Tops on the BSpline
1825     TColgp_Array1OfPnt2d  Poles( 1, NbPoles);
1826
1827     TColgp_Array1OfPnt2d TempPoles( 1, MaxDeg + 1);//to augment the degree
1828
1829     TColStd_Array1OfReal Knots( 1, NbCurves + 1);  //Nodes of the BSpline
1830
1831     Standard_Integer Compt = 1;
1832     for (i = 1; i <= NbCurves; i++) {
1833       Fit.Parameters(i, Knots(i), Knots(i+1));
1834       AppParCurves_MultiCurve MC = Fit.Value( i);       //Load the Ith Curve
1835       TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Retrieve the tops
1836       MC.Curve(1, Poles2d);
1837
1838       //Eventual augmentation of the degree
1839       Standard_Integer Inc = MaxDeg - MC.Degree();
1840       if ( Inc > 0) {
1841 //      BSplCLib::IncreaseDegree( Inc, Poles2d, PLib::NoWeights(),
1842         BSplCLib::IncreaseDegree( MaxDeg, Poles2d, PLib::NoWeights(),
1843                          TempPoles, PLib::NoWeights());
1844         //update of tops of the PCurve
1845         for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1846           Poles.SetValue( Compt, TempPoles( k));
1847           Compt++;
1848         }
1849       }
1850       else {
1851         //update of tops of the PCurve
1852         for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1853           Poles.SetValue( Compt, Poles2d( k));
1854           Compt++;
1855         }
1856       }
1857
1858       Compt--;
1859     }
1860
1861     //update of fields of  ProjLib_Approx
1862     Standard_Integer NbKnots = NbCurves + 1;
1863
1864     TColStd_Array1OfInteger   Mults( 1, NbKnots);
1865     Mults.Init(MaxDeg);
1866     Mults.SetValue( 1, MaxDeg + 1);
1867     Mults.SetValue(NbKnots, MaxDeg + 1);
1868     myProjIsDone = Standard_True;
1869     Handle(Geom2d_BSplineCurve) Dummy =
1870       new Geom2d_BSplineCurve(Poles,Knots,Mults,MaxDeg);
1871
1872     // try to smoother the Curve GeomAbs_C1.
1873
1874     Standard_Boolean OK = Standard_True;
1875
1876     for (Standard_Integer ij = 2; ij < NbKnots; ij++) {
1877       OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,Tol3d);  //OCC217
1878       //OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,myTolerance);
1879     }
1880 #ifdef OCCT_DEBUG
1881     if (!OK) {
1882       cout << "ProjLib_ComputeApproxOnPolarSurface : Smoothing echoue"<<endl;
1883     }
1884 #endif
1885     return Dummy;
1886   }
1887   return Handle(Geom2d_BSplineCurve)();
1888 }
1889
1890 //=======================================================================
1891 //function : BSpline
1892 //purpose  :
1893 //=======================================================================
1894
1895 Handle(Geom2d_BSplineCurve)
1896      ProjLib_ComputeApproxOnPolarSurface::BSpline() const
1897
1898 {
1899 //  Modified by Sergey KHROMOV - Thu Apr 18 11:16:46 2002 End
1900 //   Standard_NoSuchObject_Raise_if
1901 //     (!myProjIsDone,
1902 //      "ProjLib_ComputeApproxOnPolarSurface:BSpline");
1903 //  Modified by Sergey KHROMOV - Thu Apr 18 11:16:47 2002 End
1904   return myBSpline ;
1905 }
1906
1907 //=======================================================================
1908 //function : Curve2d
1909 //purpose  :
1910 //=======================================================================
1911
1912 Handle(Geom2d_Curve)
1913      ProjLib_ComputeApproxOnPolarSurface::Curve2d() const
1914
1915 {
1916   Standard_NoSuchObject_Raise_if
1917     (!myProjIsDone,
1918      "ProjLib_ComputeApproxOnPolarSurface:2ndCurve2d");
1919   return my2ndCurve ;
1920 }
1921
1922
1923 //=======================================================================
1924 //function : IsDone
1925 //purpose  :
1926 //=======================================================================
1927
1928 Standard_Boolean ProjLib_ComputeApproxOnPolarSurface::IsDone() const
1929
1930 {
1931   return myProjIsDone;
1932 }