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